數(shù)字信號處理實驗-譜估計

實用文檔實驗四短時傅里葉變換實驗要求信號是正弦波加正態(tài)零均值白噪聲，信噪比為10dB，信號頻率為2kHz,取樣頻率為100kHz,數(shù)據(jù)長度N=256。用周期圖法進行譜估計；采用漢明窗，分段長度L=32，用修正的周期圖求平均法進行譜估計；并分析數(shù)據(jù)長度N，分段長度M’對譜估計結(jié)果的影響。實驗原理1、周期圖法功率譜的另一定義式為：忽略上式中求統(tǒng)計平均的運算，觀測數(shù)據(jù)為：,得到周期圖法的定義：2、修正的周期圖求平均法（1）把數(shù)據(jù)長度為N的信號分成L段，每段數(shù)據(jù)長度為M，N=LM。然后把窗函數(shù)加到每一個數(shù)據(jù)段上，求出每一段的周期圖，形成修正的周期圖，再對每一個修正的周期圖進行平均。第i段的修正周期圖為式中將每一段的修正周期圖之間近似看成互不相關(guān)，最后功率譜估計為實驗結(jié)果分析與討論1、周期圖法進行譜估計仿真結(jié)果：2、用Welch算法作為修正的周期圖求平均法對信號進行譜估計（1）數(shù)據(jù)長度N對譜估計的影響，保持分段長度M不變，分別取N=256,512,1024進行譜估計由仿真結(jié)果可知，數(shù)據(jù)長度N越大，功率譜的估計對真實譜的還原度越好，不過變化程度不是非常明顯。（2）分段長度M對譜估計結(jié)果的影響，保持數(shù)據(jù)長度N不變，分別取M=32,64,128進行譜估計由仿真結(jié)果可知，分段長度M越小，功率譜的估計對真實譜的還原度越好，而且變化程度比較大；同時隨著M的減小，頻率分辨率降低了。源程序1、周期圖法N=256;%采樣點數(shù)為256

n=0:N-1;fs=100000;%采樣頻率為100kHz

f0=2000;%輸入信號的頻率為2kHz

x1=sin(2*pi*f0/fs*n);%正弦信號

x=awgn(x1,10);%加入10dB信噪比的高斯白噪聲

window=boxcar(length(x));%矩形窗[p,f]=periodogram(x,window,256,fs)%用periodogram函數(shù)計算功率譜

%FFT計算長度256，以fs的抽樣頻率繪圖

subplot(2,2,1);plot(n,x);title('時域波形');subplot(2,2,2);plot(f,p);title('功率譜曲線');p1=10*log10(p);subplot(2,2,3);plot(f,p1);title('伯德圖形式功率譜圖')2、數(shù)據(jù)長度N對譜估計的影響fs=100000;%采樣頻率為100kHzf0=2000;%輸入信號的頻率為2kHzM=32;%選定分段長度M=32noverlap=10;%子樣本交疊長度10%N=256時的功率譜估計N1=256;%采樣點數(shù)為256n=0:N1-1;x1=sin(2*pi*f0/fs*n);%正弦信號x=awgn(x1,10);%加入10dB信噪比的高斯白噪聲window=hamming(32);%漢明窗[p1,f1]=pwelch(x,window,noverlap,N1,fs);p11=10*log10(p1);%N=512時的功率譜估計N2=512;%采樣點數(shù)為512n=0:N2-1;x1=sin(2*pi*f0/fs*n);x=awgn(x1,10);[p2,f2]=pwelch(x,window,noverlap,N2,fs);p22=10*log10(p2);%N=1024時的功率譜估計N3=1024;%采樣點數(shù)為1024n=0:N3-1;x1=sin(2*pi*f0/fs*n);x=awgn(x1,10);[p3,f3]=pwelch(x,window,noverlap,N3,fs);p33=10*log10(p3);plot(f1,p11,'-',f2,p22,'--',f3,p33,'-.')legend('N=256','N=512','N=1024',4)title('分段長度M不變，數(shù)據(jù)長度N變化')3、分段長度M對譜估計結(jié)果的影響N=256;%采樣點數(shù)為256n=0:N-1;fs=100000;%采樣頻率為100kHzf0=2000;%輸入信號的頻率為2kHzx1=sin(2*pi*f0/fs*n);%正弦信號x=awgn(x1,10);%加入10dB信噪比的高斯白噪聲M1=32;%分段長度32M2=64;%分段長度64M3=128;%分段長度128noverlap=10;%交疊長度10window1=hamming(M1);window2=hamming(M2);window3=hamming(M3);[p1,f1]=pwelch(x,window1,noverlap,256,fs);p11=10*log10(p1);[p2,f2]=pwelch(x,window2,noverlap,256,fs);p22=10*log10(p2);[p3,f3]=pwelch(x,window3,noverlap,256,fs);p33=10*l