華工電信數(shù)字信號(hào)處理實(shí)驗(yàn)FFT算法應(yīng)用---頻率估計(jì)

1、數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)題目：FFT算法的應(yīng)用-頻率估計(jì)班級(jí)：2010級(jí)信息5班姓名：陳耿濤 學(xué)號(hào)：201030271709 序號(hào)：03 2012年11月11日一、實(shí)驗(yàn)題目1、對連續(xù)的單一頻率周期信號(hào) 按采樣頻率 采樣，截取長度N分別選N =20和N =16，觀察其DFT結(jié)果的幅度譜。 2、2N點(diǎn)實(shí)數(shù)序列 N=64。用一個(gè)64點(diǎn)的復(fù)數(shù)FFT程序，一次算出，并繪出。3、頻率估計(jì)1）產(chǎn)生一個(gè)單頻實(shí)信號(hào)，加上一定信噪比的噪聲。2）對含噪聲的信號(hào)進(jìn)行頻率估計(jì)。估計(jì)方法可以查閱相關(guān)文獻(xiàn)。3）統(tǒng)計(jì)估計(jì)出來的頻率和真實(shí)頻率之間的誤差。4）驗(yàn)證該頻率估計(jì)算法在不同信噪比、不同數(shù)據(jù)長度下、不同頻率時(shí)候的性能。二

2、、實(shí)驗(yàn)過程第一題：此時(shí)離散序列xa1=sin(3*pi*n1/k，即k=8。用MATLAB計(jì)算并作圖，函數(shù)fft用于計(jì)算離散傅里葉變換DFT。1、實(shí)驗(yàn)代碼：k=8;n1=0:1:19;xa1=sin(3*pi*n1/k;subplot(2,2,1;plot(n1,xa1;xlabel('t/T'ylabel('x(n'title('a'xk1=fft(xa1;xk1=abs(xk1;subplot(2,2,2;stem(n1,xk1;xlabel('k'ylabel('X(k'title('b:k=20&#

3、39;n2=0:1:15;xa2=sin(3*pi*n2/k;subplot(2,2,3;plot(n2,xa2;xlabel('t/T'ylabel('x(n'title('c'xk2=fft(xa2;xk2=abs(xk2;subplot(2,2,4;stem(n2,xk2;xlabel('k'ylabel('X(k'title('d:k=16'2、實(shí)驗(yàn)效果圖：圖 2-1 不同的截取長度的正弦信號(hào)及其DFT結(jié)果計(jì)算結(jié)果示于圖2-1，a和b分別是N=20時(shí)的截取信號(hào)和DFT結(jié)果，由于截取了兩個(gè)半周

4、期，頻譜出現(xiàn)泄漏；c和d分別是N=16時(shí)的截取信號(hào)和DFT結(jié)果，由于截取了兩個(gè)整周期，得到單一譜線的頻譜。上述頻譜的誤差主要是由于時(shí)域中對信號(hào)的非整周期截?cái)喈a(chǎn)生的頻譜泄漏。第二題：1、 實(shí)驗(yàn)代碼N=64;n=0:1:N-1; n1=2*n;n2=2*n+1;k=0:1:N-1;xn1=cos(2*pi/N*7*n1+1/2*cos(2*pi/N*19*n1;xn2=cos(2*pi/N*7*n2+1/2*cos(2*pi/N*19*n2; XK1=fft(xn1;XK2=fft(xn2; X1=XK1+exp(-j*pi*k/N.*XK2; X2=XK1-exp(-j*pi*k/N.*XK2;

5、 X1=X1 zeros(1,N;X2=zeros(1,N X2;XK=X1+X2; k=0:1:2*N-1;XK=abs(XK; stem(k,XK;xlabel('k'ylabel('|X(k|'title('X(k=DFTx(n2N'2、 實(shí)驗(yàn)結(jié)果圖 圖 2-1結(jié)果分析：由歐拉公式得：對，其2N點(diǎn)的DFT變換為：當(dāng)時(shí)， =0當(dāng)時(shí)，即由此可得當(dāng)k=14,38,90,114時(shí)有值其余為0）， 與圖2-1有相同的結(jié)論。第三題：1、算法分析在實(shí)驗(yàn)中選擇頻率f=50Hz的單頻信號(hào)，即x=sin(20*pi*t，利用函數(shù)awgn可以再信號(hào)中加入高斯白噪

6、聲，即y=awgn(x,SNR，其中信噪比的單位是分貝。對于含噪聲信號(hào)頻率的估計(jì)可以采用譜峰法，即作出y的傅里葉變換的圖象，其幅度最大的譜線對應(yīng)的頻點(diǎn)即為信號(hào)x的頻率。經(jīng)過驗(yàn)證，該算法在不同信噪比下性能均良好。2、實(shí)驗(yàn)代碼估計(jì)出含噪聲信號(hào)的頻率Fs = 1000; % 采樣頻率T = 1/Fs; % 采樣時(shí)間L = 1000; % 總的采樣點(diǎn)數(shù)fc=10;x=(0:1/fs:2;n=201;y1=sin(10*pi*fc*x; %產(chǎn)生一個(gè)信號(hào)的幅值為1頻率為50Hz的正弦信號(hào)a=0;b=0.5; %均值為a，方差為b2figure(1plot(x,y1;axis(0,0.2,-2,2;titl

7、e('未受干擾時(shí)信號(hào)的波形'ylabel('y'xlabel('x/20pi'grid;y2=awgn(y1,10; % 混入噪聲信號(hào),信噪比是10y=y1+y2; %加入高斯白噪聲之后的信號(hào)figure(2;plot(x,y;axis(0,0.2,-3,3;title('疊加了高斯白噪聲的信號(hào)波形'ylabel('y'xlabel('x/20pi'grid;NFFT = 2nextpow2(L; % 求得最接近總采樣點(diǎn)的2n,這里應(yīng)該是210=1024Y= fft(y,NFFT/L; %進(jìn)行fft

8、變換（除以總采樣點(diǎn)數(shù)，是為了后面精確看出原始信號(hào)幅值）f = Fs/2*linspace(0,1,NFFT/2+1;%頻率軸(只畫到Fs/2即可，由于y為實(shí)數(shù)，后面一半是對稱的% 畫出頻率幅度圖形，可以看出50Hz幅值是最高的figure(3;plot(f,2*abs(Y(1:NFFT/2+1 title('信號(hào)y(t的頻譜'xlabel('頻率 (Hz'ylabel('|Y(f|'統(tǒng)計(jì)估計(jì)出來的頻率和真實(shí)頻率之間的誤差for(j=1:1:10Fs = 1000; % 采樣頻率T = 1/Fs; % 采樣時(shí)間L = 1000; % 總的采樣點(diǎn)數(shù)t

9、 = (0:L-1*T; % 時(shí)間序列（時(shí)間軸）snr=10;%產(chǎn)生一個(gè)信號(hào)的幅值為1頻率為50Hz的正弦信號(hào)x = sin(2*pi*50*t;y = awgn(x,snr % 混入噪聲信號(hào)NFFT = 2nextpow2(L; % 求得最接近總采樣點(diǎn)的2n,這里應(yīng)該是210=1024Y = fft(y,NFFT/L; %進(jìn)行fft變換（除以總采樣點(diǎn)數(shù)，是為了后面精確看出原始信號(hào)幅值）f = Fs/2*linspace(0,1,NFFT/2+1;%頻率軸(只畫到Fs/2即可，由于y為實(shí)數(shù)，后面一半是對稱的% 畫出頻率幅度圖形，可以看出50Hz幅值最大%統(tǒng)計(jì)估計(jì)出來的頻率和真實(shí)頻率之間的誤差。

10、fmax=2*abs(Y(1;mn(j=0;for(i=1:1:NFFT/2+1if(2*abs(Y(i>fmaxfmax=2*abs(Y(i;mn(j=i;endendendstem(mn3、 實(shí)驗(yàn)結(jié)果圖1 未受干擾時(shí)信號(hào)的波形2 疊加了高斯白噪聲的信號(hào)波形3 估計(jì)出含高斯白噪聲的頻率從實(shí)驗(yàn)結(jié)果得知，估計(jì)出來的頻率是50HZ，跟信號(hào)的頻率相符。4 統(tǒng)計(jì)估計(jì)出來的頻率和真實(shí)頻率之間的誤差（代碼）重復(fù)進(jìn)行十次實(shí)驗(yàn)（每次加上的高斯白噪聲都是隨機(jī)的），每次得到的結(jié)果都是一樣的，估計(jì)的頻率都是52HZ，與原信號(hào)的頻率誤差很小，這是因?yàn)樗捎盟惴己茫宜x頻率比較低。5 分別改變信噪比、數(shù)據(jù)長度、頻率進(jìn)行試驗(yàn)，該頻率估計(jì)算法都可以比較精確地估計(jì)出含高斯白噪聲信號(hào)的頻率，所以這個(gè)算法的性能是優(yōu)越的。三、 心得體會(huì)：這次的實(shí)驗(yàn)內(nèi)容是FFT算法的應(yīng)用，當(dāng)然最核心的內(nèi)容是頻率的估計(jì)，前兩題都是比較常規(guī)的題目，在完成的過程中也沒有什么比較大的問題。第三題要求我們自己上網(wǎng)查詢資料，查找算法，然后進(jìn)行頻率估計(jì)。這個(gè)確