《功率谱估计性能分析及Matlab仿真》由会员分享，可在线阅读，更多相关《功率谱估计性能分析及Matlab仿真（20页珍藏版）》请在金锄头文库上搜索。
1、功率谱估计性能分析及Matlab仿真1 引言随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱1。信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方
2、差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算2。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg法等。2 经典功率谱估计经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N个样本数据估计出其功率谱1。2.1 周期图法(Periodogram)Schuster首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。取平稳随机信号x(n)的有限个观察值x(0),x
3、(1),.,x(n1)，求出其傅里叶变换N1XN(ej)x(n)ejnn0然后进行谱估计C1j2S()-XN(ej)N周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用FFT快速算法来计算。但是，这种方法需要对无限长的平稳随机序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄露，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低1o该方法基于Matlab实现的程序：clearall;loadtestx;N=4096;Fn=-0.5:1/N:0.5-1/N;px=f
4、ft(x,N);pmax=max(px);%归一化px=px/pmax;px=10*log10(px+0.000001);plot(Fn,fftshift(px);gridon;图1周期图法N409630111111115-0.3-0.2 -DI D Q 11020 30.40 5图2周期图法N128说明：3的苴(1)本报告仿真中所采用的用于功率谱估计的数据文件来自参考文献test.dat该数据为128点复序列(图3),由复数噪声加上四个复正弦组成。归一化频率分别是：f10.15,f20.16,f30.252,f40.16。图3复序列x(n),n0,127(2)从仿真图可以清晰看到，f1和f2
5、不能完全分开，仅在波形的顶部能看出是两个频率分量；止匕外，当数据长度N太大时(图1),谱曲线呈现较大的起伏；当数据长度N太小时(图2),谱的分辨率又不好。据此，周期图法不满足一致性估计条件。2.2 自相关法(BT法)自相关法的理论基础是维纳辛钦定理。1958年Blackman和Tukey给出了这一方法的具体实现。对于平稳随机信号来说，其自相关函数是确定性函数，故其功率谱也是确定的。这样可由平稳随机离散信号的有限个离散值x(0),x(1),.,x(n1)求出自相关函数Rx(m)x(n)x(nm)然后在(M,M)内对Rx(m)做傅里叶变换，得到功率谱MS()Rx(m)ejnmM该方法基于Matla
6、b实现的程序：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;-OJ-0.3-02-DIQ01020.30.40.5图4自相关法不加窗M64网i4-04-03-0.2-0.1QQ.10.20.3Q.405图5自相关法不加窗M3210-5-10-1520-25刘-3540第融-0 4-0 3 -U 2
7、4J.1 (JQ.20 30 4Q.5图6自相关法使用汉明窗(Hamming)说明：(1)该方法先由序列x(n)估计出自相关函数Rx(m),然后对Rx(m)进行傅里叶变换，便得到x(n)的功率谱估计。当延迟与数据长度之比很小时，可以有良好的估计精度。(2)图4是用自相关法(BT法)求出的功率谱，M64没有加窗；图5也是用自相关法(BT法)求出的功率谱，M32,没有加窗；图6同样是采用自相关法求出的功率谱，M32,使用了汉明窗。显然，自相关函数的延迟M越小，谱变得越平滑。2.3 Welch法段数据用一个合适的窗函数进行平滑处理,最后对各段谱求平均。这样可得功该方法的基本原理是在对随机序列分段时，
8、使每一段有部分重叠，然后对每n 0#率谱其中S()MUL1xm(n)e j(n)这里(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;1060i111111111可一50用amQ2-o.l。5.1U2Q30上图7Welch法不叠合使用汉明窗（Ha
9、mming）#10solIIIIIII-Q.5-0.4Q-3-0.2-5100102Q30.405图8Welch法叠合16点使用汉明窗（Hamming）图9Welch法叠合16点使用矩形窗（Boxcar）#图10Welch法叠合16点使用布莱克曼窗(Blackman)说明：(1)因为Welch法允许各段数据交叠，所以数据段数L会增加，使方差得到更大的改善，但是数据的交叠减小了每一段数据的不相关性，使方差的减小不会达到理论程度。另外，采用合适的窗函数可以减少信号的频谱泄露，同时也可以增加谱峰的宽度，从而提高分辨率。(2)图7是利用Welch法求出的周期图，共分四段，每段32点，没有叠合，使用了汉
10、明窗；图8也是利用Welch法求出的周期图，共分四段，每段32点，使用了汉明窗；图9是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16,且使用了矩形窗；图10是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16,使用了布莱克曼窗。从图中可以看出，由矩形窗处理的谱估计的主瓣宽度最窄，分辨率最好，但是其起伏性较大，所以其方差特性最差。由汉明窗和布莱克曼窗得到的谱估计的主瓣宽度最宽，因此其分辨率相对较差，但其旁瓣较小，大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。因此，选择不同的窗函数其主瓣宽度不一样，造成谱估计的分辨率也不相同。2.4 经典功率谱估计的性能比较
11、由以上的Matlab仿真图形和相关结果分析，我们得到了经典谱估计算法性能的直观比较:(1) 周期图法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰。(2) 自相关法(BT法)由于使用了平滑窗对周期图法估计的功率谱进行了平滑，因此方差性能较好，功率谱比周期图法估计的要平滑，但其分辨率比周期图法低。(3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。综合上述讨论，我们可以对经典谱估计的算法作大致的总结3：(1) 功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，
12、因而仍是目前较常用的谱估计方法。(2) 谱的分辨率较低，它正比于2/N，N是所使用的数据长度。(3) 方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。(4) 周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率且增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。3现代功率谱估计由前一章的讨论我们可知，经典功率谱估计方法的方差性能较差，分辨率较低。而现代谱估计技术的目标都是旨在努力改善谱估计的分辨率。参数模型法是现代谱估
13、计的主要内容，参数模型主要分为AR模型、MA模型和ARMA模型。由于AR模型具有一系列好的性能，因此是被研究最多并获得广泛应用的一种模型。本报告中现代功率谱估计的仿真基于的是AR模型。3.1 自相关法假定观察到得数据为x(0),x(1),.,x(n1)，而对于无法观察到得区间(即n0和nN1)，x(n)的样本假定为0，观测数据区间之外的数据为0，在均方误差意义下使得数据的预测误差功率最小。由于自相关矩阵R?p是Toeplitz矩阵，而且又为正定的，故可利用Levinson-Durbin递归算法高效求解，得到AR模型参数该方法基于Matlab实现的程序:clearall;loadtestx;N=
14、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;5IIIIIIII-0.5-Q.3-0.2-DID01102030.405图11自相关法p10111- -1111iiiiiiiiJ IiaI，i.lj_ _p . :L!yii.一洲TA/ .Zl-一-&叶-i i a 1Jw1二 一 ;1一 :1一JO-15-203Q3540.0.5-0.4Q-3-0.2-0100102Q30.405图12自相关法p20图13自相关法p30说明：(1)图11、12和13是用自相关法求出的AR谱曲线，阶次p分别等于10,20和30。可以看出，在阶次较低时(图11),分辨率和检测能力均不好。当p30时，片和f2处的
