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

实验三随机信号的功率谱估计方法报告人：####报告时间：2013.1.2一、实验目的1．利用自相关函数法和周期图法实现对随机信号的功率谱估计。2．观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。3．学习使用FFT提高谱估计的运算速度。4．体会非参数化功率谱估计方法的优缺点。二、实验原理与方法假设信号x(n)为平稳随机过程，其自相关序列定义为*()[()()]mExnxnm(3-1)其中E表示取数学期望，*表示共轭运算。根据定义，x(n)的功率谱密度()P与自相关序列()m存在下面关系：()()jmmPme(3-2)1()()2jmmPed(3-3)然而，实际中我们很难得到准确的自相关序列()m，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。目前，常用的线性谱估计方法有两种，即自相关函数法和周期图法。1．自相关函数法假设我们已知随机信号x(n)的M长的自相关序列{()m}，利用自相关函数法可以得到x(n)的功率谱估计：*11()()()LmimxiximLm11ˆˆ()()MjmmMPme(3-4)利用窗函数，上式又可表达为ˆˆ()()()RjmMmPWmme(3-5)其中，()RMWm为矩形窗函数，定义为1()0RMmMWmmM(3-6)因此，ˆ()P实际上是：真正功率谱()P与窗函数()RMWm傅立叶变换的卷积。矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣。为了降低旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming窗，它定义为0.540.46cos()0HMmMmWmMmM(3-7)这种窗函数可以有效的抑制旁瓣，但是这种方法使主瓣宽度增大，从而降低了谱估计的分辨率，这种主瓣大小和旁瓣干扰之间的矛盾在线性谱估计方法中是无法解决的。2.周期图方法假设已知随机信号()xn的N个样本，利用周期图方法，信号x(n)的功率谱估计为2101ˆ()()NjnnPxneN(3-8)利用上述方法得到的谱估计方差与信号()xn的功率谱平方2()P成正比，谱估计的方差如下：])sinsin(1[)()]([22wNNwwPwIVarxxN（3-9）为了减小它的方差，可以将信号序列进行分段处理，然后再求各分段结果的平均，这就是平均周期图方法，即Bartlett方法。（1）Bartlett平均周期图方法将一个随机序列()xn(0≤n≤N）分成K段，每段长度为L，各段之间互不重迭，因而N=KL，可以想到，第i段的信号序列可表示为()((1))ixnxniL，0,1,...,nLiK(3-10)对于每一段的周期图又可写成210ˆ()()LjniinIxne,1,...,iK(3-11)于是，功率谱估计定义为11ˆˆ()()KiiPIK(3-12)因此，对于固定的记录长度来讲，分段数K增大可使谱估计的方差减小，但是由于L的减小，相应的功率谱主瓣增宽，谱分辨率降低，显然，方差和分辨率也是矛盾的。除了分辨率降低以外，分段处理还会引起序列的长度有限所带来的旁瓣效应。为减小这种影响，最有效的办法是给分段序列用适当的窗函数加权，可以得到较平滑的谱估计，当然，相应的分辨率也有所下降。（2）平滑平均周期图方法这时一种改进的Bartlett周期图方法，它特别适用于FFT直接计算功率谱估值。将长度为N的平稳随机信号序列x(n)分成K段，每段长度为L，即L=N/K。但这里在计算周期图之前，先用窗函数()LWn给每段序列()ixn加权，K个修正的周期图定义为2101()()()LwjniiLnIxnWneLU,1,...,iK(3-13)其中U表示窗函数序列()LWn的能量，201()LLnUWnL(3-14)在这种情况下，功率谱估计可按下面表达式给出：11ˆ()()KwiiPIK(3-15)本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计：()1()()iiNSjniixnaewn(3-16)其中NS为正弦个数，i，i和ia分别为第i个正弦信号的数字频率、相位和幅度，i随机的分布在（0，2）之间，w(n)为零均值方差等于2w的复高斯白噪声。三、实验内容和分析1．仔细阅读有关线性谱估计的内容，根据给出的框图编制自相关函数法谱估计的程序。运行程序，输入,0,1,0,6.0,10,1002111waMN选择矩形窗。观察谱峰位置是否正确（注意：由于窗效应可能引起谱估计的非正定）。答：a.程序流图如3.1所示：（程序代码见附录一）图表3.1实验流程框图图3.1程序流程图b.根据代码运行结果如图3.2：开始输入参数：数据长度N，自相关函数个数M,平均次数K信号产生：输入正弦个数Ns，每个正弦信号的数字频率、相位和幅度，白噪声信号的方差，按照公式（3-16）产生复正弦加白噪声信号的N个采样自相关函数法由N个x(n)估计出自相关序列（M长），并对此自相关序列加矩形窗或Hamming窗，利用公式（3-5）计算[0,2*pi]之间的128个功率谱抽样点周期图方法输入FFT点数NF，按照公式（3-11）、（3-12）、（3-13）计算NF点的功率谱。结束图3.2矩形窗进行功率谱估计的结果分析：由图可知，功率谱的峰值在0.6处，与理论计算值相同。2．观察并记录参数变化对谱估计性能的影响。（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，所以主瓣宽度变宽。同样，旁瓣宽度也明显变宽。图3.3M=5时的矩形窗功率谱估计图（2）.当M=10其他条件不变，用Hanming窗进行功率谱估计的结果如图3.4所示：图3.4M=10时加hanming窗的结果图分析：很明显，经过汉明窗处理后，主瓣宽度变大，旁瓣被抑制，但是该功率谱的分辨率比用矩形窗的结果明显降低。(3).4/1，其它输入同步骤1，仿真结果如图3.5所示图3.5初始相位为45度时的功率谱估计分析：谱估计的峰值位置仍处于0.6处，并且主瓣，旁瓣宽度没有明显变化，说明相位变化对谱估计的性能没有影响。（4）.改变12w，其它输入同步骤1，仿真结果如图3.6所示图3.6噪声方差为1时的功率谱估计分析：高斯白噪声的方差为1时，对谱估计的峰值没有明显影响，但是影响了旁瓣的宽度和幅值。所以，改变信噪比对频谱的影响较小。（5）改变N=10，12w，其它输入同步骤1，其仿真结果如图3.7所示图3.7N=10，噪声方差=1，时的功率谱估计图分析：结合图3.2和图3.7可以得出，当N为100时，改变信噪比对谱估计影响很小；当N较小时（实验中N=10），改变信噪比对谱估计的影响较大。图3.7中，频谱峰值偏离了6.0，同时旁瓣的幅值增大了许多。3．运行程序，输入1,0,8.0,6.0,2,100212121aaNSN，02w，选择矩形窗，调整自相关序列长度M，使得两个正弦频率分量临界分辨出来，纪录此时的M值，并绘制此时的功率谱图。同样，在加Hamming窗的情况下，记录使两个正弦频率分量临界分开的M值，并绘制此时的功率谱图。答：运行程序二，输入1,0,8.0,6.0,2,100212121aaNSN，02w，加矩形窗，调整M，其仿真结果如图3.8,3.9所示；加Hamming窗时分别如图3.10,3.11所示：图3.8M=7时的俩正弦信号功率谱估计图图3.9M=8时的俩正弦信号功率谱估计图分析：由图观察得知，加矩形窗时，当M=8可将两个正弦频率分量临界分辨出来。当M=10，30时，利用汉明窗，得到的结果如图3.10,3.11所示：图3.10M=9时汉明窗进行功率谱估计结果图图3.11M=10时汉明窗进行功率谱估计结果图分析：当M=10时，可将俩个不同频率分量的信号分离。4．运行自相关函数法谱估计程序，输入2100,10,0,1wNMNS，选择矩形窗，观察利用自相关函数法得到的白噪声信号谱估计。改变M=3，20，观察M的变化对白噪声谱估计的影响。答：当M=3，10，20时结果如图3.12,3.13，3.14所示；图3.12M=3时，矩形窗的功率谱估计图3.13M=10时，矩形窗的功率谱估计图3.14M=20时，矩形窗的功率谱估计分析：图3.12,3.13，3.14对比发现：M越大，谱估计的各个分量变得越清晰，也就越容易分辨了。5.根据框图，编制周期图法谱估计程序。自程序FFT可以直接调用Matlab中的函数。运行程序，输入100,128,10NNFK。选择窗函数为矩形窗1NS，21110.6,0,0,0wwa，观察谱峰位置是否正确，并与步骤1结果比较。答：根据绘制的周期图法估计程序（见附录三），得结果如图3.15：图3.15加窗周期图法谱估计图分析：对比步骤1所得的图3.2，图3.15中的谱峰位置同样也出现在0.6pi处，但主瓣幅值明显降低，主瓣宽度增加，另外旁瓣数量减少。6．利用周期图方法重复步骤2、3、4的内容，这里L=N/K相当于自相关函数中的M，观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。答：利用周期图法重复2，3，4的内容，（程序见附录四），这里L=N/K相当于自相关函数法中的M，观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。N/K=10如图3.16：N/K=3如图3.17：N/K=20如图3.18：图3.16N/K=10时加窗的功率谱估计图3.17N/K=3,时加窗的功率谱估计图3.18N/K=20,时加窗的功率谱估计四、思考题1.证明式（1.5）可以表示为：-1-=0ˆˆˆ()=2Real[(m)e]-(0),MjmmP其中Real[]表示取实部。证明：101110ˆˆˆˆ()()()()MMjmjmjmmMmMmPmememek=-m由得：-11k00ˆˆˆ()(-)()MMjkjmmPkemeˆ-(0)由自相关函数的对称性得：-11k00ˆˆˆ()()()MMjkjmmPkemeˆ-(0)即是：10ˆˆ()()(+)MjkjmmPmeeˆ-(0)即得：-1-=0ˆˆˆ()=2Real[(m)e]-(0)MjmmP2.已知实信号(n)x的N个样本(0),(1),...,(1),xxxN，可以定义()AP的估计如下：11ˆˆ()(),NjmAmNPme其中||1*01ˆ()()()NmimxiximN试证明：2101ˆ()()NjnAnPxneN。证明：1||11*1011ˆˆ()()()()NmNNjmjmAmNimNPmexiximeN1||1*()011[()]()NmNjijmiimNxieximeN211*0011[()]()()NNjnjnjnnnxnexnexn