求大神帮忙看看，程序为什么不能调用函数？程序有点长。。。
主要是bp_decode函数的调用，如果ROW 64，COL 128程序是可以运行的，但是把他们的值变大，程序就出错，bp_decode函数都不能调用，求大神~错误说是stack overflow,不知道怎么改#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#define ROW 504
#define COL 1008
#define pi 3.1415926
#define times 10 //码流组数
#define LOOP 50 //每个DB值循环次数
//-------------------------------------------------------------------
//-------------------------------------------------------------------
typedef unsigned char UINT8;
typedef unsigned int UINT16;
typedef unsigned int UINT32;
typedef char INT8;
typedef int INT16;
typedef int INT32;
//-------------------------------------------------------------------
void genH(INT16 row,INT16 col);
//-------------------------------------------------------------------
void showH();
//-------------------------------------------------------------------
void bpsk(bool input[],INT16 amp,INT16 tx_waveform[]);
//-------------------------------------------------------------------
void awgn(INT16 tx_waveform[],double rx_waveform[],INT16 length,double SNR);
//-------------------------------------------------------------------
void wgn(INT16 length,double noisepower,double noise_s[]);
//-------------------------------------------------------------------
void randn(double noise[],INT16 length);
//-------------------------------------------------------------------
void bp_decode(double rx_waveform[],double SNR_db,INT16 amp,INT8 H[][COL],INT16 rearranged_cols[],bool output[]);
//-------------------------------------------------------------------
double atanh(double x);
//-------------------------------------------------------------------
INT16 sgn(double x);
struct original_bp
{
float qmn0;
float qmn1;
float dqmn;
float rmn0;
float rmn1;
float qn0;
float qn1;
float alpha;
};
//-------------------------------------------------------------------
static INT8 H[ROW][COL]={0};
//-------------------------------------------------------------------
static INT16 rearranged_cols[ROW]={0};
//-------------------------------------------------------------------
void main()
{
bool s[times*(COL-ROW)],input[times*COL],output[times*(COL-ROW)],stage_s[COL-ROW],stage_input[COL],stage_output[COL-ROW];
INT16 tx_waveform[times*COL],i,j,loop;
int error_num=0;
double rx_waveform[times*COL]={0},stage_rx[COL]={0};
srand( (unsigned)time( NULL ) );//用即时的系统时间来做随机数种子.生成随机数
clock_t start,period;
INT16 amp;//振幅
double error_p;//误码率
double SNR_db,SNR,dbnum;//SNR:信噪比
//配置信道信息
amp=1;
//循环多次提高平均性
start=clock();
//产生H矩阵
genH(ROW,COL);
period=clock();
printf("产生H矩阵用时%d\n",clock());
//showH();
for (dbnum=0;dbnum<=4;dbnum+=0.1)
{
SNR_db=dbnum;
SNR=pow(10.0,(SNR_db/10));
for (int loop=0;loop<LOOP;loop++)
{
//产生全零码字
for (int i=0;i<times*(COL-ROW);i++)
{
s[i]=bool(0);
}
//进行BPSK调制
bpsk(input,amp,tx_waveform);
//加高斯白噪声
awgn(tx_waveform,rx_waveform,times*COL,SNR);
//解码
for (int i=0;i<times;i++)
{
for (int j=0;j<COL;j++)
{
stage_rx[j]=rx_waveform[i*COL+j];
}
//原始BP译码算法
printf("ok\n");
bp_decode(stage_rx,SNR,amp,H,rearranged_cols,stage_output);
for (int j=0;j<COL-ROW;j++)
{
output[i*(COL-ROW)+j]=stage_output[j];
}
}
for (int i=0;i<times*(COL-ROW);i++)
{
if (s[i]!=output[i])
{
error_num++;
}
}
}//endofloop 总循环
printf("本次db值循环用时：%d\n",clock()-period);
period=clock();
error_p=(1.0*error_num)/(LOOP*times*(COL-ROW));
printf("%d个码，SNR=%3.3fdb,误码个数为%f,误码率为%lf\n",times*(COL-ROW),SNR_db,(1.0*error_num)/LOOP,error_p);
if(error_p==0)
break;
error_num=0;
}//endof db循环
printf("原始码字：\n");
for (i=0;i<times*(COL-ROW);i++)
{
printf("%d ",s[i]);
}
printf("\n");
printf("输出码字：\n");
for (i=0;i<times*(COL-ROW);i++)
{
printf("%d ",output[i]);
}
}
//////////////////////////////////////////////////////////////////////////
//生成校验矩阵H
void genH(INT16 row,INT16 col)
{
INT16 row_flag[ROW]={0};//每一行1的个数
INT16 bit_per_col=3;//列重3
INT16 row_num[3]={0};//一列中1的横坐标
INT16 max_ones_per_row=0;//每一行1的最多个数
INT16 i,j,k,r,newrow,loop,common,thecol,col_rearranged,looptwo;
bool ckfinish,flag;
INT16 ones_position[ROW/2]={0},ones_count=0;//每一行中1的位置及1的个数。
srand( (unsigned)time( NULL ) );//用即时的系统时间来做随机数种子.生成随机数
//每列随机产生3个1
for (j=0;j<COL;j++)
{
for (i=0;i<bit_per_col;i++)
{
row_num[i]=rand()%ROW;
//避免产生的随机行重复
while (i==1&&row_num[i]==row_num[0])
{
row_num[i]=rand()%ROW;
}
while (i==2&&(row_num[i]==row_num[0]||row_num[i]==row_num[1]))
{
row_num[i]=rand()%ROW;
}
H[row_num[i]][j]=1;
row_flag[row_num[i]]++;//记录每一行1的个数
}
}//end of 每列产生3个随机1
max_ones_per_row=ceil(COL*bit_per_col*1.0/ROW);
//在列上分散1的位置，使得每行1的个数相等或相近
loop=10;
for (loop=0;loop<10;loop++)
{
flag=1;
for (i=0;i<ROW;i++)
{
while (row_flag[i]>max_ones_per_row)
{
flag=0;//有一行中的1大于最大允许
j=rand()%COL;
if (H[i][j]==1)//随机选择该行上某一为1的列，将该列该行上的这个1分散到其他行
{
newrow=rand()%ROW;//随机查找新的行
k=0;
while ((row_flag[newrow]>=max_ones_per_row||H[newrow][j]==1)&&k<ROW)
{
newrow=rand()%ROW;
k++;
}
if (H[newrow][j]==0)
{
H[newrow][j]=1;
row_flag[newrow]++;
H[i][j]=0;
row_flag[i]--;
}
}//end of if H[][]==1
}//end of while
}//end of for i=0:row-1
if (flag==1)
break;
}//end of loop
//去除4环性
loop=10;
for (k=0;k<loop;k++)
{
ckfinish=1;//bool型变量
for (r=0;r<ROW;r++)
{
//记录每一行1的位置
ones_count=0;
for (j=0;j<COL;j++)
if (H[r][j]==1)
{
ones_position[ones_count]=j;
ones_count++;
}
for (i=0;i<ROW;i++)//遍历其他行
{
common=0;//
if (i!=r)//与除该行外的其他行进行比较
{
for (j=0;j<ones_count;j++)
{
if (H[i][ones_position[j]]==1)
{
common++;
if (common==1)
{
thecol=ones_position[j];//有重叠1的列
}
}//endof if H[][]==1
if (common==2)
{
ckfinish=0;//检测到存在4环现象，循环不结束，继续进行
common--;
if (rand()%2==0)//随机决定保留前面列还是后面列
{
col_rearranged=thecol;//保留后面的列，
thecol=ones_position[j];
}
else
col_rearranged=ones_position[j];//保留前面的列，交换后面的列
/************有待考证*******************************************/
//被交换列的1置为3
H[i][col_rearranged]=3;
newrow=rand()%ROW;
looptwo=0;
while (H[newrow][col_rearranged]!=0 && looptwo<5)//只搜索为0的
{
newrow=rand()%ROW;
looptwo++;
}
if (looptwo>=5)//超过5次后搜索范围扩大到所有0和3
{
while (H[newrow][col_rearranged]==1)
{
newrow=rand()%ROW;
}
}
H[newrow][col_rearranged]=1;
}//endof if common==2
/*******************************************************/
}//endof for j=0:ones_count-1
}//endof if r!=i
}//遍历其他行
}//endof for i=0:ROW-1
if (ckfinish==1)//如果本次循环已经不存在四环，则循环结束
{
printf("breakloop=%d\n",k);
break;
}
}//end of loop
//将所有的3变为0
for (i=0;i<ROW;i++)
for(j=0;j<COL;j++)
{
if (H[i][j]==3)
{
H[i][j]=0;
}
}
//H矩阵生成成功
}
//////////////////////////////////////////////////////////////////////////
//显示H矩阵
void showH()
{
INT16 num_ones=0;
printf("H矩阵：\n");
for (INT16 i=0;i<ROW;i++)
{
for (int j=0;j<COL;j++)
{
printf("%d",H[i][j]);
if (H[i][j]==1)
{
num_ones++;
}
}
printf("\n");
}
printf("%d\n",num_ones);
}
//////////////////////////////////////////////////////////////////////////
//用BPSK进行调制
void bpsk(bool input[],INT16 amp,INT16 tx_waveform[])
{
INT16 i;
for (i=0;i<times*COL;i++)
{
if (input[i]==1)
tx_waveform[i]=amp;
else
tx_waveform[i]=-amp;
}
}
//////////////////////////////////////////////////////////////////////////
//调制后信号通过AWGN
void awgn(INT16 tx_waveform[],double rx_waveform[],INT16 length,double SNR)
{
INT16 i,sum=0;
double sigpower,noisepower,noise_s[times*COL];
for (i=0;i<length;i++)
{
sum=sum+tx_waveform[i]*tx_waveform[i];
}
sigpower=(1.0*sum)/length;
//计算所需的噪声功率,这里用的SNR单位不是db
noisepower=sigpower/SNR;
//产生高斯白噪声
wgn(length,noisepower,noise_s);
//加噪声
for (i=0;i<length;i++)
{
rx_waveform[i]=tx_waveform[i]+noise_s[i];
}
}
//////////////////////////////////////////////////////////////////////////
//产生高斯白噪声
void wgn(INT16 length,double noisepower,double noise_s[])
{
INT16 i;
double random_s[times*COL],imp;
//产生高斯分布序列
randn(random_s,length);
imp=sqrt(1.0*noisepower);
for (i=0;i<length;i++)
{
noise_s[i]=imp*random_s[i];
}
}
//////////////////////////////////////////////////////////////////////////
//产生高斯分布序列
void randn(double random_s[],INT16 length)
{
double x1[times*COL],x2[times*COL];
INT16 i;
srand( (unsigned)time( NULL ) );
for(i=0;i<length;i++)
{
x1[i]=(1.0*rand())/RAND_MAX;
x2[i]=(1.0*rand())/RAND_MAX;
random_s[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*pi);
}
}
//////////////////////////////////////////////////////////////////////////
//用原始BP算法进行解码
void bp_decode(double rx_waveform[],double SNR,INT16 amp,INT8 H[][COL],INT16 rearranged_cols[],bool output[])
{
printf("ok\n");
struct original_bp newh[ROW][COL]={0};
float p11[COL],p10[COL],drmn,prod_rmn0,prod_rmn1,const1,const2,const3,const4,alpha_n;
INT16 i,j,k,r,iteration,sum;
INT16 ones_position[COL]={0},ones_count=0;
bool mid_out[COL],Iszero,temp;
//先验概率
for (i=0;i<COL;i++)
{
p11[i]=1/(1+exp(-2*amp*rx_waveform[i]*2*SNR));//N0/2=方差
p10[i]=1-p11[i];
}
//初始化：对特定信道预设信息比特的先验概率
for (i=0;i<ROW;i++)
for (j=0;j<COL;j++)
{
if (H[i][j]==1)
{
newh[i][j].qmn0=p10[j];
newh[i][j].qmn1=p11[j];
newh[i][j].alpha=1.0;
}
}
//showH();
//迭代100次
iteration=100;
for (r=0;r<iteration;r++)
{
//横向步骤：置信概率由信息节点传输到校验节点
for (i=0;i<ROW;i++)//计算概率差
for (j=0;j<COL;j++)
{
if (H[i][j]==1)
{
newh[i][j].dqmn=newh[i][j].qmn0-newh[i][j].qmn1;
}
}
for (i=0;i<ROW;i++)
{
ones_count=0;
for (j=0;j<COL;j++)//记录该行 为1的列数
{
if (H[i][j]==1)
{
ones_position[ones_count]=j;
ones_count++;
}
}
for (j=0;j<ones_count;j++)
{
drmn=1.0;
for (k=0;k<ones_count;k++)
{
if (k!=j)
{
drmn=drmn*newh[i][ones_position[k]].dqmn;
}
}
newh[i][ones_position[j]].rmn0=(1+drmn)/2;
newh[i][ones_position[j]].rmn1=(1-drmn)/2;
//////////////////////////////////////////////////////////////////////////
}
}
//纵向步骤:由校验节点向信息节点传输信息
for (j=0;j<COL;j++)
{
ones_count=0;
for (i=0;i<ROW;i++)//记录该列 为1的行数
{
if (H[i][j]==1)
{
ones_position[ones_count]=i;
ones_count++;
}
}
for (i=0;i<ones_count;i++)
{
prod_rmn0=1.0;
prod_rmn1=1.0;
for (k=0;k<ones_count;k++)
{
if (k!=i)
{
prod_rmn0=prod_rmn0*newh[ones_position[k]][j].rmn0;
prod_rmn1=prod_rmn1*newh[ones_position[k]][j].rmn1;
}
}
const1=p10[j]*prod_rmn0;
const2=p11[j]*prod_rmn1;
newh[ones_position[i]][j].alpha=1.0/(const1+const2);
newh[ones_position[i]][j].qmn0=newh[ones_position[i]][j].alpha*const1;
newh[ones_position[i]][j].qmn1=newh[ones_position[i]][j].alpha*const2;
//更新伪后验概率，因为每列都是一样的qn0qn1,所以只计算每列最后一个就可以.
//经过上面的纵向循环,i-1为ones_position的最后一项
const3=const1*newh[ones_position[i]][j].rmn0;
const4=const2*newh[ones_position[i]][j].rmn1;
alpha_n=1.0/(const3+const4);
newh[ones_position[i]][j].qn0=alpha_n*const3;
newh[ones_position[i]][j].qn1=alpha_n*const4;
//译码尝试，对信息节点的后验概率作硬判决
if (newh[ones_position[i]][j].qn1>0.5)
mid_out[j]=1;
else
mid_out[j]=0;
}
}//endof 纵向
//如果判决条件满足，则停止译码
for (i=0;i<ROW;i++)
{
Iszero=1;
sum=0;
for (j=0;j<COL;j++)
{
sum=sum+H[i][j]*mid_out[j];
}
if (sum%2!=0)
{
Iszero=0;
}
if (Iszero==0)//若有一行作乘后不为0，则判决不满足
{
// printf("判决不满足\n");
break;
}
}
if (Iszero==1)
{
// printf("迭代%d次结束\n",r);
break;
}
}//迭代结束
//获得译码后的输出
for (i=0;i<ROW;i++)
{
if (rearranged_cols[i]!=0)
{
temp=mid_out[i];
mid_out[i]=mid_out[rearranged_cols[i]];
mid_out[rearranged_cols[i]]=temp;
}
}
for (i=ROW;i<COL;i++)
{
output[i-ROW]=mid_out[i];
}
}
//////////////////////////////////////////////////////////////////////////
//sgn函数
INT16 sgn(double x)
{
if (x>0)
return 1;
if (x<0)
return -1;
if (x==0)
return 0;
}
[ 本帖最后由 aialors 于 2013-4-30 21:51 编辑 ]
