多種功率譜估計(jì)鄭州大學(xué)隨機(jī)信號(hào)處理大作業(yè)

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