多種功率譜估計鄭州大學隨機信號處理大作業(yè)

隨機信號處理大作業(yè)多種功率譜估計的算法實現(xiàn)及性能比較一、引言二、原理及過程Xx(ejw)即ep(n)=ep-1(n-1)+k,eg-1(n)=1,2,…,p(3.4)①由初始條件ed(n)=eg(n)=x(n)和式(3.4)求出k?③由k?和式(3.3)求出ef(n),e{(n),再由式(3.4)估計k?即e=[1eiw.…ei(N-1)w]T若令a;=1(j=M+1,M+2,…,N),則計算功率譜人為的設置信噪比范圍為-15db~15db,通過分別調用上面的求功率譜的函數(shù)，求特定信噪比下1000次功率譜估計的均方誤差，分別繪制采用上述功率譜估計時，均方誤差①求古典法在-15db下的1000次功率譜估計的均方誤差②同上求-14db~15db下的1000次功率譜估計的均方誤差③繪制均方誤差與信噪比的關系④同理，求出采用上文中現(xiàn)代譜估計時，均方誤差與信噪比的關系結果分析原始信號為信噪比和均方誤差的關系，并分析各種功率譜估計方法的性能優(yōu)劣。1,古典法N不小于130時可以識別出兩個波峰n*nnn=128f1=0.20002,直接解yole-walker方程n=130f1=0.2000f2=2.2126P=10時，無論N取多少，均不能識別兩個波峰P=40N=128f1=0.1998fP=20N=256f1=0.1995p=40n=256f1=0.1995f2=0.2135P=20N=128f1=0.1992P=40N=128f1=0.1998f2=0.2140P=404Burg算法N=256f1=0.1995f2=0.2135P=10時，N取任何值都不能識別出兩個波峰P=10N=128f1=0.2031P=10N=256f1=0.2027P=20時，N不小于55時可識別出兩個波峰P=20N=128f1=0.2003f2=0.2138P=20N=256f1=0.2005f2=0.2125P=40時，N大于50時可識別出兩個波峰P=40N=128f1=0.2005f2=0.2138P=40N=256f1=0.2008f2=0.2130在相同的條件下，f1=0.2,f2=0.213,N=128時，觀察圖形發(fā)現(xiàn)，用古典法只能估計一個頻率點，用AR模型估計時，當階次P達到一定時，仍能分辨出兩頻率點，所以只要選擇合適的階次，就能用AR模型估計兩接近的頻率點，而N越大，圖像越精細，越易于分辨兩頻率點，對于古典法，當N合適時，仍能分辨相近的頻率點。由圖易知，分辨率的關系為：古典法<直接解Yole-walker方程法<快速遞推法<Burg算法。基于AR模型的現(xiàn)代功率譜估計質量明顯優(yōu)于經(jīng)典譜估計，基于AR模型的現(xiàn)代譜估計是逐步改進，性能遞增，精度提高的，直接解Yole-Walker方程運算量大，且不一定有解，采用Levinson遞推可以大大減小運算量，但由于自相關函數(shù)的估計默認取樣點以外的值為零，引入了一定的誤差，而Burg算法引入前后向預測誤差，來估計功率譜，避免了自相關函數(shù)帶來的誤差，因而精度較高，因此，Burg算法求得的AR模型最穩(wěn)定，而且Burg算法不需要自相關函數(shù)，所以性能優(yōu)于自相關法，由圖可知，在P=20,N=55時即可分辨出兩個峰值，所以在短數(shù)據(jù)時，Burg算法優(yōu)勢明顯，具有較高但是Burg算法進行譜估計會出現(xiàn)譜線分裂、譜峰偏移等問題，而且對于低信噪比情況，用AR模型很難準確估計出淹沒在噪聲中的正弦波頻率，而改進的MUSIC算法采用特征分解技術，在估計正弦信號和噪聲疊加的信號中明顯優(yōu)于AR模型功率譜估計。6信噪比與均方誤差的關系RR起初在剛接觸這門課時，甚至有些反感老師的嚴格要求，大概是習慣了其他老師的“溫柔鄉(xiāng)”吧，現(xiàn)在回過頭，在開始寫這篇報告時，對老師有一種說不出的感覺，我想那就是感謝的羞澀表達吧。雖然老師的表達有時不靠譜，但老師比那些看似靠譜的老師要強很多，他牢記灌輸知識的使命，又深諳人心之道，不時來點心靈雞湯，努力的鞭撻我們前進，將知識灌輸給我們，這讓我不由得想起了：“小時候，不愛吃飯的我被媽媽訓斥著吃飯的事”,想到這，對老師有一種莫名的親切，簡直深得我心。天啊，我也不知道要寫這些，也許，這tm是發(fā)自肺腑，不由而衷吧。首先，讓我對Matlab和Matlab的語法有了更深的了解，提高了我的編程能力，雖然以前開設過Matlab的實驗課，卻從沒有真正入門Matlab,可以說，陳老師是我Matlab的啟蒙老師啊!然后，讓我真正理解了一部分信號處理信號分析的知識，加深了對概念的理解，以前學信號分析信號處理時，只是記憶公式和概念，很少去分析為什么,甚至也不去管為什么,陳老師是一位真正有學識的人，他不照本宣科，總能說出自己的理解，而這些很多是其他老師未曾提到的(或許是我沒聽到其他老師提),這對于融會貫通整個學科至關重要。最后，讓我對一些問題有了更深層次的看法。程序附錄1經(jīng)典法(周期圖法)N=130;w=n/N;wn=randn(1,N)xn=10*sin(2*pi*0.2*n+(pi/3))+5*sin(2*pi*0.213*n+(pi/4))+wnXk=fft(xn);'兩個正弦信號與白噪聲疊加的時域波形);f=(locs-1)/(2*N)2直接解Yole-Walker方程N=256;p1=40;x=10*sin(2*pi*f1*n+pi/3)+5*sin(2*pi*f2*n+pi/4)+randn(1,N);sum=sum+x(n)*x(n);sum=0;sum=x(n)*x(n+m)+sum;sum=0;%求自相關函數(shù)form=1:p1a(i,k)=R(i-k);end%a為p行p列的矩陣，對角線及以上部分為0B(i,i)=RO,G=(RO+b)^(1/2);f=(locs-1)/(2*2000)3Levinson-Durbin快速遞推法N=256;n=1:N;p1=40;x=10*sin(2*pi*f1*n+pi/3)+5*sin(2*pi*f2*n+pi/4)+randn(1,N);'兩個正弦信號與白噪聲疊加的時域波形);%產(chǎn)生兩個正弦信號與白噪聲的sum=sum+x(n)*x(n);sum=0;sum=x(n)*x(n+m)+sum;sum=0;%求自相關函數(shù)form=2:p1forj=1:m-1sum=sum+aa(m-1,i)*R(m-i);aa(m,m)=-(R(m)+sum)/p(m-1);%定義levinson-durbin遞推式中的km%由H(z)用freqz()求解功率譜G=p(p1)^(1/2);A=[1,aa(p1;:)];[H,w]=freqz(B,A,2000);%將pi分成2000份f=(locs-1)/(2*2000)4Brug算法xn=10*sin(2*pi*n*f1+pi/3)+5*sin(2*pi*n*f2+pi/4)wn=randn(1,N);%產(chǎn)生高斯白噪聲中數(shù)組索引從1開始forj=1:Nsum=sum+xn(i)*xn(i);r(1)=sum/N;%因matlab中數(shù)組索引從1開始suma=0;sumb=0;suma=suma+ef(m-1,n)*eb(m-1,n-1);%初步構造KM中的分子sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2;%構造KM中的分母k(m-1)=-2*suma/sumb;%完全構造KM,根據(jù)前后向預測均方誤差之和最小來求反射系suma=0;sumb=0;%每次循環(huán)前賦零值forj=1:m-1aa(m,m)=k(m);p(m)=p(m-1)*(1-k(m)G=p(P)^(1/2);A=[1,aa(P;)];title('Burg算法功率譜)f=(locs-1)/(2*2000)x=10*sin(2*pi*0.2*n+pi/3)+5*sum=sum+x(n)*x(n);end%求自相關函數(shù)Rx(m,n)=R(m-n);%將值賦值給矩陣Rx(m,m)=R0;%將主對角線的值變?yōu)镽(0)列是對應的右特征向量，使得A*V=V*Dp=0;%特偵知的個數(shù)pev=e*V(:1:N-p);title('MUSIC算法求功率譜)xlabel('頻率)axis([00.5010])%e矩陣%功率譜信噪比與均方誤差的關系將求功率譜的各種算法改寫成函數(shù)以便直接調用，由與采用不同功率譜估計方法時，僅需替換代碼中的一部分，關系圖譜已附，以下僅提供采用Burg時的代碼和子程序代碼：MSE1=0;MSE2=0;V=26;s=[1:26];c=-15;10;M1=zeros(1,v);%定義1行V列的全零矩陣forc=1:v%求30次信噪比和均方誤差fori=1:500%循環(huán)500次求均方誤差x=awgn(x1,c,measured');%把信噪比為c的高斯白噪聲加到信號x1中pw1=Q(x);%我將burg算法改寫成了Q函數(shù)，此處直接調用。對應的位置f1=(locs(1,1)-1)/50f2=(locs(1,2)-1)/500;M1(1,s)=MSE1;M2(1,s)=MSE2S(1,s)=c;xlabel('SNR/dB');xlabel('SNR/dB');Q函數(shù)(burg子程序)的代碼：N=512;%定義采樣點數(shù)P=40;%濾波器階數(shù)的最大取值x1=10*sin(2*pi*n*f1+pi/3)+5*sin(2*pi*n*f2+pi/4);x=awgn(x1,10,'measured');%把信噪比為c的高斯白噪聲加到信號x1中ef(1,i)=x(i);中數(shù)組索引從1開始fori=1:Nr(1)=sum/N;%因matlab中數(shù)組索引從1開始%Burg遞推suma=0;sumb=0form=2:(P+1);%循環(huán)的階次，因matlab中數(shù)組索引從1開始suma=suma+ef(m-1,n)*eb(m-1,n-1);%初步構造KM中的分子sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2;%構造KM中的分母k(m-1)=-2*suma/sumb;%完全構造KM,根據(jù)前后向預測均方誤差之和最小來求反射系suma=0;sumb=0;%每次循環(huán)前賦零值forn