窄带高斯白噪声

窄带高斯随机信号功率谱的产生

摘要

基于Visual Studio的开发环境，本文首先通过线性同余法产生高斯白噪声。接着利用FIR

数字滤波器的原理，产生低通滤波器。然后让高斯信号过滤波器，最后通过莱斯表达式的方

式表示窄带高斯白噪声，求出其对应的自相关函数，利用维纳—辛钦公式，对自相关函数进

行傅里叶变换便可以得出功率谱密度函数。

一、问题重述、分析

画出窄带高斯白噪声的功率谱。

要画出窄带高斯白噪声的功率谱，最主要的问题是求出它的自相关函数，然后通过维纳

—辛钦公式便可求出。而对于自相关函数的求解主要是得出窄带高斯信号，所以本文基于莱

斯表达式，求出它的窄带高斯信号。具体流程如图（1）：

sinwt

0

Nt

a



低通滤

波器



加法

器

减

高斯白噪声

Nt

b



低通滤

波器

窄带高斯随

机信号



加

coswt

0

图(1)算法流图

二、问题的解决与算法分析

1、产生高斯白噪声

基于前一次大作业，我们知道产生白噪声的方法有线性同余法。

void gauss(float b[])

{

SYSTEMTIME sys;

GetLocalTime(&sys);

int ed;

float a[100];

ed=conds; //随机种子输入,产生范围为0—2^31-1

int i,j,w,k;

float x,y;

for(w=0;w<5000;w++)

for(k=0;k<20;k++) //5000*20个均匀分布的统计近似服从正态分布

{

x=0;y=0;

for(i=0;i<100;i++)

{

ed=16807*(ed%127773); //利用乘同余法产生其他随机数,共100个

if (ed<0)

ed += 2147483647; //处理溢出的情况

a[i]=ed/2147483647.0; //a即为0-1均匀分布数

}

for(j=20;j<32;j++)

{

x=a[j]+x; //取其中12的个数的和，方便计算

}

y=x-6; //产生标准正态分布随机数

for(j=0;j

{

if(-3+6.0/L1*jy)

b[j]=b[j]+1; //统计各个范围内随机数个数

}

}

int gdriver=DETECT,gmode,errorcode;

initgraph(&gdriver,&gmode,"");

initgraph(1000, 600);

int X0=100,Y0=500;

tcolor(RED);

for(i=0;i

{

moveto(X0+i/2,Y0);

lineto(X0+i/2,Y0-b[i]);

}

getch();

clograph();

}

2、低通滤波器的产生

对于低通滤波器的选取我们这里选取FIR（有限脉冲响应滤波器）滤波器。我们这里采用

窗化法的方法。对于窗化法我们首先要先确定截止频率，然后通过阻带最小衰减的分贝

w

d

值确定窗函数的类型，依据主瓣宽度小于过渡带宽度的指标，确定点数N，然后求出M（位

移量）。最后确定窗函数。

窗函数的产生在数字信号处理的实验课中已经做过，所以此处相对比较简单。

代码：

//*************低通滤波器

int firlpf(float h[],float as,float wp,float ws,float fs)//h[]为滤波器，as为 wp ws 单位是Hz

{

float wd,bt;

int m,n,i;

wp=wp/fs*2*PI;

ws=ws/fs*2*PI;

wd=(wp+ws)/2;

bt=ws-wp;

if(as<21)//选择矩形窗

{

n=ceil(4*PI/bt);

3、高斯信号过低通滤波器

高斯信号过低通滤波器之后变成窄带高斯信号，但是是低频的。然后通过莱斯表达式

产生窄带高斯噪声（高频）----目标信号。

根据莱斯表达式的方式，我们将两个高斯信号通过低通滤波器后分别乘以正弦和余弦函

数，将信号搬移到高频区域。最后求出信号的自相关函数，对自相关函数做傅里叶变换得

出功率谱密度函数。

(1)利用公式可求出窄带高斯信号的自相关函数Rx（m）；

1

R[m]xx

根据课本第180页自相关函数估计公式，，求出自相关函数后利

N

用FFT可求出其功率谱密度函数；

代码：

for(m=0;m<L1+NN-1;m++) //求出自相关函数

{

for(i=0;i<L1+NN-1-m;i++)

{

r[m]=r[m]+N[i]*N[i+m];

}

r[m]=r[m]/(L1+NN-1);

}

N|m|1



i0

ii|m|

（2）因为自相关函数Rx（m）与功率谱密度函数Gx（w）互为傅里叶变换对，所以可以利用

FFT求出功率谱密度函数；

fft(r,xi,1,L1+NN-1);//维纳—辛钦公式产生功率谱密度

//FFT的实现程序

void fft(float xr[],float xi[],int m,int n)

{

int loop1,loop2,loop3,le,le1,ip,i=1;

float wr,wi,tr,ur1,ur,ui,ls,ti;

ls=log((float)(n+1))/log(2.0);

ss(xr,xi,n);

for(loop1=1;loop1<=ls;loop1++)

{

ur=1.0;

ui=0.0;

le=pow(2.0,loop1);

le1=le/2;

wr=cos(PI*m/le1);

wi=-sin(PI*m/le1);

for(loop2=0;loop2

{

for(loop3=loop2;loop3

{

ip=loop3+le1;

tr=ur*xr[ip]-ui*xi[ip];

ti=ur*xi[ip]+ui*xr[ip];

xr[ip]=xr[loop3]-tr;

xi[ip]=xi[loop3]-ti;

xr[loop3]=tr+xr[loop3];

xi[loop3]=ti+xi[loop3];

}

ur1=wr*ur-ui*wi;

ui=wr*ui+ur*wi;

ur=ur1;

}

}

if(m==-1)

{

for(i=0;i

{

xr[i]=xr[i]/n;

xi[i]=xi[i]/n;

}

}

}

void ss(float xr[],float xi[],int n) //比特逆序

{

int i=0,j,s;

float a,bj;

for(j=1;j

{

for(s=n/2;s<=i;s=s/2)

{

i=i-s;

}

i=i+s;

if(i>j)

N2(t)的时域图

窄带高斯噪声的功率谱密度函数

结果分析：通过上面的图形我们可以看出我们产生的高斯白噪声与标准正态分布较为接

近，窗函数的时域和频域的图形都较为接近理论。如果想要进一步提高精度可以

1、增加高斯白噪声的准确度。

2、减小窗函数的过渡带。

3、若要提高分辨率，则提高fft点数。

三、心得体会

通过本次大作业，对知识点的理解，还有编程能力的提高都有很大帮助。比如对于随机

信号产生，高斯信号的产生，还有如何去处理时域与频域的关系，如何将这个问题转换到

另一个问题，通过这种联想的思想去解决生活中的问题。但是，从宏观角度来看，我最大