功率谱估计性能分析及Matlab仿真

1功率谱估计性能分析及Matlab仿真1引言随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg法等。2经典功率谱估计经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N个样本数据估计出其功率谱[1]。2.1周期图法(Periodogram)Schuster首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。取平稳随机信号()xn的有限个观察值(0),(1),...,(1)xxxn，求出其傅里叶变换10()()NjjnNnXexne然后进行谱估计221()()jNSXeN周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用FFT快速算法来计算。但是，这种方法需要对无限长的平稳随机序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄露，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低[1]。该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;Fn=-0.5:1/N:0.5-1/N;px=fft(x,N);pmax=max(px);%归一化px=px/pmax;px=10*log10(px+0.000001);plot(Fn,fftshift(px));gridon;图1周期图法4096N3图2周期图法128N说明：(1)本报告仿真中所采用的用于功率谱估计的数据文件来自参考文献[3]的test.dat。该数据为128点复序列（图3），由复数噪声加上四个复正弦组成。其归一化频率分别是：12340.15,0.16,0.252,0.16ffff。图3复序列(),[0,127]xnn(2)从仿真图可以清晰看到，1f和2f不能完全分开，仅在波形的顶部能看出4是两个频率分量；此外，当数据长度N太大时（图1），谱曲线呈现较大的起伏；当数据长度N太小时（图2），谱的分辨率又不好。据此，周期图法不满足一致性估计条件。2.2自相关法(BT法)自相关法的理论基础是维纳—辛钦定理。1958年Blackman和Tukey给出了这一方法的具体实现。对于平稳随机信号来说，其自相关函数是确定性函数，故其功率谱也是确定的。这样可由平稳随机离散信号的有限个离散值(0),(1),...,(1)xxxn求出自相关函数101()()()NmxnRmxnxnmN然后在(,)MM内对()xRm做傅里叶变换，得到功率谱()()MjnxmMSRme该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;Fn=-0.5:1/N:0.5-1/N;Mlag=64;rx=xcorr(x,Mlag,'unbiased');px=fft(rx,N);pmax=max(px);%归一化px=px/pmax;px=10*log10(px+0.000001);plot(Fn,fftshift(px));gridon;5图4自相关法不加窗64M图5自相关法不加窗32M6图6自相关法使用汉明窗(Hamming)说明：(1)该方法先由序列()xn估计出自相关函数()xRm，然后对()xRm进行傅里叶变换，便得到()xn的功率谱估计。当延迟与数据长度之比很小时，可以有良好的估计精度。(2)图4是用自相关法（BT法）求出的功率谱，64M没有加窗；图5也是用自相关法（BT法）求出的功率谱，32M，没有加窗；图6同样是采用自相关法求出的功率谱，32M，使用了汉明窗。显然，自相关函数的延迟M越小，谱变得越平滑。2.3Welch法该方法的基本原理是在对随机序列分段时，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均。这样可得功率谱21101()()LMijminSxneMUL其中10()MnUn7这里()n为窗函数。该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;Fn=-0.5:1/N:0.5-1/N;xpsd=pwelch(x,hamming(33),16,N,'whole');mmax=max(xpsd);%归一化xpsd=xpsd/mmax;xpsd=10*log10(xpsd+0.000001);plot(Fn,fftshift(xpsd));gridon;图7Welch法不叠合使用汉明窗(Hamming)8图8Welch法叠合16点使用汉明窗(Hamming)图9Welch法叠合16点使用矩形窗(Boxcar)9图10Welch法叠合16点使用布莱克曼窗(Blackman)说明：(1)因为Welch法允许各段数据交叠，所以数据段数L会增加，使方差得到更大的改善，但是数据的交叠减小了每一段数据的不相关性，使方差的减小不会达到理论程度。另外，采用合适的窗函数可以减少信号的频谱泄露，同时也可以增加谱峰的宽度，从而提高分辨率。(2)图7是利用Welch法求出的周期图，共分四段，每段32点，没有叠合，使用了汉明窗；图8也是利用Welch法求出的周期图，共分四段，每段32点，，使用了汉明窗；图9是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16，且使用了矩形窗；图10是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16，使用了布莱克曼窗。从图中可以看出，由矩形窗处理的谱估计的主瓣宽度最窄，分辨率最好，但是其起伏性较大，所以其方差特性最差。由汉明窗和布莱克曼窗得到的谱估计的主瓣宽度最宽，因此其分辨率相对较差，但其旁瓣较小，大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。因此，选择不同的窗函数其主瓣宽度不一样，造成谱估计的分辨率也不相同。2.4经典功率谱估计的性能比较由以上的Matlab仿真图形和相关结果分析，我们得到了经典谱估计算法性能的直观比较：10(1)周期图法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰。(2)自相关法（BT法）由于使用了平滑窗对周期图法估计的功率谱进行了平滑，因此方差性能较好，功率谱比周期图法估计的要平滑，但其分辨率比周期图法低。(3)Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。综合上述讨论，我们可以对经典谱估计的算法作大致的总结[3]：(1)功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而仍是目前较常用的谱估计方法。(2)谱的分辨率较低，它正比于2/N，N是所使用的数据长度。(3)方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。(4)周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率且增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。3现代功率谱估计由前一章的讨论我们可知，经典功率谱估计方法的方差性能较差，分辨率较低。而现代谱估计技术的目标都是旨在努力改善谱估计的分辨率。参数模型法是现代谱估计的主要内容，参数模型主要分为AR模型、MA模型和ARMA模型。由于AR模型具有一系列好的性能，因此是被研究最多并获得广泛应用的一种模型。本报告中现代功率谱估计的仿真基于的是AR模型。3.1自相关法假定观察到得数据为(0),(1),...,(1)xxxn，而对于无法观察到得区间（即01nnN和），()xn的样本假定为0，观测数据区间之外的数据为0，在均方误差意义下使得数据的预测误差功率最小。由于自相关矩阵ˆpR是Toeplitz矩阵，而且又为正定的，故可利用Levinson-Durbin递归算法高效求解，得到AR模型11参数。该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;fn=-0.5:1/N:0.5-1/N;xpsd=pyulear(x,20,N);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plot(fn,fftshift(xpsd));gridon;图11自相关法10p12图12自相关法20p图13自相关法30p说明：(1)图11、12和13是用自相关法求出的AR谱曲线，阶次p分别等于10，20和30。可以看出，在阶次较低时（图11），分辨率和检测能力均不好。当30p时，1f和2f处的两个正弦刚刚可以分开，在3f和4f处的两个正弦也可以检出。13因此必须通过提高阶次来达到分辨出间隔较小的频率点的效果。(2)AR模型的自相关法等效于对前向预测的误差序列前后加窗，加窗的结果是使得自相关法的分辨率降低。数据越短，分辨率越不好。3.2协方差法协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取得不同。由于协方差法对于观察区间0,1N外()xn样本并未假定为0，故预测误差功率表达式中的()xnk总是落在观察区间0,1N中，为此预测误差功率的求和上下限必须在,1pN之间。但由此得到的自相关矩阵ˆpR是对此的半正定矩阵，且不具有Toeplitz性质，故不能采用Levinson-Durbin递归算法求解，因此得到的AR模型可能不稳定。该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;tn=-0.5:1/N:0.5-1/N;xpsd=pcov(x,10,N);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plot(tn,fftshift(xpsd));gridon;14图14协方差法10p说明：(1)可以看到谱图在信号源频率处：12340.15,0.16,0.252,0.16ffff谱线狭窄突出，其他处谱线起伏较为剧烈。(2)采用协方差法对信号进行建模，能够较好地反映出信号真正的模型。3.3修正的协方差法AR谱估计的协方差算法基于的是最小化前向预测误差。而修正的协方差算法基于的是最小化前向和后向预测误差。这样使得它的误差功率的计算是在相对于协方差法多一倍的数据点上进行，这在观察数据长度很短的情况下，是非常有利的，但这要求信号在正反两个方向上呈现相同的特性。此外，由此得到的自相关矩阵ˆpR不具有Toeplitz性质，故其正则方程不能采用Levinson-Durbin递归算法求解。该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;fn=-0.5:1/N:0.5-1/N;xpsd=pmcov(x,10,N);15pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plo