1、參數(shù)法功率譜估計(jì)一、 信號的產(chǎn)生（一）信號組成在本實(shí)驗(yàn)中，需要事先產(chǎn)生待估計(jì)的信號，為了使實(shí)驗(yàn)結(jié)果較為明顯，我產(chǎn)生了由兩個不同頻率的正弦信號（頻率差相對較大）和加性高斯白噪聲組成的信號。（二） 程序N=1024;n=0:N-1;xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+randn(1,1024);這樣就產(chǎn)生了加有白噪聲的兩個正弦信號其波形如下二、參數(shù)模型法功率譜估計(jì)（一）算法原理簡介1參數(shù)模型法是現(xiàn)代譜估計(jì)的主要內(nèi)容，思路如下： 假定所研究的過程是由一個白噪聲序列激勵一個因果穩(wěn)定的可逆線性系統(tǒng)的輸出； 由已知的，或其自相關(guān)函數(shù)估計(jì)的參數(shù)； 由的參數(shù)來估計(jì)的

2、功率譜。2自回歸模型，簡稱AR模型，它是一個全極點(diǎn)的模型?！白曰貧w”的含義是：該模型現(xiàn)在的輸出是現(xiàn)在的輸入和過去p個輸出的加權(quán)和。此模型可以表現(xiàn)為以下三式： ； ； 。3AR模型的正則方程建立了參數(shù)和的自相關(guān)函數(shù)的關(guān)系，公式如下： 時， 時。（二）兩種AR模型階次的算法1Yule-Walker算法(自相關(guān)法)（1）算法主要思想Yule-Walker算法通過解Yule-Walker方程獲得AR模型參數(shù)。從低階開始遞推，直到階次p，給出了在每一個階次時的所有參數(shù)。公式如下： ； ； 。（2）運(yùn)算簡要框圖估計(jì)功率譜密度（2M-1）點(diǎn)FFTX(n)求自相關(guān)估計(jì)參數(shù)輸出Yule-Walker法譜估計(jì)運(yùn)算

3、簡要框圖(3) 程序示例clear all;close all;N=512;n=0:N-1;xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+2*randn(1,512);Rx=zeros(1,N+1);%從課本上的公式來看，Rx（m）中的m屬于（0，m），即共有m+1個，故在這里設(shè)Rx是一個一行，N+1列的向量figure(1)plot(n,xn);title(（a）兩正弦信號加白噪聲波形)%下面用書中所講自相關(guān)函數(shù)估計(jì)中的漸進(jìn)無偏估計(jì)來估計(jì)自相關(guān)函數(shù)for m=1:N+1;%由于在matlab中，下角標(biāo)不能是0，m屬于（0，m），在此只能從1到N+1 sum=

4、0; for n=1:(N+1-m);%同樣道理，把書中公式里m換成m-1,N換成N+1，求和下限變?yōu)? sum=sum+xn(n).*xn(n+m-1); end Rx(m)= sum/N;%切記，這里的Rx(1)才是自相關(guān)函數(shù)在0點(diǎn)的取值。Rx（m）只是一個存儲數(shù)據(jù)的代號，為了跟書中公式一致，才叫Rxend%下面估計(jì)各參數(shù)P=50;a=zeros(P,P);%a中有兩個變量m,i，所以設(shè)a是P行P列的向量km=zeros(1,P);%因?yàn)殡A次是P，故反射系數(shù)有P個p=zeros(1,P+1);%由于matlab中沒有，故用p來代替表示，的范圍是（0，P）共有P+1個%下面計(jì)算AR模型參數(shù)的

5、各個初始化值p(1)=Rx(1);a(1,1)=-Rx(2)/Rx(1);km(1)=a(1,1);p(2)=Rx(1).*(1-abs(a(1,1).2);for m=2:P %由于m=1時的各個值在上面已經(jīng)給出，故從m=2開始求 sum1=0; for i=1:m-1 sum1=sum1+a(m-1,i).*Rx(m-i+1); end a(m,m)=-(Rx(m+1)+sum1)/p(m);%km(m)=a(m,m);求出km for i=1:m-1 a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i); end p(m+1)=p(m)*(1-abs(a(m,m).2);e

6、ndz=1,a(P,:);G=sqrt(p(P);H w=freqz(G,z,512);%調(diào)用計(jì)算數(shù)字濾波器頻響的函數(shù)figure(2)plot(w/(2*pi),10*log10(abs(H).2);title(自相關(guān)法)ylabel(10log(PSD)title(b) yule-walker法估計(jì)功率譜密度 )（4）結(jié)果分析從波形圖中可以十分清晰的分辨出兩個不同頻率的正弦波2Burg法（1）算法主要思想Burg法不是直接估計(jì)AR模型的參數(shù)，而是先估計(jì)反射系數(shù)。使用線性預(yù)測的方法來計(jì)算不同階數(shù)下的反射系數(shù)，其同時使用前向和后向線性預(yù)測，使前向和后向預(yù)測誤差的平均功率相對各階反射系數(shù)最小，將

7、反射系數(shù)代入Levinson-Durbin公式即可求解。（2）運(yùn)算簡要框圖 用L-D公式求解估計(jì)反射系數(shù)X(n) 輸出（3）程序示例clear all;close all;N=512;n=0:N-1;xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+randn(1,512);Rx=zeros(1,N+1);%從課本上的公式來看，Rx（m）中的m屬于（0，m），即共有m+1個，故在這里設(shè)Rx是一個一行，N+1列的向量figure(1)plot(n,xn);title(（a）兩正弦信號加白噪聲波形)%下面用書中所講自相關(guān)函數(shù)估計(jì)中的漸進(jìn)無偏估計(jì)來估計(jì)自相關(guān)函數(shù)for

8、m=1:N+1;%由于在matlab中，下角標(biāo)不能是0，m屬于（0，m），在此只能從1到N+1 sum=0; for n=1:(N+1-m);%同樣道理，把書中公式里m換成m-1,N換成N+1，求和下限變?yōu)? sum=sum+xn(n).*xn(n+m-1); end Rx(m)= sum/N;%切記，這里的Rx(1)才是自相關(guān)函數(shù)在0點(diǎn)的取值。Rx（m）只是一個存儲數(shù)據(jù)的代號，為了跟書中公式一致，才叫RxendP=50;a=zeros(P,P);%a中有兩個變量m,i，所以設(shè)a是P行P列的向量p=zeros(1,P+1);%由于matlab中沒有，故用p來代替表示，的范圍是（0，P）共有P+

9、1個%下面計(jì)算AR模型參數(shù)a(1,1)=-Rx(2)/Rx(1);ef=zeros(P,N);eb=zeros(P,N);ef(1,:)=xn;eb(1,:)=xn;for m=2:P+1; km1=0;km2=0; for n=m:N km1=km1+ef(m-1,n).*eb(m-1,n-1); km2=km2+(ef(m-1,n).2+(eb(m-1,n-1).2; end a(m,m)=(-2)*km1./km2; for n=m:N ef(m,n)=ef(m-1,n)+a(m,m).*eb(m-1,n-1); eb(m,n)=eb(m-1,n-1)+a(m,m).*ef(m-1,n)

10、;end endp(1)=Rx(1);p(2)=Rx(1).*(1-abs(a(1,1).2);a=a(2:P+1,2:P+1);for m=2:P %由于m=1時的各個值在上面已經(jīng)給出，故從m=2開始求 for i=1:m-1 a(m,i)=a(m-1,i)+a(m,m)*a(m-1,m-i); end p(m)=p(m-1)*(1-abs(a(m,m).2);endz=1,a(P,:);G=sqrt(p(P);H, w=freqz(G,z,512);%調(diào)用計(jì)算數(shù)字濾波器頻響的函數(shù)figure(2)plot(w/(2*pi),10*log10(abs(H).2);ylabel(10log(P

11、SD)title(b)Burg法估計(jì)功率譜密度 ) （4）結(jié)果分析下面是程序運(yùn)行后的結(jié)果從上圖中我們同樣可以看到信號中有兩個頻率分量，一個在0.2處，一個在0.213左右處，與產(chǎn)生的信號相一致。（三）、兩種算法的比較由上面兩圖我們可以看出Burg算法得出的功率譜更平坦些。（四）與古典法功率譜估計(jì)的比較512點(diǎn)在上面這兩張圖中我們就不能夠很好的分辨兩個靠的比較近的譜峰，在同樣條件下參數(shù)法功率譜估計(jì)卻能夠?qū)崿F(xiàn)譜峰的分辨，充分體現(xiàn)了參數(shù)法功率譜估計(jì)的優(yōu)越性。下面減小一下點(diǎn)數(shù)，再進(jìn)行分析256點(diǎn)古典法 256點(diǎn)參數(shù)法功率譜估計(jì)由上面兩圖的對比我們可以看到，當(dāng)采樣的點(diǎn)數(shù)減少的時候，古典法的估計(jì)就十分不準(zhǔn)確了，而參數(shù)法功率譜估計(jì)卻還可以讓人清晰地分辨兩個挨得較近的譜峰。優(yōu)點(diǎn)十分突出。一些體會在持續(xù)了半月的做作業(yè)期間內(nèi)，我自己可以很清楚的看到自己的進(jìn)步，從什么都不太懂到現(xiàn)在通過自己的努力可以編出相應(yīng)的程序