数字信号处理II——随机信号的功率谱估计方法一、实验目的1.利用自相关函数法和周期图法实现对随机信号的功率谱估计。2.观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。3.学习使用FFT提高谱估计的运算速度。4.体会非参数化功率谱估计方法的优缺点。二、实验原理与方法假设信号()xn为平稳随机过程,其自相关序列定义为:*()()()mExnxnm(0.1)其中E表示取数学期望,*表示取共轭。根据定义,()xn的功率谱密度()Pw与自相关序列()m存在如下关系:()()jmmPme(0.2)1()()2jmmPed(0.3)然而,实际中我们很难得到准确的自相关序列()m,只能通过随机信号的一段样本序列来估计信号的自相关序列,进而得到信号的功率谱估计。目前常用的线性谱估计方法有两种:自相关函数法和周期图方法,本实验将对这两种方法分别予以讨论。1.自相关函数法假设已知随机信号()xn的N个观测样本,则其自相关序列可以用下式进行估计:||1*01ˆ()()()||1||NmnmxnxnmmNNm(0.4)当仅使用长度为2M-1的自相关序列时,对其进行傅立叶变换即可得到功率谱估计如下:11ˆˆ()()MjmmMPme(0.5)其中M为加窗长度,Re()cMWm为矩形窗函数,定义如下:Re1,||()0,||cMmMWmmM(0.6)因此,ˆ()Pw在一定程度上可以看作是“真正的功率谱()Pw”与窗函数傅立叶变换的卷积。矩形窗函数不仅降低了谱估计的分辨率,而且使谱估计产生了旁瓣,旁瓣效应使那些处于旁瓣附近功率较小的频率分量被淹没掉。为了降低旁瓣影响,可以采用具有较小旁瓣的窗函数,如Hamming窗,它定义为:0.540.46cos,||()0,||HamMmmMWmMmM(0.7)这种窗函数可以有效的抑制旁瓣,但此时主瓣宽度增大,从而降低了谱估计的分辨率。这种主瓣和旁瓣之间的矛盾在非参数化功率谱估计方法中是无法解决的。2.周期图方法假设已知随机信号()xn的N个样本,利用周期图方法,信号()xn的功率谱估计为:1201ˆ()|()|NjmnPxneN(0.8)利用上述方法得到的谱估计方差与信号()xn的功率谱平方2[()]Pw成正比,为了减小它的方差,可以将信号序列进行分段处理,然后再求各分段结果的平均,这就是平均周期图方法,即Bartlett方法,它特别适用于FFT直接计算功率谱估值。(1)Bartlett平均周期图方法将一个随机序列()xn(0)nN分成K段,每段长度为L,各段之间互不重叠,因而NLK。可以想到,第i段的信号序列可表示为:()((1))0n,1ixnxniLLiK(0.9)每一段的周期图又可写成:1201ˆ()|()|1LjmiinIxneiKL(0.10)于是,功率谱估计定义为:11ˆˆ()()KiiPIK(0.11)对于固定的记录长度来讲,分段数K增大可使谱估计的方差减小,但是由于L的减小,相应的功率谱主瓣增宽,谱分辨率降低,显然,方差和分辨率也是矛盾的。除了分辨率降低以外,分段处理还会引起序列的长度有限所带来的旁瓣效应。为减小这种影响,最有效的办法是给分段序列用适当的窗函数加权,可以得到较平滑的谱估计,当然,相应的分辨率也有所下降。(2)平滑平均周期图方法这是一种改进的Bartlett周期图方法,将长度为N的平稳随机信号序列()xn分成K段,每段长度为L,即/LNK。这里在计算周期图之前,先用窗函数LW(n)给每段序列(n)ix加权,K个修正的周期图定义为:2-1-=01ˆ()=(n)W(n)e1LWjmiiLnIxiKLU(0.12)其中U表示窗函数序列每个样本的平均能量,即:2-1=01=W(n)LLnUL(0.13)在这种情况下,功率谱估计可按下面表达式给出:=11ˆˆ()=()KWiiPIK(0.14)本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计:(n+)=1(n)=+(n)siiNjiixaew(0.15)其中Ns为正弦个数,i,i和ia分别为第i个正弦信号的数字频率、相位和幅度,i随机分布在[0,2)之间,(n)w为零均值,方差等于2w的复高斯白噪声。三、实验内容及步骤1、仔细阅读有关功率谱估计的内容,按照图1给出的框图编制自相关函数法谱估计程序。运行程序,输入2=100,=10,=1,=0.6,=0,=0,siiwNMN选择矩形窗。观察谱峰位置是否正确(注意:由于窗效应可能引起谱估计的非正定)。开始参数输入:数据长度N,自相关函数个数M,平均次数K信号产生:输入正弦信号个数及其数字频率、相位、振幅,方差为的白噪声并产生复正弦加白噪声信号N个采样自相关函数法由x(n)的样本估计自相关序列并对此序列加入矩形窗和Hamming窗再求出在内点128个抽样点周期图法输入FFT点数并由公式计算点的功率谱结束2w[0,2)FNFN图1实验流程框图图2M=10时,矩形窗的谱估计从图中可以看出谱峰位置在0.6*pi处,这个位置是正确的。2、观察并记录参数变化对谱估计性能的影响。(1)改变M=5,其它输入同步骤1,观察功率谱估计的主瓣宽度和旁瓣大小随自相关序列长度的变化情况。图3M=5时,矩形窗的谱估计观察M=10和M=5的仿真结果,可以看出M越小,则主瓣宽度越大,分辨率越低,谱峰高度也降低了;旁瓣越少,旁瓣最大幅度也减小;主瓣与旁瓣的幅值比、宽度比基本没有改变。(2)选择窗函数为Hamming窗,其它输入同步骤1,观察不同的窗函数对谱估计性能的影响。(一般来讲,与矩形窗相比,加Hamming窗,主瓣宽度增加一倍,但旁瓣减小许多)。图4M=10时,hamming窗的谱估计与矩形窗相比,选用汉明窗,功率谱的主瓣宽度增加了近一倍,旁瓣也减少了许多,分辨率降低了。(3)改变1=/4,其它输入同步骤1,观察初始相位的变化对谱估计性能的影响。图54/1,矩形窗的谱估计由图可知改变初始相位,对频谱估计的结果没有影响。(4)改变2=1w,其它输入同步骤1,观察信噪比的变化对谱估计性能的影响。图621w,矩形窗的谱估计当噪声增加,信噪比减小时,会使谱估计得到的结果的精确性减小,而且增加的噪声越大,对谱估计影响就越大,(5)改变2=10,=1wN,其它输入同步骤1,结合(4)的内容,观察数据长度及信噪比对谱估计性能的影响。图7N=10,21w,矩形窗的谱估计由图知,在N较大时,改变信噪比对谱估计影响比较小;而当N较小时,改变信噪比对谱估计的影响较大。在上图中我们可以看出此时频谱峰值已偏离了0.6pi,同时旁瓣的幅值增大了许多。3、运行程序,输入2121212=100,=2,=0.6,=0.8,==0,==1,=0,swNNaa选择矩形窗,调整自相关序列长度M,使得两个正弦频率分量临界分辨出来,记录此时M值,并绘出功率谱图。同样,在加Hamming窗的情况下,记录使两个正弦频率分量临界分辨出来的M值,并绘出此时的功率谱图。图8当M=7时,矩形窗的谱估计图9当M=8时,矩形窗的谱估计由上面两图可知,在M=8时,加矩形窗可以将两个正弦频率分量临界分辨出来。图10当M=9时,汉明窗的谱估计图11当M=10时,汉明窗的谱估计由上两图可知,在M=10时,加汉明窗可以将两个正弦频率临界分辨出来。4、运行自相关函数法谱估计程序,输入2=100,=0,=1,M=10,swNN选择矩形窗,观察利用自相关函数法得到白噪声信号谱估计。改变M=3、20,观察M的变化对白噪声谱估计的影响。图12当M=3时,加矩形窗的白噪声谱估计图13当M=10时,加矩形窗的白噪声谱估计图14当M=20时,加矩形窗的白噪声谱估计由上面三幅图知,M越大,加矩形窗的谱估计各个分量变得越清晰,更容易分辨。5.根据图1给出的框图,编制周期图法谱估计程序,子程序FFT可以直接调用Matlab中的函数。运行程序,输入=100,=128,=10,FNNK选择窗函数为矩形窗,输入2111=1,=0.6,=0,=1,=0,swNa观察谱峰位置是否正确,并与步骤1结果比较。图15N/K=10时的周期图法谱估计此图峰值在0.6*pi处,位置正确,并且由对比可知,主瓣峰值明显下降了,而且其宽充增加,旁瓣数目有所减少。6.利用周期图法重复2,3,4的内容,这里=/KLN相当于自相关函数法中的M,观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。(1)重复步骤2(1),选矩形窗可得如下图所示的周期图法谱估计图16N/K=5时,矩形窗的Bartlett周期图谱估计对比每2(1)步中结果可知,用Bartlett周期图谱估计的主瓣在N/K减小时,主瓣宽度变大,分辨率降低,峰值下降,而且旁瓣数目减少,其规律和自相关函数法中减小M时,规律一致。(2)重复步骤2(2),选Hamming窗可得如下图所示的周期图法谱估计图17N/K=5时,汉明窗的Bartlett周期图谱估计对比第2(2)步中结果可知与矩形窗相比,汉明窗知主瓣宽度加倍,旁瓣减小许多。Bartlett周期图谱估计和自相关谱出现的规律相似。(3)重复步骤2(3)得到如下图所示的周期图法谱估计图184/1,N/K=10时矩形窗的Bartlett周期图谱估计对比每2(3)可知,变化对谱估计没有太大影响。(4)重复步骤2(4)得如下图所示的周期图法谱估计图19取12w,N/K=10时矩形窗的Bartlett周期图谱估计对比第2(4)步可知,增大高斯白噪声的方差时,对谱估计的峰值没有影响,但却使旁瓣增大。(5)重复步骤2(5)得如下图所示的周期图法谱估计图20N=50,12w时矩形窗的Bartlett周期图谱估计由图可知,当N减小,同时增大噪声方差时,矩形窗的Bartlett周期图谱估计的主瓣宽度增加。(5)用周期图法重复步骤3,输入100,2,NNS10.6,20.8,12120,1aa,02w,加矩形窗,调整段数K,其仿真结果如下图所示;加Hamming窗时,仿真如下图所示图21K=20时,矩形窗的Bartlett周期图谱估计图22K=25时,矩形窗的Bartlett周期图谱估计由上两幅图可知,加矩形窗时,K=20时可以将两个正弦频率分量临界分辨出来,但峰值位置有所变化。图23K=20时,汉明窗的Bartlett周期图谱估计图24K=25时,汉明窗的Bartlett周期图谱估计同样由上面两图可知,加汉明窗时,当K=20时,可以将两个正弦频率分量临界分辨出来。同样其频率中心也发生了变化。(6)利用周期图方法重复步骤4,加矩形窗,2100,0,1wNNS,当N/K=10、4、20,仿真结果分别如下图所示图25N/K=4时,矩形窗的Bartlett周期图谱估计图26N/K=10时,矩形窗的Bartlett周期图谱估计图27N/K=20时,矩形窗的Bartlett周期图谱估计由上图可知,N/K越大,谱估计的各个分量变得越清晰,更容易分辨。综上可右,对比自相关函数法的谱估计,由步骤5、6的仿真结果,我们发现周期图法谱估计减小了峰值的幅度,加宽主瓣宽度,但是周期图使方差减小,从而得到一致的谱估计。周期图法在分辨率和稳定性方面是优于自相关函数法的谱估计。四、思考题1.证明式(1.5)可以表示为:-1-=0ˆˆˆ()=2Real[(m)e]-(0),MjmmP其中Real[]表示取实部。证明:101110ˆˆˆˆ()()()()MMjmjmjmmMmMmPmememek=-m由得:-11k00ˆˆˆ()(-)()MMjkjmmPkeme