高斯白噪声 C语言实现

    技术2022-05-20  46

    今天编了一个用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]; }}


    最新回复(0)