自己編寫算法的功率譜密度的三種matlab實(shí)現(xiàn)方法

1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)功率譜密度的三種matlab實(shí)現(xiàn)方法一：實(shí)驗(yàn)?zāi)康模?1)掌握三種算法的概念、應(yīng)用及特點(diǎn)；(2)了解譜估計(jì)在信號分析中的作用；(3) 能夠利用burg法對信號作譜估計(jì)，對信號的特點(diǎn)加以分析。二;實(shí)驗(yàn)內(nèi)容：簡單說明三種方法的原理。用三種方法編寫程序，在matlab中實(shí)現(xiàn)。將計(jì)算結(jié)果表示成圖形的形式，給出三種情況的功率譜圖。比較三種方法的特性。寫出自己的心得體會(huì)。三：實(shí)驗(yàn)原理：1.周期圖法：周期圖法又稱直接法。它是從隨機(jī)信號x(n)中截取N長的一段，把它視為能量有限x(n)真

2、實(shí)功率譜的估計(jì)的抽樣.認(rèn)為隨機(jī)序列是廣義平穩(wěn)且各態(tài)遍歷的，可以用其一個(gè)樣本x(n)中的一段來估計(jì)該隨機(jī)序列的功率譜。這當(dāng)然必然帶來誤差。由于對采用DFT，就默認(rèn)在時(shí)域是周期的，以及在頻域是周期的。這種方法把隨機(jī)序列樣本x(n)看成是截得一段的周期延拓，這也就是周期圖法這個(gè)名字的來歷。相關(guān)法（間接法):這種方法以相關(guān)函數(shù)為媒介來計(jì)算功率譜，所以又叫間接法。這種方法的具體步驟是：第一步：從無限長隨機(jī)序列x(n)中截取長度N的有限長序列列第二步：由N長序列求（2M-1）點(diǎn)的自相關(guān)函數(shù)序列。 (2-1) 這里，m=-(M-1),-1,0,1,M-1，MN，是雙邊序列，但是由自相關(guān)函數(shù)的偶對稱性式，只要

3、求出m=0,。,M-1的傅里葉變換，另一半也就知道了。第三步：由相關(guān)函數(shù)的傅式變換求功率譜。即以上過程中經(jīng)歷了兩次截?cái)?，一次是將x(n)截成N長，稱為加數(shù)據(jù)窗，一次是將x(n)截成（2M-1）長，稱為加延遲窗。因此所得的功率譜僅是近似值，也叫譜估計(jì)，式中的代表估值。一般取MN，因?yàn)橹挥挟?dāng)M較小時(shí)，序列傅式變換的點(diǎn)數(shù)才較小，功率譜的計(jì)算量才不至于大到難以實(shí)現(xiàn)，而且譜估計(jì)質(zhì)量也較好。因此，在FFT問世之前，相關(guān)法是最常用的譜估計(jì)方法。三：Burg法：AR模型功率譜估計(jì)又稱為自回歸模型，它是一個(gè)全極點(diǎn)的模型，要利用AR模型進(jìn)行功率譜估計(jì)須通過levinson_dubin遞推算法由Yule-Walke

4、r方程求得AR的參數(shù):2，12p。計(jì)算中，預(yù)測系數(shù)必須滿足Lenvinson-Durbin遞推關(guān)系，并且可直接計(jì)算而無需首先計(jì)算自相關(guān)系數(shù)。這種方法的優(yōu)點(diǎn)就是對未知數(shù)據(jù)不需要做任何假設(shè)，估計(jì)精度較高。其缺點(diǎn)是在分析噪聲中的正弦信號時(shí)，會(huì)引起譜線分裂，且譜峰的位置和正弦信號的相位有很大的關(guān)系。Burg算法是使前向預(yù)測誤差和后向預(yù)測誤差均方誤差之和最小來求取Km的，它不對已知數(shù)據(jù)段之外的數(shù)據(jù)做認(rèn)為假設(shè)。計(jì)算m階預(yù)測誤差的遞推表示公式如下：求取反射系數(shù)的公式如下：對于平穩(wěn)隨機(jī)過程，可以用時(shí)間平均代替集合平均，因此上式可寫成：這樣便可求得AR模型的反射系數(shù)。將m階AR模型的反射系數(shù)和m-1階AR模型的

5、系數(shù)代入到Levinson關(guān)系式中，可以求得AR模型其他的p-1個(gè)參數(shù)。Levinson關(guān)系式如下：m階AR模型的第m+1個(gè)參數(shù)G，其中是預(yù)測誤差功率，可由遞推公式求得。易知為進(jìn)行該式的遞推，必須知道0階AR模型誤差功率，可知該式由給定序列易于求得。完成上述過程，即最終求得了表征該隨機(jī)信號的AR模型的p+1個(gè)參數(shù) 。然后根據(jù) 即可求得該隨機(jī)信號的功率譜密度。實(shí)驗(yàn)內(nèi)容：實(shí)驗(yàn)程序及實(shí)驗(yàn)圖像周期法：Fs=1000;nfft=10000; %2nn=0:Fs;x=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n);X=fft(x,nfft);

6、Pxx=abs(X).2/length(n); %求解PSDt=0:round(nfft/2-1);f=t/nfft;P=10*log10(Pxx(t+1); %縱坐標(biāo)的單位為dBplot(f,P);grid on nfft=200 nfft=1024nfft=10000相關(guān)法：clear;Fs=1000; %采樣頻率n=0:Fs;%產(chǎn)生含有噪聲的序列nfft=1024;xn=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n);cxn=xcorr(xn,unbiased); %計(jì)算序列的自相關(guān)函數(shù)CXk=fft(cxn,nfft);

7、%求出功率譜密度Pxx=abs(CXk);index=0:round(nfft/2-1);f=index/nfft;plot_Pxx=10*log10(Pxx(index+1);plot(f,plot_Pxx);xlabel(頻率);ylabel(功率/DB);grid on;nfft=256nfft=1024Burg法：clearFs=1000 %設(shè)置關(guān)鍵變量，可通過調(diào)節(jié)這些變量觀察不同效果f1=0.2;f2=0.213;nfft=128; %取樣點(diǎn)數(shù)p=50; %階數(shù)p應(yīng)該選擇在N/3pN/2delta=1;m=sqrt(-1);f=0:1/1000:0.5;n=1:Fs;xn=sin(2

8、*pi*f1*n)+sqrt(2)*sin(2*pi*f2*n)+randn(size(n);figure;plot(n,xn);title(burg時(shí)域);xn= xn(:);N=length(xn);ef = xn;eb = xn;a = 1;for l=1:p efp = ef(2:end);%m-1階前向預(yù)測誤差 ebp = eb(1:end-1);%m-1階后向預(yù)測誤差 num = -2.*ebp*efp;%1km分子多項(xiàng)式 den = efp*efp+ebp*ebp;%1km的分母多項(xiàng)式 k(l) = num ./ den;%計(jì)算反射系數(shù) % 更新前向和后向預(yù)測誤差 ef = ef

9、p + k(l)*ebp;%各階前向預(yù)測誤差 eb = ebp + k(l)*efp;%各階后向預(yù)測誤差 % 計(jì)算模型參數(shù) a=a;0 + k(l)*0;conj(flipud(a);%AR模型參數(shù)aenda1=a(2:p+1); for i=1:length(f) %循環(huán)遞推 sum=0; for k=1:p sum=sum+a1(k)*exp(-m*2*pi*f(i)*k); end Pbrg(i)=delta/(abs(1+sum)2; Pbrg_f(i)=10*log10(Pbrg(i);%求出功率譜endfigureplot(f,Pbrg_f);title(burg頻域); nfft

10、=128 nfft=256五：結(jié)果比較分析(1) 在采樣點(diǎn)相同的時(shí)，周期圖法的特點(diǎn)是離散性大，曲線粗糙，方差較大，但是分辨率較高；采用周期突發(fā)估計(jì)得出的功率譜很不平滑，相應(yīng)的估計(jì)協(xié)方差比較大。而且采用增加采樣點(diǎn)的辦法也不能吃周期圖變得更加平滑，這是周期圖法的缺點(diǎn)。周期圖法得出的估計(jì)譜方差特性不好：當(dāng)數(shù)據(jù)長度N太大時(shí)，撲線的起伏加??；N太小時(shí)譜的分辨率又不好。對其改進(jìn)的主要方法有二種，即平均和平滑，平均就是將截取的數(shù)據(jù)段再分成L個(gè)小段，分別計(jì)算功率譜后取功率譜的平均，這種方法使估計(jì)的方差減少，但偏差加大，分辨率下降。平滑是用一個(gè)適當(dāng)?shù)拇昂瘮?shù)與算出的功率譜進(jìn)行卷積，使譜線平滑。這種方法得出的譜估計(jì)

11、是無偏的，方差也小，但分辨率下降。(2)相關(guān)法當(dāng)延遲與數(shù)據(jù)長度之比很小時(shí)，可以有良好的估計(jì)精度，相關(guān)法的收斂性較好，曲線平滑，方差較小，但是功率譜主瓣較寬，分辨率低。 （3）用Burg算法進(jìn)行功率譜估計(jì)時(shí)令前后向預(yù)測誤差功率之和最小,即對前向序列誤差和后向序列誤差前后都不加窗。Burg算法是建立在數(shù)據(jù)基礎(chǔ)之上的，避免了先計(jì)算自相關(guān)函數(shù)從而提高計(jì)算速度。是較為通用的方法，計(jì)算不太復(fù)雜 并且分辨率優(yōu)于自相關(guān)法。但對于白噪聲加正弦信號有時(shí)會(huì)出現(xiàn)譜線分裂現(xiàn)象，并且從上兩個(gè)圖中可以看出burg法產(chǎn)生的功率譜曲線比較平滑即方差小，分辨率高，可以明顯的觀察到兩個(gè)譜峰，在降低模型階次后譜的分辨率降低(兩個(gè)譜峰

12、幾乎變成一個(gè)譜峰)，但是曲線的平滑性更好。并且采樣點(diǎn)數(shù)越大，譜圖的分辨率就越高。對比nfft=128和nfft=256即可發(fā)現(xiàn)。除此之外還發(fā)現(xiàn)對于上面三種情況采樣點(diǎn)數(shù)越大，其功率譜密度也越大。還有就是階數(shù)p應(yīng)該選擇在N/3pN/2比較合適，這個(gè)通過測試可以得到驗(yàn)證。六：心得體會(huì) 第一次作業(yè)在上課的時(shí)候被老師點(diǎn)到，當(dāng)時(shí)老師問到burg法產(chǎn)生的圖像是否正確，當(dāng)時(shí)我覺得應(yīng)該沒錯(cuò)誤啊，不過因?yàn)樽约簩φ麄€(gè)過程的理解有限，沒聽出老師要表達(dá)什么意思，所以就只能沉默了。不過這篇實(shí)驗(yàn)報(bào)告確實(shí)不是復(fù)制粘貼得到的，是經(jīng)過查詢很多資料自己寫出來的。后來下課后自己看了下程序，發(fā)現(xiàn)老師說的應(yīng)該是burg法產(chǎn)生的圖像沒有兩個(gè)信號的譜峰，這是因?yàn)閄軸沒有進(jìn)行單位的調(diào)節(jié)導(dǎo)致譜圖在最左側(cè)位置所以看不到。并且自己上次是直接調(diào)用matlab中的Pburg函數(shù)，這雖然簡單省事但是對burg算法無法進(jìn)行較為深入的認(rèn)識(shí)，并且后來問了同學(xué)才知道老師是不允許直接調(diào)用burg函數(shù)的。雖然上課時(shí)被老師點(diǎn)到，出了點(diǎn)小丑，不過我覺得這倒是對自己的一個(gè)激勵(lì)，否則的話我到以為自己做的還行，實(shí)際上自己做的根本沒有達(dá)到老師的要求，在下課后自己為了用burg算法進(jìn)行求解功率譜自己又在網(wǎng)上查找了很多資料并結(jié)合課本上的介紹，網(wǎng)上很多資料