数字信号处理――谱估计实验

1、实验三 随机信号的功率谱估计方法报告人: # 报告时间:2013.1.2一、实验目的1.利用自相关函数法和周期图法实现对随机信号的功率谱估计。2.观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。3.学习使用FFT 提高谱估计的运算速度。4.体会非参数化功率谱估计方法的优缺点。二、实验原理与方法假设信号x (n 为平稳随机过程,其自相关序列定义为*(m E x n x n m =+ (3-1其中E 表示取数学期望,*表示共轭运算。根据定义,x (n 的功率谱密度(P 与自相关序列(m 存在下面关系: (j m m P m e -=-=

2、(3-21(2j m m P e d -= (3-3然而,实际中我们很难得到准确的自相关序列(m ,只能通过随机信号的一段样本序列来估计信号的自相关序列,进而得到信号的功率谱估计。目前,常用的线性谱估计方法有两种,即自相关函数法和周期图法。1. 自相关函数法假设我们已知随机信号x(n的M 长的自相关序列(m ,利用自相关函数法可以得到x(n的功率谱估计:*11(L m i m x i x i m L m -=+- 11(M j m m M P m e -=-+= (3-4 利用窗函数,上式又可表达为(R j m M m PW m m e -=-= (3-5 其中,(R MW m 为矩形窗函数,

3、定义为 1(0RM m M W m m M <= (3-6因此,(P 实际上是:真正功率谱(P 与窗函数(R MW m 傅立叶变换的卷积。矩形窗函数不仅降低了谱估计的分辨率,而且使谱估计产生了旁瓣。为了降低旁瓣影响,可以采用具有较小旁瓣的窗函数,如Hamming 窗,它定义为0.540.46cos (0H M m M m W m M m M <+= (3-7这种窗函数可以有效的抑制旁瓣,但是这种方法使主瓣宽度增大,从而降低了谱估计的分辨率,这种主瓣大小和旁瓣干扰之间的矛盾在线性谱估计方法中是无法解决的。2.周期图方法假设已知随机信号(x n 的N 个样本,利用周期图方法,信号x(n

4、的功率谱估计为2101(N j n n P x n eN -= (3-8利用上述方法得到的谱估计方差与信号(x n 的功率谱平方2(P 成正比, 谱估计的方差如下: sin sin (1(22wN Nw w P w I Var xx N + (3-9 为了减小它的方差,可以将信号序列进行分段处理,然后再求各分段结果的平均,这就是平均周期图方法,即Bartlett 方法。(1Bartlett 平均周期图方法将一个随机序列(x n (0n N 分成K 段,每段长度为L ,各段之间互不重迭,因而N=KL ,可以想到,第i 段的信号序列可表示为(1i x n x n i L =+- , 0,1,.,n

5、 L i K <= (3-10对于每一段的周期图又可写成210(L j n i in I x n e -=, 1,.,i K = (3-11于是,功率谱估计定义为11(Kii P I K = (3-12 因此,对于固定的记录长度来讲,分段数K 增大可使谱估计的方差减小,但是由于L 的减小,相应的功率谱主瓣增宽,谱分辨率降低,显然,方差和分辨率也是矛盾的。除了分辨率降低以外,分段处理还会引起序列的长度有限所带来的旁瓣效应。为减小这种影响,最有效的办法是给分段序列用适当的窗函数加权,可以得到较平滑的谱估计,当然,相应的分辨率也有所下降。(2平滑平均周期图方法这时一种改进的Bartlett 周

6、期图方法,它特别适用于FFT 直接计算功率谱估值。将长度为N 的平稳随机信号序列x(n分成K 段,每段长度为L ,即L=N/K 。但这里在计算周期图之前,先用窗函数(L W n 给每段序列(i x n 加权,K 个修正的周期图定义为2101(L w j n i i Ln I x n W n e LU -=,1,.,i K = (3-13其中U 表示窗函数序列(L W n 的能量,201(LL n U W n L = (3-14 在这种情况下,功率谱估计可按下面表达式给出:11(K w i i P IK = (3-15本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计

7、:(1(i i NSj n i i x n a e w n +=+ (3-16其中NS 为正弦个数,i ,i 和i a 分别为第i 个正弦信号的数字频率、相位和幅度,i 随机的分布在(0,2之间,w(n为零均值方差等于2w的复高斯白噪声。 三、实验内容和分析1. 仔细阅读有关线性谱估计的内容,根据给出的框图编制自相关函数法谱估计的程序。运行程序,输入,0,1,0,6.0,10,1002111=w a M N 选择矩形窗。观察谱峰位置是否正确(注意:由于窗效应可能引起谱估计的非正定。 答:a.程序流图如3.1所示:( 程序代码见附录一图表3.1 实验流程框图 图3.1程序流程图b.根据代码运行结

8、果如图3.2:开始输入参数:数据长度 N ,自相关函数个数 M,平均次数 K信号产生:输入正弦个数Ns ,每个正弦信号的数字频率、相位和幅度,白噪声信号的方差,按照公式( 3-16 产生复正弦加白噪声信号的N 个采样 自相关函数法由N 个x(n估计出自相关序列(M 长,并对此自相关序列加矩形窗或Hamming 窗,利用公式(3-5 计算0,2*pi 之间的128个功率谱抽样周期图方法 输入FFT 点数N F ,按照公式(3-11、(3-12、(3-13 计算N F 点的功率谱。结束 图3.2 矩形窗进行功率谱估计的结果分析:由图可知,功率谱的峰值在0.6处,与理论计算值相同。2. 观察并记录参

9、数变化对谱估计性能的影响。(1改变M=5,其它输入同步骤1,观察功率谱估计的主瓣宽度和旁瓣大小随自相关序列长度的变化情况。(2选择窗函数为Hamming 窗,其它输入同步骤1,观察不同的窗函数对谱估计性能的影响。(3改变4/1=,其它输入同步骤1,观察初始相位的变化对谱估计性能的影响。(4改变12=w,其它输入同步骤1,观察信噪比变化对谱估计性能的影响。 (5改变N=10,12=w,其它输入同步骤1,结合(4的内容,观察数据长度及信噪比对谱估计性能的影响。答:(1.M=5时(如图3.3,其他条件不变,经观察得知主瓣宽度约在(0.38-0.82,M=10时为0.48-0.70,所以主瓣宽度变宽。

10、同样,旁瓣宽度也明显变宽。 图3.3 M=5时的矩形窗功率谱估计图(2.当M=10其他条件不变,用Hanming 窗进行功率谱估计的结果如图3.4所示: 图3.4 M=10时加hanming 窗的结果图分析:很明显,经过汉明窗处理后,主瓣宽度变大,旁瓣被抑制,但是该功率谱的分辨率比用矩形窗的结果明显降低。(3 . 4/1=,其它输入同步骤1,仿真结果如图3.5所示 图3.5 初始相位为45度时的功率谱估计分析:谱估计的峰值位置仍处于0.6处,并且主瓣,旁瓣宽度没有明显变化,说明相位变化对谱估计的性能没有影响。(4.改变12=w,其它输入同步骤1,仿真结果如图3.6所示 图3.6 噪声方差为1时

11、的功率谱估计分析:高斯白噪声的方差为1时,对谱估计的峰值没有明显影响,但是影响了旁瓣的宽度和幅值。所以,改变信噪比对频谱的影响较小。(5改变N=10,12=w,其它输入同步骤1,其仿真结果如图3.7所示 图3.7 N=10,噪声方差=1,时的功率谱估计图分析:结合图3.2和图3.7可以得出,当N 为100时,改变信噪比对谱估计影响很小;当N 较小时(实验中N=10,改变信噪比对谱估计的影响较大。图3.7中,频谱峰值偏离了6.0,同时旁瓣的幅值增大了许多。3.运行程序,输入1,0,8.0,6.0,2,100212121=a a NS N ,02=w ,选择矩形窗,调整自相关序列长度M ,使得两个

12、正弦频率分量临界分辨出来,纪录此时的M 值,并绘制此时的功率谱图。同样,在加Hamming 窗的情况下,记录使两个正弦频率分量临界分开的M 值,并绘制此时的功率谱图。 答:运行程序二,输入1,0,8.0,6.0,2,100212121=a a NS N ,02=w ,加矩形窗,调整M ,其仿真结果如图3.8,3.9所示;加Hamming 窗时分别如图3.10,3.11所示: 图3.8 M=7时的俩正弦信号功率谱估计图 图3.9 M=8时的俩正弦信号功率谱估计图分析:由图观察得知,加矩形窗时,当M=8可将两个正弦频率分量临界分辨出来。当M=10,30时,利用汉明窗,得到的结果如图3.10,3.1

13、1所示: 图3.10 M=9时 汉明窗进行功率谱估计结果图 图3.11 M=10时 汉明窗进行功率谱估计结果图分析:当M=10时,可将俩个不同频率分量的信号分离。4.运行自相关函数法谱估计程序,输入2100,10,0,1wN M NS =,选择矩形窗,观察利用自相关函数法得到的白噪声信号谱估计。改变M=3,20,观察M 的变化对白噪声谱估计的影响。答:当M=3,10,20时结果如图3.12, 3.13,3.14所示; 图3.12 M=3时,矩形窗的功率谱估计 图3.13 M=10时,矩形窗的功率谱估计 图3.14 M=20时,矩形窗的功率谱估计分析:图3.12, 3.13,3.14对比发现:M

14、 越大,谱估计的各个分量变得越清晰,也就越容易分辨了。5. 根据框图,编制周期图法谱估计程序。自程序FFT 可以直接调用Matlab 中的函数。运行程序,输入100,128,10N NF K =。选择窗函数为矩形窗1NS =,21110.6,0,0,0w w a =,观察谱峰位置是否正确,并与步骤1结果比较。答:根据绘制的周期图法估计程序(见附录三,得结果如图3.15: 图3.15 加窗周期图法谱估计图分析:对比步骤1所得的图3.2,图3.15中的谱峰位置同样也出现在0.6pi处,但主瓣幅值明显降低,主瓣宽度增加,另外旁瓣数量减少。6.利用周期图方法重复步骤2、3、4的内容,这里L=N/K相当

15、于自相关函数中的M,观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。答:利用周期图法重复2,3,4 的内容,(程序见附录四,这里L = N /K相当于自相关函数法中的M ,观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。N/K =10 如图3.16:N/K =3 如图3.17:N/K =20 如图3.18: 图3.16 N/K=10时加窗的功率谱估计 图3.17 N/K=3,时加窗的功率谱估计 图3.18 N/K=20,时加窗的功率谱估计四、思考题1. 证明式(1.5可以表示为:-1-=0(=2Real(me -(0,M j m m P其中Real表示取实部

16、。 证明:10111(M M j m j m j mm M m M m Pm e m e m e-=-+=-+=+k=-m 由得:-11k 00(-(M M j kj m m Pk e m e -=+-(0 由自相关函数的对称性得:-11k 0(M M j k j m m P k e m e -=+-(0即是: 10(+M j kj m m Pm e e -=-(0即得: -1-=0(=2Real(me -(0M j m m P 2. 已知实信号(nx 的 N 个样本(0,(1,.,(1,x x x N -,可以定义(A P 的估计如下:11(,N j m Am N P m e -=-+=其中

17、|1*01(N m i m x i x i m N -=+试证明:211(N j nAn P x n eN-=。证明:1|11*111(N m N N j m j m A m N i m N P m e x i x i m e N -=-+=-+=+=1|1*(11(N m N j i j m i i m N x i e x i m e N-+=-+211*11(N N j n j nj nn n x n ex n e x n eNN-= #附录一.N=100;%样本数从0到N-1M=10;%是矩形窗长度变量,2M-1是矩形窗长度w1=0.6*pi;fai1=0;a1=1;sigma_w2=0

18、;%方差为0Ns=1;I=256;%截取样本为256个%信号产生xn=zeros(N,1;wn=sqrt(sigma_w2*( randn(N,1;%噪声for n=1:Nxn(n=a1*exp(j*(w1*n+fai1+wn(n;%定义序列endr_xx=xcorr(xn,'unbiased'%产生自相关序列long=2*M-1;%矩形窗Wrect=rectwin(long;%矩形窗p1=r_xx(N-M+1:N+M-1.*Wrect;p2=fft(p1,I;p3=abs(p2(1:I/2;%显示128 个功率谱采样点figure(1,clfstem(0:I/2-1/(I/2

19、,p3,'r'grid on; %归一化title('矩形窗功率谱估计'%hanming窗Whanmin=hamming(long;p11=r_xx(N-M+1:N+M-1.*Whanmin;p22=fft(p11,I;p33=abs(p22(1:I/2;figure(2,clfstem(0:I/2-1/(I/2,p33,'b'grid on;title('加hamming 窗功率谱估计'附录二.clc;clear all;n=100;m=20;ns=0;w1=0.6*pi;w2=0.8*pi;fi1=0;fi2=0;a1=1;a

20、2=1;sigma2=1;l=256;sn1=zeros(n,1;sn2=zeros(n,1;xn=zeros(;sn1=a1*exp(i*(w1*(1:n'+fi1;%定义W1的信号函数sn2=a2*exp(i*(w2*(1:n'+fi2;%定义W2的信号函数wn=sqrt(sigma2.*(randn(n,1;%定义噪声xn=sn1+sn2+wn;%俩个不同频率分量的信号与噪声叠加xxn=xcorr(xn,'unbiased' %产生信号long=2*m-1; %矩形窗的定义长度为2M-1; Rectangle=window(rectwin,long;p1=

21、xxn(n-m+1:n+m-1.*Rectangle;p2=fft(p1,l;p3=abs(p2(1:l/2;long1=2*m-1; %汉明窗的定义长度为2M-1; Hammwin=window(hamming,long1;p11=xxn(n-m+1:n+m-1.*Hammwin;p22=fft(p11,l;p33=abs(p22(1:l/2;figure(1;stem(0:l/2-1/(l/2,p3,'r'title('矩形窗功率谱估计'grid on; figure(2;stem(0:l/2-1/(l/2,p33,'.'title('

22、;汉明窗功率谱估计'grid on; 附录三.N=100;Nf=128;K=10;Ns=0;w1=0.6*pi;fai1=0;a1=1;sigma_w2=0;I=N/K;wn=sqrt(sigma_w2*( randn(N,1; x=zeros(N,1;for i=1:Nx(i=a1*exp(j*(w1*i+fai1+wn(i; end%加矩形窗rectangular=window(rectwin,I;u=sum(rectangular.2/I;%计算窗函数序列的能量for i=1:Kxrectwin=x(i-1*I+1:i*I,:.*rectangular; xrectwindft=

23、fft( xrectwin,Nf;ixrec(i,:=xrectwindft.2/u/I;endpwrec=abs(sum(ixrec/K;%加hamming 窗Hamming=window(hamming,I;u=sum(Hamming.2/I;%计算窗函数序列能量for i=1:Kxhamm=x(i-1*I+1:i*I,:.*Hamming;xhammdft=fft(xhamm,Nf;ixhamm(i,:=xhammdft.2/Nf/I;endpwhamm=abs(sum(ixhamm/K;figure(1,clfstem(0:Nf/2-1/(Nf/2-1,pwrec(1:Nf/2,'r'grid on ; title('加矩形窗功率谱估计'figure(2,clfstem(0:Nf/2-1/(Nf/2-1,pwhamm(1:Nf/2,'b'title('加 hamming 窗功率谱估计' 附