今天编了一个用C语言实现高斯白噪声的程序,一方面自己的项目要用,另一方面可以开源的,希望可以给大家带来方便。请大家使用时也注明是baiguangyu001 白栎旸写的,这样我就满足了,呵呵。
做出来的序列*n,其方差比预先规定的要小10的-2次方量级,但是总体来说是比较准确的。可以通过改变量化精度Stride来增加概率估计的准确性。
/* 本程序用于产生SNR=db(dB)的均值为0的白噪声序列,编写人:白栎旸,baiguangyu001 *//* db是信噪比, length是n的长度,n是产生出的白噪声 */# include # include # define PI (3.14)# define Stride (0.005)
short interleaver[length] = {...........};注:interleaver是个交织器,长度与参数length一致就行,里面是乱序排列,这个在matlab里用randperm实现就行。例如:short interleaver[10] = {7,1,9,6,3,2,5,8,4,0};
void awgn(float db, int length, float *n){ float pdf[10000] = {0}; float *amp; float delta=0, amplitude = 0, sum = 0, temp1=0, temp2=0; int *num; int counter, counter1, counter2, counter3; int half_length;
delta = 1 / ( 2 * pow(10,0.1*db) ); // 计算噪声功率 counter = 0; // 通过概率密度计算幅度集合和对应的概率密度 while( sum < (0.5/Stride) ) { pdf[counter] = 1/sqrt(2*PI*delta)*exp(-0.5*pow(amplitude,2)/delta); sum += pdf[counter]; amplitude += Stride; counter++; } // 最后得到的counter是幅度集合的长度 num = (int *)calloc(counter,sizeof(int)); // num每个幅度在序列中出现的次数,幅度放在amp中 amp = (float *)calloc(counter,sizeof(float)); amplitude = 0; sum = 0; temp1 = length*Stride; temp2 = 0; for(counter1=0;counter1<counter;counter1++) { temp2 = temp1*pdf[counter1]; if( (int)(temp2+0.5)>(int)temp2 ) *(num+counter1) = (int)(temp2+1); else *(num+counter1) = (int)temp2; sum += *(num+counter1); *(amp+counter1)=amplitude; amplitude += Stride; } temp1 = 0; half_length = length>>1; if(sum<half_length) // 使num的和恰好是length的一半 { temp1 = half_length-sum; for(counter1=0;counter1<temp1;counter1++) (*(num+counter-1-counter1))++; } else if(sum>half_length) { temp1 = sum-half_length; *num -= temp1; } counter3 = 0; for(counter1=0;counter1<counter;counter1++) // 往*n里填入幅度值 { for(counter2=0;counter2<*(num+counter1);counter2++) { *(n+counter3) = *(amp+counter1); *(n+length-1-counter3) = -*(amp+counter1); counter3++; } }
for(counter1=0;counter1<length;counter1++) { temp1 = *(n+counter1); *(n+counter1) = *(n+interleaver[counter1]); *(n+interleaver[counter1]) = temp1; for(counter2=0;counter2<length-counter1-1;counter2++) if(interleaver[counter1+counter2+1]==counter1) interleaver[counter1+counter2+1] = interleaver[counter1]; }}