譜估計(jì)實(shí)驗(yàn)報告

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

2、、協(xié)方差法、修正協(xié)方差法估計(jì)功率譜，分析數(shù)據(jù)長度和模型階次對估計(jì)結(jié)果的影響（可采用MATLAB自帶的功率譜分析函數(shù)）。（2）、調(diào)整正弦信號信噪比，分析信噪比的降低對估計(jì)效果的影響。報告內(nèi)容一、實(shí)驗(yàn)題目一1、問題分析（1）、w(n)與x(n)的產(chǎn)生w(n)產(chǎn)生：均值為0，方差為1白噪聲利用matlab中randn函數(shù)即可。表達(dá)如下：w=sqrt(1)*randn(1,N); sqrt(1)表示方差為1。x(n)產(chǎn)生：第一種思路：利用迭代的方法由，其中，然后利用上述公式依次向后遞推即可得。matlab代碼實(shí)現(xiàn)如下，注意到matlab中元素下標(biāo)都是從1開始的：x=w(1) zeros(1,N-1);

3、for i=2:N x(i)=0.8*x(i-1)+w(i);end此方法簡單，可以很容易地產(chǎn)生所需數(shù)目的數(shù)據(jù)。第二種思路：利用卷積的方法對線性時不變系統(tǒng)，輸入輸出滿足卷積關(guān)系：。由，可得，從而可得系統(tǒng)的沖擊響應(yīng)：。然后進(jìn)行卷積運(yùn)算即可。Matlab代碼實(shí)現(xiàn)如下：n=1:N;h(n)=(0.8).(n-1); %要注意n-1不是n，因?yàn)闆_擊響應(yīng)是從0開始的y=conv(w(n),h(n)； %共2*N-1項(xiàng)，取前N項(xiàng)即可需注意：實(shí)際h(n)是從0開始的，matlab處理元素從下標(biāo)1開始，因此，公式中應(yīng)是n-1不是n。而且，計(jì)算完成后卷積結(jié)果是為2*N-1項(xiàng)，取前N項(xiàng)即可。兩種方法結(jié)果為方便觀察

4、，令N=5時，實(shí)驗(yàn)結(jié)果如下：x = 0.6232 1.2976 1.9790 0.5911 0.6849 y = 0.6232 1.2976 1.9790 0.5911 0.6849 0.3437 0.0131 -0.2978 0.0868取卷積的前N項(xiàng)，可以看出兩種方式結(jié)果是相同的。（2）、信號x(n)的理想功率譜系統(tǒng)為AR模型，理想功率譜為：因此，對h(n)進(jìn)行傅里葉變換后，取模的平方即可。（3）、周期圖法譜估計(jì)根據(jù)周期圖法譜估計(jì)原理：先對觀測數(shù)據(jù)x(n)進(jìn)行傅里葉變換，后平方，最后除以N即可。（4）、Bartlett平均周期圖譜估計(jì)為了減小估計(jì)的方差，將數(shù)據(jù)分為L段，則每段有M=N/L個

5、數(shù)據(jù)，分別用周期圖法進(jìn)行估計(jì)后求平均。具體公式如下：將得到的L個周期圖進(jìn)行平均，作為信號x(n)的功率譜估計(jì)，公式如下：經(jīng)過平均處理，可使估計(jì)方差減小。2、實(shí)驗(yàn)結(jié)果與分析（1）、w(n)與x(n)的產(chǎn)生圖1 白噪聲w(n)圖2 平穩(wěn)隨機(jī)信號x(n)（2）、信號x(n)的理想功率譜信號理想功率譜如下圖3所示：圖3 信號的理想功率譜從圖中可以看出，理想功率譜是平滑的。下圖4是功率譜的分貝形式：圖4 信號的理想功率譜（dB）（3）、周期圖法譜估計(jì)當(dāng)數(shù)據(jù)長度N為1024時，實(shí)驗(yàn)結(jié)果如下圖5所示：圖5 N=1024周期圖譜估計(jì)結(jié)果當(dāng)數(shù)據(jù)長度N為256時，實(shí)驗(yàn)結(jié)果如下圖6所示：圖6 N=256周期圖譜估計(jì)

6、結(jié)果對比圖如下圖7：圖7 N=1024/256周期圖譜估計(jì)對比從以上結(jié)果可以看出N=1024時頻譜分辨率明顯要高于N=256時的頻譜分辨率。（4）、Bartlett平均周期圖譜估計(jì)當(dāng)數(shù)據(jù)長度N=1024及256時，分段數(shù)L分別為2和8時信號的功率譜為便于對比，將結(jié)果表示如下圖8一幅圖中：圖8 Bartlett平均周期圖譜估計(jì)結(jié)果分析：1、首先數(shù)據(jù)的長度對分辨率有影響，數(shù)據(jù)長度N=1024時的頻譜分辨率比N=256時的頻譜分辨率高；2、分段數(shù)L對頻譜的分辨率和平滑性（方差）也有很大影響。當(dāng)數(shù)據(jù)數(shù)目N一定時，L加大，每一段的數(shù)據(jù)量M就會減小，因此估計(jì)方差減小，曲線越平滑，但此時偏移加大，分辨率降低

7、，即估計(jì)量的方差和分辨率是一對矛盾，二者的效果可以通過合理選取L互換。3、收獲體會1、通過實(shí)驗(yàn)，對經(jīng)典譜估計(jì)法有了更深刻的理解，根據(jù)課堂中經(jīng)典譜估計(jì)的理論，掌握了周期圖法，分段周期圖法的具體實(shí)現(xiàn)，更加深刻地體會到了理論的原理以及這些估計(jì)方法的優(yōu)缺點(diǎn)，對這些估計(jì)方法獲得了真正的理解；2、很小的細(xì)節(jié)也要注意，也需要認(rèn)真思考，如信號x(n)的產(chǎn)生，通過兩種方法的對比，對卷積有了更深入的認(rèn)識，通過不同方法得到相同結(jié)果，學(xué)會了從不同的角度分析問題；3、對matlab掌握更好，如自己利用子函數(shù)調(diào)用的方式實(shí)現(xiàn)了對不同分段平均周期圖法的實(shí)現(xiàn)。函數(shù)傳遞參數(shù)分段數(shù)L即數(shù)據(jù)x(n)，子函數(shù)代碼如下：function

8、 Px_L=Px_mean(L,x)N=length(x); %計(jì)算數(shù)據(jù)長度M=N/L; %每段有M個數(shù)據(jù)tmp=zeros(1,M); %初始化中間變量，用于保存每段觀測數(shù)據(jù)的功率 for i=1:L %觀每段有M=N/L個數(shù)據(jù) for j=1:M y(i,j)=x(j+i*M-M); %取出L段數(shù)據(jù) end xk=fft(y(i,:),M); %fft變換 Ix=(abs(xk).2)./M; %計(jì)算每段的功率 Px=tmp+Ix; %累加每段的功率 end Px_L=Px./L; %平均功率譜二、實(shí)驗(yàn)題目二1、問題分析（1）、信號產(chǎn)生首先根據(jù)題目給定的信噪比計(jì)算正弦信號的幅度，產(chǎn)生所需的

9、正弦信號后與噪聲信號疊加形成觀測信號。若正弦信號幅度為A，則其功率為：；信噪比：，其中P為噪聲功率；由以上公式便可計(jì)算出兩個正弦信號的幅度。（2）利用MATLAB功率譜分析自相關(guān)法首先由觀測數(shù)據(jù)估計(jì)出其自相關(guān)函數(shù)，再解該方程得到模型參數(shù)，進(jìn)而求得信號的功率譜，一般利用Levenson-Durbin遞推法求解。matlab中的自相關(guān)法譜估計(jì)如下： pyulear(xn,p,nfft,fs); 其中， xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點(diǎn)，fs是采樣頻率。Burg法Burg法直接利用觀測數(shù)據(jù)計(jì)算AR模型參數(shù)，適用于觀測數(shù)據(jù)長度較短的時候。matlab中Burg法譜估

10、計(jì)如下： pburg(xn,p,nfft,fs);其中，px1是返回的功率譜，f1是返回的頻率分布；xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點(diǎn)，fs是采樣頻率。協(xié)方差法協(xié)方差法利用到的數(shù)據(jù)完全一致，相對自相關(guān)法沒有數(shù)據(jù)為0的假設(shè)。matlab中協(xié)方差法譜估計(jì)如下： pcov(xn,p,nfft,fs);其中，px1是返回的功率譜，f1是返回的頻率分布；xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點(diǎn)，fs是采樣頻率。修正協(xié)方差法修正協(xié)方差法使用到了前向和后向預(yù)測誤差對模型參數(shù)進(jìn)行估計(jì)。matlab中的修正協(xié)方差法譜估計(jì)如下： pmcov(xn,p,

11、nfft,fs);其中，px1是返回的功率譜，f1是返回的頻率分布；xn是輸入觀測信號，p是階數(shù)，nfft是觀測數(shù)據(jù)兩倍長度的頻率采樣點(diǎn)，fs是采樣頻率。2、實(shí)驗(yàn)結(jié)果與分析（1）數(shù)據(jù)長度的影響：將四種方法產(chǎn)生的功率譜顯示在了一幅圖中。L=1024時，如下圖9所示。L=256時，結(jié)果如下圖10所示。 圖9 L=1024圖10 L=256分析：對比兩幅圖可以發(fā)現(xiàn)：隨著數(shù)據(jù)長度的減小功率譜的分辨率在逐漸降低。即數(shù)據(jù)長度越大，功率譜分辨率越高，數(shù)據(jù)長度越小，功率譜分辨率越低。（2）模型階次的影響：當(dāng)階次發(fā)生變化時，如p=100時結(jié)果如圖11所示，p=30時，結(jié)果如圖12所示。圖11 階次p=100圖1

12、2 階次p=30分析：對比階次p=100與階次p=30時的結(jié)果可以看出，當(dāng)階次太小時，分辨率比較低，無法估計(jì)出小信噪比的信號，即小信噪比信號無法在圖中觀察出來。（3）信噪比的影響：將110Hz正弦信號信噪比降低到10，結(jié)果如下圖13所示：圖13 SNR1=10 SNR2=10對比發(fā)現(xiàn)：變化最明顯的是自相關(guān)法和Burg法，當(dāng)兩個正弦信號信噪比相差不多時，自相關(guān)法的估計(jì)效果還是較好的；而Burg法效果很差，每次實(shí)驗(yàn)的結(jié)果不同，經(jīng)常分辨不出一個正弦頻率，出現(xiàn)了譜峰偏移和譜線分裂的現(xiàn)象。而對于自相關(guān)法，并不適合于分辨信噪比相差較多的兩個正弦信號混合的情況，結(jié)果會譜峰偏移，出現(xiàn)觀測不到一個正弦信號的頻率

13、。如修改信噪比SNR1=30 SNR2=10的情況，結(jié)果如下圖14所示，仍然分辨不出信噪比較低的信號，即110Hz ，SNR2=10的正弦信號。圖13 SNR1=30 SNR2=103、收獲體會1、通過實(shí)驗(yàn)，加深了對AR模型譜估計(jì)的流程，掌握了結(jié)合matlab利用自相關(guān)法、Burg法、協(xié)方差法、修正協(xié)方差法估計(jì)功率譜，體會了各種方法的性能，優(yōu)缺點(diǎn)等；2、能更加熟練地通過matlab幫助實(shí)現(xiàn)對函數(shù)的查找，使用；3、由于考試最近考試較多，時間有限，沒有自己編程實(shí)現(xiàn)各種現(xiàn)代功率譜估計(jì)的方法，在時間空閑時，也會結(jié)合課堂講授的流程圖對各種方法嘗試自行編程實(shí)現(xiàn)。附 matlab源代碼實(shí)驗(yàn)題目一clc; %

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

15、conv(w(n),h(n) %共2*N-1項(xiàng)，取前N項(xiàng)即可% 令N=5時，實(shí)驗(yàn)結(jié)果，取卷積的前N項(xiàng)，可以看出兩種方式結(jié)果是相同的% x =% 0.6232 1.2976 1.9790 0.5911 0.6849% y =% 0.6232 1.2976 1.9790 0.5911 0.6849 0.3437 0.0131% -0.2978 0.0868% 當(dāng)數(shù)據(jù)長度為256時，取前256個數(shù)據(jù)為x2N2=256;x2=x(1:256);% 理想功率譜計(jì)算n=1:N;h(n)=(0.8).(n-1); %要注意n-1不是n，因?yàn)闆_擊響應(yīng)是從0開始的H=fft(h,N); %對沖擊響應(yīng)做傅里葉變換

16、Px_ideal=(abs(H).2;Px_ideal_db=10*log10(Px_ideal); %以分貝為單位表征% 周期圖法估計(jì)功率譜% 數(shù)據(jù)長度為1024時xk=fft(x,N);%傅里葉變換Px_zqt=(abs(xk).2)./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=2

17、x(1024個數(shù)據(jù))時Px_1024_8=Px_mean(8,x); %參數(shù)L=8 x(1024個數(shù)據(jù))時Px_256_2=Px_mean(2,x2); %參數(shù)L=2 x(256個數(shù)據(jù))時Px_256_8=Px_mean(8,x2); %參數(shù)L=8 x(256個數(shù)據(jù))時Px_1024_2_db=10*log10(Px_mean(2,x); %參數(shù)L=2 x(1024個數(shù)據(jù))時Px_1024_8_db=10*log10(Px_mean(8,x); %參數(shù)L=8 x(1024個數(shù)據(jù))時Px_256_2_db=10*log10(Px_mean(2,x2); %參數(shù)L=2 x(256個數(shù)據(jù))時Px_2

18、56_8_db=10*log10(Px_mean(8,x2); %參數(shù)L=8 x(256個數(shù)據(jù))時% 顯示噪聲w(n)與平穩(wěn)信號x(n)figure;plot(n,w(n); xlabel('n'); ylabel('w(n)');legend('隨機(jī)噪聲w(n)')figure;plot(n,x(n); xlabel('n'); ylabel('x(n)');legend('平穩(wěn)隨機(jī)信號x(n)')% 顯示理想功率譜figure;plot(f,Px_ideal);xlabel('頻率f（H

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

20、x');legend('1024');subplot(2,2,2);plot(f,Px_zqt_db);xlabel('頻率f（Hz）');ylabel('功率譜Px（dB）');legend('1024');% 數(shù)據(jù)長度為256時subplot(2,2,3);plot(f,Px_zqt2);xlabel('頻率f（Hz）');ylabel('功率譜Px');legend('256');subplot(2,2,4);plot(f,Px_zqt_db2);xlabel('

21、;頻率f（Hz）');ylabel('功率譜Px（dB）');legend('256');% 平均周期圖法結(jié)果figure;subplot(2,2,1);plot(1:length(Px_1024_2),Px_1024_2);xlabel('頻率f（Hz）');ylabel('功率譜Px');legend('L=2 N=1024');subplot(2,2,2);plot(1:length(Px_1024_8),Px_1024_8);xlabel('頻率f（Hz）');ylabel('

22、;功率譜Px');legend('L=8 N=1024');subplot(2,2,3);plot(1:length(Px_256_2),Px_256_2);xlabel('頻率f（Hz）');ylabel('功率譜Px');legend('L=2 N=256');subplot(2,2,4);plot(1:length(Px_256_8),Px_256_8);xlabel('頻率f（Hz）');ylabel('功率譜Px');legend('L=8 N=256');子程序P

23、x_mean.mfunction Px_L=Px_mean(L,x)N=length(x); %計(jì)算數(shù)據(jù)長度M=N/L; %每段有M個數(shù)據(jù)tmp=zeros(1,M); %初始化中間變量，用于保存每段觀測數(shù)據(jù)的功率 for i=1:L %觀每段有M=N/L個數(shù)據(jù) for j=1:M y(i,j)=x(j+i*M-M); %取出L段數(shù)據(jù) end xk=fft(y(i,:),M); %fft變換 Ix=(abs(xk).2)./M; %計(jì)算每段的功率 Px=tmp+Ix; %累加每段的功率 end Px_L=Px./L; %平均功率譜實(shí)驗(yàn)題目二clc;clear; fs=1000; % 采樣頻率 n=0:1/fs:2; %取樣數(shù)據(jù)點(diǎn) snr1=30; %正弦波d的信噪比snr2=10; %正弦波d2的信噪比dB A1=sqrt(2*10(snr1/10); %正弦波x1的幅度 A2=sqrt(2*10(snr2/10); %正弦波x2的幅度 An=sqrt(1); %噪聲w(n)的幅度 x1=A1*sin(2*pi*100*n); %正弦信號x1 x2=A2*sin(2*pi*110*n); %正弦信號x2 w=An*randn(size(n); %噪聲信號 w(n) xn=x1+x2+w;