功率谱估计浅谈讲解

1、功率谱估计浅谈摘要：介绍了几种常用的经典功率谱估计与现代功率谱估计的方法原理，并利用Matlab对随机信号进行功率谱估计，对两种方法做出比较，分别给出其优缺点。关键词：功率谱；功率谱估计；经典功率谱估计；现代功率谱估计前言功率谱估计是从频率分析随机信号的一种方法,一般分成两大类:一类是经典谱估计;另一类是现代谱估计。由于经典谱估计中将数据工作区以外的未知数据假设为零，这相当于数据加窗，导致分辨率降低和谱估计不稳定。现代谱估计则不再简单地将观察区外的未知数据假设为零，而是先将信号的观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。周期图、自相关法及

2、其改进方法(Welch)为经典(非参数)谱估计方法，其以相关和傅里叶变换为基础,对于长数据记录较适用,但无法根本解决频率分辨率低和谱估计稳定性的问题，特别是在数据记录很短的情况下，这一问题尤其突出。以随机过程的参数模型为基础的现代参数法功率谱估计具有更高的频率分辨率和更好的适应性，可实现信号检测或信噪分离，对语音、声纳雷达、电磁波及地震波等信号处理具有重要意义，并广泛应用于通信、自动控制、地球物理等领域。在现代参数法功率谱估计方法中，比较有效且实用的是AR模型法，Burg谱估计法,现代谱估计避免了计算相关，对短数据具有更强的适应性，从而弥补了经典谱估计法的不足，但其也有一些自身的缺陷。下面就给

3、出这两类谱估计的简单原理介绍与方法实现。经典谱估计法经典法是基于传统的傅里叶变换。本文主要介绍一种方法：周期图法。周期图法由于对信号做功率谱估计，需要用计算机实现，如果是连续信号，则需要变换为离散信号。下面讨论离散随机信号序列的功率谱问题。连续时间随机信号的功率谱密度与自相关函数是一对傅里叶变换对，即：S(0)=JxR(t)ejdxg若R(m)是R(0)的抽样序列，由序列的傅里叶变化的关系，可得xxS(ej)二工R(m)e-jnxxmg即S(e；o)与R(m)也是一对傅里叶变换对。显然，由序列傅里叶的频谱特性可知xxS(e；o)是以2“为周期的。而实际计算只能从离散随机信号序列x(n)的有限长

4、x(长度为N)的数据来对S(ej)与R(m)进行估计。设有限长离散序列为x(n)，则:xx1R(m)二x(m)*x(m)xNNNNS(ej)二DFTR(m)xNxN由DFT的下列卷积特性:若X(ej)DFTX(m)，则:NX*(ej)DFTX(m)N从而:DTFTR(m)丄DFTx(m)DFTx(m)xNNNN即S(ej)=X(ej)X*(ej)丄X(ej)2xnNN综上所述，先用FFT求出随机离散信号N点的DFT，再计算幅频特性的平方，然后除以N,即得出该随机信号的功率谱估计。由于这种估计方法在把R(t)离散x化的同时，使其功率谱周期化，故称之为“周期图法”，也称为经典谱估计法。周期图法进行

5、谱估计，是有偏估计，由于卷积的计算过程会导致功率谱真实值的尖峰附近产生泄漏，相对地平滑了尖峰值，因此造成谱估计的失真。另外，当NT*时，功率谱估计的方差不为零，所以不是一致性估计。并且功率谱估计在等于2“/N整数倍的各数字频率点互不相关。其谱估计的波动比较显著，特别是当N越大、2兀/N越小时，波动越明显。但如果N取得太小，又会造成分辨率的下降。图1.原始信号1频率用402020-HS5SS图2原始信号1的功率谱估计originalsignal图3.原始信号2周期團N=10241002003004005006007008009001000频率屁ODD2-2mp、B&曰图4.原始信号2的功率谱估计

6、oO32OOLOO图5.平均周期图法(4*256)平均周期團(重1J2)N=102450100150200250300350400450500ooooO4321mp、無掛肓图6.平均周期图法(重叠一半)图1所示的信号为xn=sin(2*兀*50*t)+2sin(2*兀*120*t)+randn(1,N)，其中,randn是正态分布随机数组，N为256,t是从0到1,dt为1/256。图2为该信号的功率谱估计。图2所示的信号为xn=sin(2*兀*50*t)+2sin(2*兀*120*t)+randn(1,N)，其中，randn是正态分布随机数组，N为1024,t是从0到1,dt为1/1024。

7、图4为该信号的功率谱估计。图5是将图2所示的信号分为四段，每段的范围分别为(1,256),(257,512),(513,768),(768,1024).每一道都没有重叠。然后对分段分别作傅里叶变换，再把功率谱加起来做平均，得到图5。图6是将图2所示的信号分为六段，分别为(1:256)，(129:384)，(257:512)，(385:640)，(513:768)，(641:896)，(769:1024)。每两段之间都重合一半。图1和图3相比，图1较为平滑，相应的，图1的功率谱也比较平滑。图5和图6比,图6较为平滑，这是因为图6的谱是六段的平均。图8.Bartlett平均周期图法现代谱估计法现代

8、参数法功率谱估计方法中，比较有效且实用的是AR模型法，Burg谱估计法,在本文中介绍的是AR模型法。AR模型法经典谱的主要缺点是频率分辨率低。这是由于周期图法在计算中把观测到的有限长的N个数据以外的数据认为是零，这显然与事实不符。如果把已观测到的数据估计出一白噪声激励，就不必认为N个以外的数据全为零，就有可能克服经典谱估计的缺点。一个实际中的随机过程总是可以用以下模型很好的表示：Ybz-iii=11+Yaz-knk=1当除b外的所有b均为零时的形式称为p阶自回归模型即AR模型，又称为全极点0i模型。P（）=xx当方差为b2的白噪声通过AR模型时，输出的功率谱密度为:m:1+艺ae-jmkk=1

9、若已知参数a,a,a及b2，就可以得到信号的功率谱估计。它们之间是12pm图10.自相关函数的无偏估计1000so1012026303640图11.AR模型求出的功率谱图7所示的信号长度为155个点，长度较周期图法里的信号短，信号为X=sin(2*兀*20*n)+sin(2*兀*30*n)+(1/20)*randn(size(n)，其中，n为0155/100,间隔为1/100。图8为自相关函数的无偏估计。图9为功率谱。从图9可以看出，这种方法具有极高的分辨能力。只是在20和30处有两个峰值，在其它地方的值为零。将信号改变成以下信号：X=sin(2*兀*F1*n)+sin(2*兀*F2*n)+(

10、1/20)*randn(size(n)与经典谱估计方法进行一个对比（基于第一个信号）：相比与前面几种方法得到的结果来看，相差非常大，不仅仅是分辨率的提高，其余无效噪音或者说扰动都是非常小的。结论通过实验可以直接看出以下特性：1）经典功率谱估计的方差大，谱分辨率差。采用经典的傅里叶变换及窗口截断，对长序列有良好估计。2）现代谱的分辨率较高。这是由于在时域的开窗，使得在频域发生“频谱泄漏”，即功率谱的主瓣能量泄漏到旁瓣中，导致弱信号的主瓣被强信号的旁瓣所湮没，造成谱的模糊。3）平均周期图法的收敛性较好，曲线平滑，但是功率谱主瓣较宽，分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gib

11、bs现象而造成的。参考文献刘志刚，李录明，赵冬梅.现代谱估计法及应用效果J.石油地球物理勘探，2009,S1:5-9+167+9.樊剑，刘铁，胡亮.基于现代时频分析技术的地震波时变谱估计J.振动与冲击，2007,11:79-82+98+184.蔡希玲，刘学伟，吕英梅，曹锡娜.统计F-X谱估计方法及应用J.勘探地球物理进展，2008,03:181-186+163.李刚.宽带信号空间谱估计算法研究D.哈尔滨工程大学,2007.姚武川，姚天任.经典谱估计方法的MATLAB分析J.华中理工大学学报,2000,04:45-47.王凤瑛，张丽丽.功率谱估计及其MATLAB仿真J.微计算机信息,2006,3

12、1:287-289.王福杰，潘宏侠.MATLAB中几种功率谱估计函数的比较分析与选择J.电子产品可靠性与环境试验，2009，06:28-31.附录：谱估计的Matlab实现程序1：经典法谱估计clf;Fs=1000;N=256;Nfft=256;%数据的长度和FFT所用的数据长度n=O:N-l;t二n/Fs;%采用的时间序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1，N);Pxx=10*loglO(abs(fft(xn,Nfft).“2)/N);%Fourier振幅谱平方的平均值，并转化为dBf=(0:length(Pxx)T)*Fs/length(

13、Pxx);%给出频率序列subplot(2,l,l),plot(f,Pxx);%绘制功率谱曲线xlabel(频率/Hz);ylabel(功率谱/dB);title(周期图N=256);gridon;程序2：经典谱加窗分析Fs=1000;%采样频率n=0:1/Fs:1;%产生含有噪音的信号xn=sin(2*pi*50*n)+2*sin(2*pi*120*n)+randn(size(n);plot(n,xn);window二boxcar(length(xn);%矩形窗nfft=512;Pxx,f=periodogram(xn,window,nfft,Fs);%直接法figure(1)plot(f,

14、10*log10(Pxx);window二boxcar(length(xn);%矩形窗nfft=1024;Pxx,f=periodogram(xn,window,nfft,Fs);%直接法figure(2)plot(f,10*log10(Pxx);程序3：Bartlett平均周期图法fs=1000;n=0:1/fs:1;xn=sin(2*pi*50*n)+2*sin(2*pi*120*n)+randn(size(n);nfft=1024;window=hamming(nfft);noverlap=0;p=0.9;Pxx,Pxxc=psd(xn,nfft,fs,window,noverlap,p

15、);index=0:round(nfft/2-1);k=index*fs/nfft;plot_Pxx=10*log10(Pxx(index+1);plot_Pxxc=10*log10(Pxx(index+1);figure(1)plot(k,plot_Pxx)figure(2)plot(k,plot_Pxxplot_Pxx-plot_Pxxcplot_Pxx+plot_Pxxc);程序4：现代谱估计法clearall%产生信号FS=100;%设采样频率为100,则根据Fl/FS=0.2,F2/FS=0.3OR0.25，可以得到F1,F2便于计算；N=155;%数据长度改变数据长度会导致分辨率的

16、变化；n=0:1/FS:N/FS;F1=20;%第一个SIN信号的频率；F2=30;%第二个SIN信号的频率，取25或者30；X=sin(2*pi*F1*n)+sin(2*pi*F2*n)+(1/20)*randn(size(n);%产生信号,由信噪比为10dB推出噪声功率；信号长度从X(1)至【X(N+1)；%产生自相关函数数组form=1:N+1%初始化R(m),R(m)用来存放自相关函数；由于R(0)在MATLAB里无效，所以从1开始到N+1；R(m)=0;endform=1:N+1%做自相关函数的无偏估计；S=0;forn=1:N+2-mH=X(n)*X(n+m-1);S=S+H;en

17、dR(m)=S/N;end%估计完毕subplot(3,1,1);plot(X);subplot(3,1,2);plot(R);%Levinson算法主体部分fork=1:N+1%初始化后面算法中要用到的数组,其中a(k,k)用来存放AR模型参数，FPE(k)是最终预测误差准则函数，a(k,k)=0;%U2(k)是AR模型中的另一参数即误差功率；FPE(k)=0;U2(k+1)=0;endU2(1)=R(1);a(l,l)=-R(2)/R(l);%由Levinson算法，计算一阶模型参数a11；U2(2)=(l-(abs(a(l,l)厂2)*U2(1);%由Levinson算法，计算一阶模型参

18、数误差功率；FPE(1)=U2(2)*(N+2)/N;%计算一阶模型的最终预测误差准则函数；S=0;fork=2:N%由Levinson梯归公式计算K阶模型参数和FPE函数；forl=1:k-1M=a(k-1,l)*R(k-l+1);S=S+M;enda(k,k)=-(R(k+1)+S)/U2(k);fori=1:k-1a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);endU2(k+l)=(l(abs(a(k,k)厂2)*U2(k);FPE(k)=U2(k+1)*(N+k+1)/(N-k+1);S=0;end%梯归计算完毕；%确定阶数Pmin=FPE(1);%求出使得FPE函数取最小值的阶数P；fork=2:NifFPE(k)minmin=FPE(k);p=k;endend%p=2%为了调整效果可以在这里自行指定阶数；