功率谱分析

随机信号处理之功率谱估计姓名:***学号：************专业：电子信息工程一．原理古典谱估计之相关函数法相关法谱估计是以相关函数为媒介来计算功率谱,又叫做间接法它的理论基础是维纳-辛钦定理,其具体步骤如下:第一步,由获得的N点数据构成的有限长序列xn(n)来估计自相关函数,即:第二步,由自相关函数的傅里叶变换求功率谱,即以上两步经历了两次截断,一次是估计时仅利用了x(n)的N个观测值,这相当于对x(n)加矩形窗截断.该窗是加在数据上的,一般称为加数据窗.另一次是估计时仅利用了从-(M-1)到(M-1)的,这相当于对加矩形窗截断,将截成(2M-1)长,这称为加延时窗.式中和分别表示对和的估值.由于M≪N,使得上式的运算量不是很大,因此在FFT问世之前,相关法是最常用的谱估计方法古典谱估计之周期图法首先，数据加窗得到有限长序列𝑥𝑁(n)直接求傅里叶变换，得频谱𝑋𝑁(eiω),即𝑋𝑁(𝑒𝑖𝜔)=∑𝑥𝑁(𝑛)𝑒−𝑖𝜔𝑛𝑁−1𝑛=0然后，取频谱幅度的平方，并除以N，以此作为对x(n)真实功率谱Sx(eiω)的估计，即Sx(eiω)=1𝑁|𝑋𝑁(𝑒𝑖𝜔)|2事实上，周期图法谱估计与相关法谱估计的差异只是估计自相关函数的方法不同。参数模型法谱估计之Yule-Walker方程矩阵估计高斯消去法是解尤勒-沃克方程的一种算法,这种方法直接求解线性方程组,运算量较大。在解之前,先大概了解一下尤勒-沃克方程:N1xNNn01R(m)x(n)x(nm)N1(1)ˆˆS()()MjjmxxmMeRmexR(m)()NxnˆS()jxexR(m)xR()mxR()mxR(m)ˆS()jxexR()mS()jxe如前所述,P阶ＡＲ模型的系统函数为可以看出，Ｐ阶ＡＲ模型有Ｐ＋１个待定系数，由上式，可得白噪声激励得到的系统输出可以理解为，用ｎ时刻之前的ｐ个值的线性组合来预测ｎ时刻的值预测误差为．在均方误差最小准则下，组合系数的选择应使预测误差的均方值最小．经过一系列的运算，最终可以得到ＡＲ模型的正则方程也就是尤勒－沃克（Yule-Walker）方程．由P个方程可以求出P个参数𝑎𝑖.有了参数𝑎𝑖（i=1,2,3,4,…,p），就可以根据自相关函数和参数𝑎𝑖求系统增益G。再由公式𝐒𝐲(𝐞𝐢𝛚)=𝐒𝐱(𝐞𝐢𝛚)∗|𝐇(𝐞𝐢𝛚)|𝟐可以得出功率谱。参数模型法谱估计之levinson-durbin快速递推法Levinson-durbin递推算法是解尤勒-沃克方程的快速有效的算法,这种方法利用自相关矩阵具有的一些好的性质,是运算量大大减少,在介绍Levinson-durbin递推算法之前,先大概了解一下尤勒-沃克方程:如前所述,P阶ＡＲ模型的系统函数为可以看出，Ｐ阶ＡＲ模型有Ｐ＋１个待定系数，由上式，可得白噪声激励得到的系统输出1G()1miiiHzaz1()()()piixnaxniGn1()piiaxni()xn()Gn12,,,paaa()Gn121(),1,2,()(),0pixixpixiaRmimpRmaRiGm1G()1miiiHzaz1()()()piixnaxniGn可以理解为，用ｎ时刻之前的ｐ个值的线性组合来预测ｎ时刻的值预测误差为．在均方误差最小准则下，组合系数的选择应使预测误差的均方值最小．经过一系列的运算，最终可以得到ＡＲ模型的正则方程也就是尤勒－沃克（ｙｕｌｅ－ｗａｌｋｅｒ）方程．在上面方程中令ｐ＝１，得一阶ＡＲ模型的尤勒-沃克方程为可以解得进而求得在一阶的基础上进行递推，将阶次为ｍ时的第ｍ个系数定义为反射系数，用表示，于是可以将计算ｍ阶ＡＲ模型参数的Levinson-durbin递推算法表示如下：式中，．1()piiaxni()xn()Gn12,,,paaa()Gn121(),1,2,()(),0pixixpixiaRmimpRmaRiGm111(0)(1)(1)(1)(1)(0)0xxxxRaRRaR1(1)(1)(0)xxRaR2211(1)(0)(0)[1(1)](0)xxxxRRRaR()mammk11111121[()()()]()()(),1,2,3,,1(1)mxmximmmmmmmmmRmaiRmikaiaikamiimk1,2,mp20[()](0)xExnR参数模型法谱估计之Burg算法Burg算法不是直接估计AR模型的参数，而是先估计反射系数Km，再利用Levinson关系式求的AR模型的参数。估计反射系数所依据的准则是使前向预测误差FPE功率𝜀𝑓和后向预测误差BPE功率𝜀𝑏的平均值ε最小，即ε=12(𝜀𝑓+𝜀𝑏)=12{∑{[𝑒𝑚−1𝑓(𝑛)]2𝑁−1𝑛=𝑚+∑[𝑒𝑚−1𝑏(𝑛−1)]2𝑁−1𝑛=𝑚}式中各阶前向预测误差和后向预测误差由递推公式求出，即{𝑒𝑚𝑓(𝑛)=𝑒𝑚−1𝑓(𝑛)+𝑘𝑚𝑒𝑚−1𝑏(𝑛−1)𝑒𝑚𝑏(𝑛)=𝑒𝑚−1𝑏(𝑛−1)+𝑘𝑚𝑒𝑚−1𝑓(𝑛)𝑒0𝑓(𝑛)=𝑒0𝑏(𝑛)=𝑥(𝑛)得到反射系数的估计公式为：𝑘𝑚=2∑[𝑒𝑚−1𝑓(𝑛)𝑒𝑚−1𝑓𝑏(𝑛−1)]𝑁−1𝑛=𝑚∑{[𝑒𝑚−1𝑓(𝑛)]2+[𝑒𝑚−1𝑏(𝑛−1)]2}𝑁−1𝑛=𝑚上式中预测误差的求和范围表明，Burg算法所采用的数据加窗方法是协方差法，不含有对已知数据段之外数据的人为假设。将m阶AR模型的反射系数和m-1阶AR模型的系数代入到Levinson关系式中，可以求得AR模型其他的p-1个参数。Levinson关系式如下：1-m1,2,ii),-(m(i)(i)1-m1-mm，akaamm阶AR模型的第m+1个参数G，m2G其中m是预测误差功率，可由递推公式)-(12m1-mmk求得。易知为进行该式的递推，必须知道0阶AR模型误差功率00=(0)(n)Ex2Rx可知该式由给定序列易于求得。完成上述过程，即最终求得了表征该随机信号的AR模型的p+1个参数。然后根据)(ejxS=2j2)H(e即可求得该随机信号的功率谱密度。二．程序，图像及分析时域波形图i.古典谱估计之自相关法程序：Fs=1000;N=1024;n=0:N-1;xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);cxn=xcorr(xn);CXk=fft(cxn,N);Pxx=abs(CXk);index=0:round(N/2-1);k=index*Fs/N;figure(1);plot(k,Pxx(index+1));图像：图一古典谱估计之周期图法程序：Fs=1000;N=1024;n=0:N-1;xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);x1n=fft(xn);x2n=abs(x1n);x3n=x2n.*x2n;Pxx2=abs(x3n);index=0:round(N/2-1);k=index*Fs/N;figure(2);plot(k,Pxx2(index+1));图像：图二由图一和图二可以看出，我们利用这两种方法可以分辨出两根谱线，分别位于f=200HZ和f=213HZ附近，根据Fs=1000HZ可知，谱线位于w/pi=0.4及w/pi=0.426处，而且周期图法优于自相关法。这两种方法都是用获得的N个数据对随机过程进行功率谱估计，它隐含着对无限长数据加了一个矩形窗截断。时域与矩形窗相乘对应于频域中与矩形窗频谱卷积，故估计谱就相当于真实谱与矩形窗频谱卷积的结果。对功率谱有影响的是矩形窗频谱的幅度谱，它存在两个问题：一是主瓣不是无限窄，二是有旁瓣。由于主瓣不是无限窄，真实的功率谱与主瓣卷积后将使功率向附近频域扩散，是信号模糊，降低了谱分辨率，主瓣越宽分辨率越低。N减小时，主瓣宽度变宽，以上两图是取N=1024时得出的，此时我们假设N=128可以看出无法分辨出两根谱线。（如下图三四）由于存在旁瓣，又会产生两个后果，一是功率谱主瓣内的能量“泄露”到旁瓣使谱估计的方差增大，二是与旁瓣卷积后得到的功率谱完全属于干扰。严重情况下，强信号与旁瓣的卷积有可能大于弱信号与主瓣的卷积，使弱信号淹没在强信号的干扰中，无法检测出来。我们将w/pi=0.426处谱线的信号强度调大为w/pi=0.4处的20倍可以看出此时分辨不出两根谱线（如图五六）。可见两种方法都不是对Sx(eiω)的一致估计。于是我们利用平滑和平均的方法改进周期图法。图三图四图五图六平滑是用一个适当的窗函数W(eiω)与算出的功率谱进行卷积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小但分辨率下降。平均是将截取的数据段𝑥𝑁(n)再分成L个小段，分别计算功率谱后取平均。因为L个平均的方差比随机变量单独的方差小L倍，所以当L→∞时，L个平均方差趋于0，可以达到一致估计的目的。平均的程序：N=1024;n=0:N-1;xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);Nfft=1024;K=4;L=N/K;Sx=zeros(1,Nfft/2);fork=1:Kks=(k-1)*L+1;ke=ks+L-1;X=fft(x(ks:ke),Nfft);X=abs(X).*abs(X);fori=1:Nfft/2Sx(i)=X(i)+Sx(i);endendfori=1:Nfft/2Sx(i)=Sx(i)/N;endindex=0:round(Nfft/2-1);f=index/Nfft;plot(f,Sx(index+1));图像：图七图八显然在平均的时候（如图七），分辨率降低，且分段越多分辨率越低，如图八为分段数是14时，此时无法分辨出两根谱线。所以才采取重叠分段法，同时改变窗函数的种类观察，如下面程序即图九N=1024;n=0:N-1;xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);Nfft=1024;K=4;D=fix(N/(K+1));L=2*D;Sx=zeros(1,Nfft/2);wn=hamming(L);w=wn';fork=1:Kks=(k-1)*D+1;ke=ks+L-1;xk=x(ks:ke).*w;X=fft(xk,Nfft);X=abs(X).*abs(X);fori=1:Nfft/2Sx(i)=X(i)+Sx(i);endendfori=1:Nfft/2Sx(i)=Sx(i)/N;endindex=0:round(Nfft/2-1);f=index/Nfft;plot(f,Sx(index+1));图九在k=14时这种方法可以分辨出两根谱线，分辨率提高，如下图十图十参数模型法谱估计之Yule-Walker方程矩阵估计程序：N=1024;n=0:N-1;Fs=500;x=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);form=1:Ns=0;forn=1:N-(m-1)s=s+x(n)*x(n+m-1);endR(m)=s/N;endp=100;A=ones(1,p);fori=1:p-1r(1,i)=A(1,i)*R(i);endT=toeplitz(r);B=ones(p,1);fori=1:pY(i,1)=B(i,1)*R(i+1);enda=-T\Y;G=r(1,1);fori=1:pG=G+a(i,1)*R(i+1);endX=1;K=[1a'];[H,w]=freqz(X,K);Sx=G*(abs(H).*abs(H));figu