1、功率譜估計的古典算法與現代算法的比較選取周期圖法與Burg算法為例現代信號分析中,對于常見的具有各態歷經的平穩隨機信號,不可能用清楚的數學關系式來描述，但可以利用給定的N個樣本數據估計一個平穩隨機信號的功率譜密度叫做功率譜估計(PSD)。功率譜估計可以分為經典功率譜估計(非參數估計)和現代功率譜估計(參數彳*計)。一、古典功率譜估計古典功率譜估計是將數據工作區外的未知數據假設為零，相當于數據加窗經典功率譜估計方法分為：相關函數法(BT法)、周期圖法以及兩種改進的周期圖估計法。1、相關法相關法是以相關函數為媒介來計算功率譜的，所以又叫間接法，它的理論基礎是維納-辛欽定理。先對數據工作區外的未知數

2、據賦值為零，再由序列x(n)估計出自相關函數R(n),最后對R(n)進行傅立葉變換,便得到x(n)的功率譜估計。2、周期圖法周期圖法是由獲得的N點數據構成的有限長序列直接求fft得其頻譜，取頻譜幅度的平方再除以N,以此作為對x(n)真實功率譜的估計。3、改進的周期圖法改進的周期圖法的主要途徑是平滑和平均。平滑是用一個適當的窗函數與算出的功率譜進行卷積，使譜線平滑，這種方法得出的譜估計是無偏的，方差也小，但分辨率下降；平均就是將截取的數據段再分成L個平均的小段，分別計算功率譜后取功率譜的平均，當L趨于無窮大的時候，L個平均的方差趨于零，可以達到一致譜估計的目的。由于存在旁瓣，會產生兩個后果：一是

3、功率譜主瓣能量泄露到旁瓣使譜估計的方差增大，二是與旁瓣卷積后得到的功率譜完全屬于干擾，嚴重情況下，強信號與旁瓣的卷積可能大于弱信號與主瓣的卷積，使弱信號淹沒在強信號的干擾中無法檢測出來。這是古典法譜估計的主要缺點，即便是改進的周期圖法也無法克服分辨率低的缺點。我們從中選取周期圖法作比較，其算法實現如下：Fs=600;%采樣頻率n=0:1/Fs:1;%產生含有噪聲的序列xn=cos(2*pi*40*n)+cos(2*pi*90*n)+0.1*randn(size(n);n=1:length(xn);figure。)；subplot(2,1,1);plot(n,xn);window=boxcar(

4、length(xn);%矩形窗nfft=1024;Pxx,f=periodogram(xn,window,nfft,Fs);subplot(2,1,2);plot(f,10*log10(Pxx);得到的圖形為：44上i卜'01002003004005006007D0二、現代譜估計參數模型法是現代譜估計中的主要內容，AR模型參數的求解有三種方法：自相關法、Burg遞推算法和改進協方差法。Burg算法不是直接估計AR模型的參數，而是先估計反射系數Km,再利用Levinson關系式求得AR模型的參數。Burg算法采用的數據加窗方法是協方差法，不含有對已知數據段之外的數據做人為的假設。1.其原

5、理如下：Burg算法是使前向預測誤差和后向預測誤差均方誤差之和最小來求取Km的，它不對已知數據段之外的數據做認為假設。計算m階預測誤差的遞推表示公式如下：ffbem(n)=em-i(n)kmem-i(n-1)bbf6m=em-i(n-1)kmem-ifbeo(n)=a(n)=x(n)求取反射系數的公式如下：fb._2Eem-i(n)em-i(n-1)mf2b2Eem-i(n)em-i(n-i)對于平穩隨機過程，可以用時間平均代替集合平均，因此上式可寫成:N-12、em-1(n)em-1(n-1)ikm=-N-1n,m=1,2，pem-1(n/em-1(n-1)inm這樣便可求得AR模型的反射系

6、數。將m階AR模型的反射系數和m-1階AR模型的系數代入到Levinson關系式中，可以求得AR模型其他的p-1個參數。Levinson關系式如下:am(i)=am-1(i)+kmam-1(m-i),i=1,2，m-1m階AR模型白第m+1個參數G,2G:m其中p是預測誤差功率，可由遞推公式mdkm2)求得。易知為進行該式的遞推，必須知道。階AR模型誤差功率P。；?0=Ex2(n) L Rx(0)可知該式由給定序列易于求得。完成上述過程，即最終求得了表征該隨機信號的AR模型的p+1個參數O然后根據Sx(ej0)=<T2H(ej°)即可求得該隨機信號的功率譜密度。2.Burg算法

7、的Matlab運算程序如下：clearallN=512;f1=0.19;f2=0.21;n=1:N;x(n)=sin(2*pi*n*f1)+sin(2*pi*n*f2)+2*randn(size(n);subplot(2,1,1);plot(n,x(n);xlabel('n');ylabel('x(n)');title('兩個正弦信號與白噪聲疊加的時域波形');p=input('InputaNumber>')ef=zeros(p,N);eb=zeros(p,N);de=zeros(p,p);ef(1,:)=x(n);eb(1

8、,:)=x(n);cov(1)=x*x'/N;k(1)=0;form=2:p+1mol=0;den=0;forn=m:Nmol=mol+(-2)*ef(m-1,n)*eb(m-1,n-1);den=den+(ef(m-1,n)A2+(eb(m-1,n-1)A2;endk(m)=mol./den;de(m,m)=k(m);h(1)=cov(1)*(1-de(1,1)A2);forn=m:Nef(m,n)=ef(m-1,n)+k(m).*eb(m-1,n-1);eb(m,n)=eb(m-1,n-1)+k(m).*ef(m-1,n);endendk=k(1,2:p+1);de=de(2:p+

9、1,2:p+1);form=2:pfori=1:m-1de(m,i)=de(m-1,i)+k(m)*de(m-1,m-i);h(m)=h(m-1)*(1-k(m)A2);endendz=de(p,:);s=1,z;n=(0:511)/512;Hw=fft(s,512);subplot(2,1,2);plot(n,h(p)./abs(Hw).A2);xlabel('頻率(Hz)');ylabel('1010g(PSD)');title('基于burg算法的功率譜估計');3.其運行結果如下：P=10兩個正弦信號與白噪聲疊加的時域波形10020030

10、0400500600EX基于Mg算法也功率譜估計3020100nsd)J°40P=30陽帆桐神州1D0200300400500600基于hurg算法幅功率譜估計0.10.20.3D.J0506070.80.9頻率(Hz)ofl-oo321QSCL)號0LQsd60_0_6000100200300400500n基于bug算法的功率譜估計 10050Qsdwa已 L9.,8.n-,6.5LI4 0.30.,20.1 0.P=50105©0x-5-1DP=60基于bug算法的功率譜估計ogoo一P=80基于burg算法的功率譜估計150rrIiiir616口豆5000,1020,

11、3040.50.6070.80.91頻率(Hz)P=100兩個正弦信號與白噪聲疊加的時域波形琳M曬怖視仙100200300400500600基于burg算法懦功率譜估計1U6口亙50P=120基于bu項算法的功率譜估計01020.3040.50.B070.80.91頻率(Hz).JLooQoooooo4321nsa.)號0LP=130基于hurg算法的功率譜估計(asd)Bo豆P=160兩個正弦信號與白噪聲疊加的時域波形100200300400500600n基于bu嗚算法的功率譜估計5 0 石EX10O100.10.20.30.40.50.607 D.B 0.911頻率（Hz）6000QsdwQoL結果分析：1、從圖中我們可以清晰的看到Burg算法的優越性，Burg算法求解AR模型的過程是非常穩定的，而且具有很高的分辨率。當然對于Burg算法來說，P即階數的選擇是至關重要的，我們從實驗結果圖中可以看出，當P介于50和80之間時，得到的頻譜圖是較優越的，P在130左右時頻譜圖是最優越的，這也符合了經驗定理，