譜估計試驗報告匯總

1、第四章上機作業(yè)實驗報告實驗題目1、假設(shè)一平穩(wěn)隨機信號為x(n)=0.8x(n-1)+w(n),其中w(n)是均值為0,方差為1的白噪聲，數(shù)據(jù)長度為1024。(1)、產(chǎn)生符合要求的w(n)和x(n)；(2)、給出信號x(n)的理想功率譜；(3)、編寫周期圖譜估計函數(shù)，估計數(shù)據(jù)長度N=1024及256時信號功率譜，分析估計效果。(4)、編寫B(tài)artlett平均周期圖函數(shù)，估計當(dāng)數(shù)據(jù)長度N=1024及256時，分段數(shù)L分別為2和8時信號x(n)的功率譜，分析估計效果。2、假設(shè)均值為0,方差為1的白噪聲w(n)中混有兩個正弦信號，該正弦信號的頻率分別為100Hz和110Hz,信噪比分別為10dB和30

2、dB,初始相位都為0,采樣頻率為1000Hz。(1)、采用自相關(guān)法、Burg法、協(xié)方差法、修正協(xié)方差法估計功率譜，分析數(shù)據(jù)長度和模型階次對估計結(jié)果的影響(可采用MATLAB自帶的功率譜分析函數(shù))。(2)、調(diào)整正弦信號信噪比，分析信噪比的降低對估計效果的影響。報告內(nèi)容一、實驗題目一1、問題分析(1)、w(n)與x(n)的產(chǎn)生w(n)產(chǎn)生：均值為0,方差為1白噪聲w(n)利用matlab中randn函數(shù)即可。表達如下：w=sqrt(1)*randn(1,N);sqrt(1)表示方差為1。x(n)產(chǎn)生：第一種思路：利用迭代的方法由x(n)=0.8x(n-1)+w(n),其中x(0)=w(0),然后利

3、用上述公式依次向后遞推即可得x(n)。matlab代碼實現(xiàn)如下，注意到matlab中元素下標(biāo)都是從1開始的：x=w(1)zeros(1,N-1);fori=2:Nx(i)=0.8*x(i-1)+w(i);end此方法簡單，可以很容易地產(chǎn)生所需數(shù)目的數(shù)據(jù)。第二種思路：利用卷積的方法對線性時不變系統(tǒng)，輸入輸出滿足卷積關(guān)系：x(n)=w(n)*h(n)。1由x(n)=0.8xn-1+wn,可得H(z)=y,從而可得系統(tǒng)的沖擊響應(yīng)：1-0.8zh(n)=0.8nu(n)。然后進行卷積運算即可。Matlab代碼實現(xiàn)如下：n=1:N;h(n)=(0.8).A(n-1);%要注意n-1不是n,因為沖擊響應(yīng)是

4、從0開始的y=conv(w(n),h(n);%共2*N-1項，取前N項即可需注意：實際h(n)是從0開始的，matlab處理元素從下標(biāo)1開始，因此，公式中應(yīng)是n-1不是n。而且，計算完成后卷積結(jié)果是為2*N-1項，取前N項即可。兩種方法結(jié)果為方便觀察，令N=5時，實驗結(jié)果如下：x=0.62321.29761.97900.59110.6849y=0.62321.29761.97900.59110.68490.34370.0131-0.29780.0868取卷積的前N項，可以看出兩種方式結(jié)果是相同的。(2)、信號x(n)的理想功率譜一-.八.2.系統(tǒng)為AR模型，理想功率譜為：Pxx(ejw)=Q2

5、H(ejw)=H(ejw)因此，對h(n)進行傅里葉變換后，取模的平方即可。(3)、周期圖法譜估計根據(jù)周期圖法譜估計原理:A1NJLP/)=NFn)n先對觀測數(shù)據(jù)x(n)進行傅里葉變換，后平方，最后除以N即可。(4)、Bartlett平均周期圖譜估計為了減小估計的方差，將數(shù)據(jù)分為L段，則每段有M=N/L個數(shù)據(jù)，分別用周期圖法進行估計后求平均。具體公式如下：,、1Ii(w)二MM1x-Xi(n)e-jwnn=e將得到的L個周期圖進行平均，作為信號x(n)的功率譜估計，公式如下:jw、_1:Pxx(e)=-1i(w)Lt經(jīng)過平均處理，可使估計方差減小2、實驗結(jié)果與分析(1)、w(n)與x(n)的產(chǎn)

6、生圖1白噪聲w(n)6圖2平穩(wěn)隨機信號x(n)平穗隨機信號4n)1200(2)、信號x(n)的理想功率譜信號理想功率譜如下圖3所示：302520151050_22f020040060080010001200頻率f(Hz)圖3信號的理想功率譜從圖中可以看出，理想功率譜是平滑的下圖4是功率譜的分貝形式:-L181-16LI4，O42O8642O24B111-23qpBBP一圖4信號的理想功率譜(dB)(3)、周期圖法譜估計當(dāng)數(shù)據(jù)長度N為1024時，實驗結(jié)果如下圖5所示:圖5N=1024周期圖譜估計結(jié)果當(dāng)數(shù)據(jù)長度N為256時，實驗結(jié)果如下圖6所示:20200400600頻率f(H?)80010001

7、200001一xd軀Ms80604LI圖6N=256周期圖譜估計結(jié)果對比圖如下圖7：O1Er00Ijrow50頻從以上結(jié)果可以看出N=1024時頻譜分辨率明顯要高于N=256時的頻譜分辨率。(4)、Bartlett平均周期圖譜估計當(dāng)數(shù)據(jù)長度N=1024及256時，分段數(shù)L分別為2和8時信號x(n)的功率譜為便于對比，將結(jié)果表示如下圖8一幅圖中：60L=2N-1O24-Irl46420KdHS耕s00650:lrl50頻率f(He)050100150頻率f(He)0010203040432-xd軀科集頻率f(Hz)圖8Bartlett平均周期圖譜估計結(jié)果分析:1、首先數(shù)據(jù)的長度對分辨率有影響，數(shù)

8、據(jù)長度N=1024時的頻譜分辨率比N=256時的頻譜分辨率高；2、分段數(shù)L對頻譜的分辨率和平滑性(方差)也有很大影響。當(dāng)數(shù)據(jù)數(shù)目N一定時，L加大，每一段的數(shù)據(jù)量M就會減小，因此估計方差減小，曲線越平滑，但此時偏移加大，分辨率降低，即估計量的方差和分辨率是一對矛盾，二者的效果可以通過合理選取L互換。3、收獲體會1、通過實驗，對經(jīng)典譜估計法有了更深刻的理解，根據(jù)課堂中經(jīng)典譜估計的理論，掌握了周期圖法，分段周期圖法的具體實現(xiàn)，更加深刻地體會到了理論的原理以及這些估計方法的優(yōu)缺點，對這些估計方法獲得了真正的理解；2、很小的細節(jié)也要注意，也需要認真思考，如信號x(n)的產(chǎn)生，通過兩種方法的對比，對卷積有

9、了更深入的認識，通過不同方法得到相同結(jié)果，學(xué)會了從不同的角度分析問題；3、對matlab掌握更好，如自己利用子函數(shù)調(diào)用的方式實現(xiàn)了對不同分段平均周期圖法的實現(xiàn)。函數(shù)傳遞參數(shù)分段數(shù)L即數(shù)據(jù)x(n),子函數(shù)代碼如下：functionPx_L=Px_mean(L,x)N=length(x);%計算數(shù)據(jù)長度M=N/L;%每段有M個數(shù)據(jù)tmp=zeros(1,M);%初始化中間變量，用于保存每段觀測數(shù)據(jù)的功率fori=1:L%觀每段有M=N/L個數(shù)據(jù)forj=1:My(i,j)=x(j+i*M-M);%取出L段數(shù)據(jù)endxk=fft(y(i,:),M);%fft變換Ix=(abs(xk)A2)./M;%

10、計算每段的功率Px=tmp+Ix;%累加每段的功率endPx_L=Px./L;%平均功率譜二、實驗題目二1、問題分析(1)、信號產(chǎn)生首先根據(jù)題目給定的信噪比計算正弦信號的幅度，產(chǎn)生所需的正弦信號后與噪聲信號疊加形成觀測信號。A2若正弦信號幅度為A,則其功率為：S=A;2s信噪比：SNT=1Ologio一,其中P為噪聲功率；P由以上公式便可計算出兩個正弦信號的幅度。(2)利用MATLAB功率譜分析自相關(guān)法首先由觀測數(shù)據(jù)估計出其自相關(guān)函數(shù)，再解該方程得到模型參數(shù)，進而求得信號的功率譜，一般利用Levenson-Durbin遞推法求解。matlab中的自相關(guān)法譜估計如下：pyulear(xn,p,n

11、fft,fs);其中，xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點，fs是采樣頻率。Burg法Burg法直接利用觀測數(shù)據(jù)計算AR模型參數(shù)，適用于觀測數(shù)據(jù)長度較短的時候。matlab中Burg法譜估計如下：pburg(xn,p,nfft,fs);其中，px1是返回的功率譜，f1是返回的頻率分布；xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點，fs是采樣頻率。協(xié)方差法協(xié)方差法利用到的數(shù)據(jù)完全一致，相對自相關(guān)法沒有數(shù)據(jù)為0的假設(shè)。matlab中協(xié)方差法譜估計如下：pcov(xn,p,nfft,fs);其中，px1是返回的功率譜，f1是返回的頻率分布；xn

12、是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點，fs是采樣頻率。修正協(xié)方差法修正協(xié)方差法使用到了前向和后向預(yù)測誤差對模型參數(shù)進行估計。matlab中的修正協(xié)方差法譜估計如下：pmcov(xn,p,nfft,fs);其中，px1是返回的功率譜，f1是返回的頻率分布；xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點，fs是采樣頻率。2、實驗結(jié)果與分析(1)數(shù)據(jù)長度的影響：將四種方法產(chǎn)生的功率譜顯示在了一幅圖中。L=1024時,如下圖9所示。L=256時，結(jié)果如下圖10所示自相關(guān)法功率譜估計協(xié)方差法功率譜估計10sMd胭ss100120140頻率f<Hz)

13、口oo215is修正協(xié)方差法功率譜估計_05010015080864290100120140頻率f(Hz)burg法功率譜估計15100120頻率f(He)140105xats頻率f(Hz)L=1Q24圖9L=1024自相關(guān)法功率譜估計1lAr80100120140頻率f(Hz)修正協(xié)方差法功率譜估計50100150頻率f(Hz)協(xié)方差法功率譜估計,5登里博s080100120140頻率f(Hz)burg法功率譜估計8100120140頻率f(Hz)642Md出端存80L=256分析：對比兩幅圖可以發(fā)現(xiàn)：隨著數(shù)據(jù)長度的減小功率譜的分辨率在逐漸降低。即數(shù)據(jù)長度越大，功率譜分辨率越高，數(shù)據(jù)長度越小

14、，功率譜分辨率越低。(2)模型階次的影響：當(dāng)階次發(fā)生變化時，如p=100時結(jié)果如圖11所示，p=30時，結(jié)果如圖12所示。自相關(guān)法功率譜估計協(xié)方差法功率譜估計1500(£觸得叁80100120140頻!率f£Hz)burg法功率譜估計050100頻率f(Hz)八80100120140頻率f(Hz)修正協(xié)方差法功率譜估計15n胭隔離08010012U頻率f(Hz)140.5CL他搏吞p=100圖11階次p=100自相關(guān)法功率譜估計80100120140頻率"H?)修正協(xié)方差法功率譜估計5協(xié)方差法功率譜估計080100120140頻率ffH?)50100160頻率f(

15、Hz)001100120140頻率f(Hz)burg法功率譜估計oO5(14班解集BO自相關(guān)法功率譜估計00例100120140頻率f(Hz)修正協(xié)方差法功率譜估計上粗MSXCL觸新年協(xié)方差法功率譜估計460100150頻率f(Hz)32180100120140頻率f(Hz)00500分析：對比階次p=100與階次p=30時的結(jié)果可以看出，當(dāng)階次太小時，分辨率比較低，無法估計出小信噪比的信號，即小信噪比信號無法在圖中觀察出來。(3)信噪比的影響：將110Hz正弦信號信噪比降低到10,結(jié)果如下圖13所示:321、口100120140頻率f(Hz)bwg法功率譜估計圖13SNR1=10SNR2=1

16、0對比發(fā)現(xiàn)：變化最明顯的是自相關(guān)法和Burg法，當(dāng)兩個正弦信號信噪比相差不多時，自相關(guān)法的估計效果還是較好的；而Burg法效果很差，每次實驗的結(jié)果不同，經(jīng)常分辨不出一個正弦頻率，出現(xiàn)了譜峰偏移和譜線分裂的現(xiàn)象。而對于自相關(guān)法，并不適合于分辨信噪比相差較多的兩個正弦信號混合的情況，結(jié)果會譜峰偏移，出現(xiàn)觀測不到一個正弦信號的頻率。如修改信噪比SNR1=30SNR2=10的情況，結(jié)果如下圖14所示，仍然分辨不出信噪比較低的信號,即110Hz,SNR2=10的正弦信號。自相關(guān)法功率譜估計080100120140頻率f(Hz)修正協(xié)方差法功率譜估計ooOo-JO505協(xié)方差法功率謂估計080100120

17、140頻率f(Hz)burg法功率譜估計42Kd!S耕542xdssus642/疆»繇005010015080100120140頻率fHG頻率f(W)圖13SNR1=30SNR2=103、收獲體會1、通過實驗，加深了對AR模型譜估計的流程，掌握了結(jié)合matlab利用自相關(guān)法、Burg法、協(xié)方差法、修正協(xié)方差法估計功率譜，體會了各種方法的性能，優(yōu)缺點等；2、能更加熟練地通過matlab幫助實現(xiàn)對函數(shù)的查找，使用；3、由于考試最近考試較多，時間有限，沒有自己編程實現(xiàn)各種現(xiàn)代功率譜估計的方法，在時間空閑時，也會結(jié)合課堂講授的流程圖對各種方法嘗試自行編程實現(xiàn)。附matlab源代碼實驗題目一c

18、lc;%清除命令窗口clear;%清空變量N=1024;%數(shù)據(jù)長度df=1;%設(shè)定頻率分辨率fs=(N-1)*df;%計算采樣率t=0:1/fs:(N-1)/fs;%截取信號的時間段f=0:df:fs;%功率譜估計的頻率分辨率和范圍w=sqrt(1)*randn(1,N);%平穩(wěn)隨機信號的獲取%準(zhǔn)確獲得觀測信號是很關(guān)鍵的%處理方法一：迭代的方式x=wzeros(1,N-1);fori=2:Nx(i)=0.8*x(i-1)+w(i);end%處理方法二：卷積的方式%n=1:N;%h(n)=(0.8).A(n-1);%要注意n-1不是n,因為沖擊響應(yīng)是從0開始的%y=conv(w(n),h(n)%

19、共2*N-1項，取前N項即可%令N=5時，實驗結(jié)果，取卷積的前N項，可以看出兩種方式結(jié)果是相同的%x=%0.62321.29761.97900.59110.6849%y=%0.62321.29761.97900.59110.68490.34370.0131%-0.29780.0868%當(dāng)數(shù)據(jù)長度為256時，取前256個數(shù)據(jù)為x2N2=256;x2=x(1:256);%理想功率譜計算h(n)=(0.8).A(n-1);%要注意n-1不是n,因為沖擊響應(yīng)是從0開始的H=fft(h,N);%對沖擊響應(yīng)做傅里葉變換Px_ideal=(abs(H).A2;Px_ideal_db=10*log10(Px_

20、ideal);%以分貝為單位表征%周期圖法估計功率譜%數(shù)據(jù)長度為1024時xk=fft(x,N);%傅里葉變換Px_zqt=(abs(xk).A2)./N);Px_zqt_db=10*log10(Px_zqt);%數(shù)據(jù)長度為256時xk2=fft(x2,N);%傅里葉變換Px_zqt2=(abs(xk2)，2)./N2);Px_zqt_db2=10*log10(Px_zqt2);%Px_mean(L,x):平均周期圖法，分段數(shù)目L,數(shù)據(jù)為xPx_1024_2=Px_mean(2,x);%參數(shù)L=2x(1024個數(shù)據(jù))時Px_1024_8=Px_mean(8,x);%參數(shù)L=8x(1024個數(shù)據(jù)

21、)時Px_256_2=Px_mean(2,x2);%參數(shù)L=2x(256個數(shù)據(jù))時Px_256_8=Px_mean(8,x2);%參數(shù)L=8x(256個數(shù)據(jù))時Px_1024_2_db=10*log10(Px_mean(2,x);%參數(shù)L=2x(1024個數(shù)據(jù))時Px_1024_8_db=10*log10(Px_mean(8,x);%參數(shù)L=8x(1024個數(shù)據(jù))時Px_256_2_db=10*log10(Px_mean(2,x2);%參數(shù)L=2x(256個數(shù)據(jù))時Px_256_8_db=10*log10(Px_mean(8,x2);%參數(shù)L=8x(256個數(shù)據(jù))時%顯示噪聲w(n)與平穩(wěn)信號

22、x(n)figure;plot(n,w(n);xlabel('n');ylabel('w(n)');legend('隨機噪聲w(n),)figure;plot(n,x(n);xlabel('n');ylabel('x(n)');legend('平穩(wěn)隨機信號x(n),)%顯示理想功率譜figure;plot(f,Px_ideal);xlabel('頻率f(Hz)');ylabel('功率譜legend('理想功率譜Px_ideal');figure;plot(f,Px_idea

23、l_db);xlabel('頻率f(Hz)');ylabel('功率譜legend('理想功率譜Px_ideal_db');Px_ideal');Px_ideal_db(dB),)%周期圖法估計功率譜結(jié)果%數(shù)據(jù)長度為1024時figure;subplot(2,2,1);plot(f,Px_zqt);xlabel('頻率f(Hz)legend('1024');subplot(2,2,2);plot(f,Px_zqt_db);xlabel('頻率f(Hz)legend('1024');%數(shù)據(jù)長度為256

24、時subplot(2,2,3);plot(f,Px_zqt2);xlabel('頻率f(Hz)legend('256');');ylabel('功率譜');ylabel('功率譜');ylabel('功率譜Px');Px(dB)');Px');subplot(2,2,4);plot(f,Px_zqt_db2);xlabel('頻率f(Hz)');ylabel('功率譜Px(dB)');legend('256');%平均周期圖法結(jié)果figure;subp

25、lot(2,2,1);plot(1:length(Px_1024_2),Px_1024_2);xlabel('頻率f(Hz)');ylabel('功率譜Px');legend('L=2N=1024');subplot(2,2,2);plot(1:length(Px_1024_8),Px_1024_8);xlabel('頻率f(Hz)');ylabel('功率譜Px');legend('L=8N=1024');subplot(2,2,3);plot(1:length(Px_256_2),Px_256_

26、2);xlabel('頻率f(Hz)');ylabel('功率譜Px');legend('L=2N=256');subplot(2,2,4);plot(1:length(Px_256_8),Px_256_8);xlabel('頻率f(Hz)');ylabel('功率譜Px');legend('L=8N=256');子程序Px_mean.mfunctionPx_L=Px_mean(L,x)M=N/L;%每段有M個數(shù)據(jù)tmp=zeros(1,M);%初始化中間變量，用于保存每段觀測數(shù)據(jù)的功率fori=1:L%觀每段有M=N/L個數(shù)據(jù)forj=1:My(i,j)=x(j+i*M-M);%取出L段數(shù)據(jù)endxk=fft(y(i,:),M);%fft變換Ix=(abs(xk).A2)./M;%計算每段的功率Px=tmp+Ix;%累加每段的功率endPx_L=Px./L;%平均功率譜實驗題目二clc;clear;fs=1000;%采樣頻率n=0:1/fs:2;%取樣數(shù)據(jù)點snr1=30;%正弦波d的信噪比snr2=10;%正弦波d2的信噪比dBA1=sqrt(2*10A(snr1/10);A2=sqr