實(shí)驗(yàn)三DFT和FFT頻譜分析

1、實(shí)驗(yàn)三 DFT 和 FFT 頻譜分析一、實(shí)驗(yàn)?zāi)康?.掌握 DFT 頻譜分析的原理與編程方法。2.理解 FFT 算法的編程思想。2.熟練掌握利用FFT 對(duì)信號(hào)作頻譜分析，包括正確地進(jìn)行參數(shù)選擇、畫頻譜及讀頻譜圖。3. 利用 FFT 頻譜分析進(jìn)行快速卷積和太陽(yáng)黑子周期性檢測(cè)。二、實(shí)驗(yàn)環(huán)境1.Windows xp 以上操作系統(tǒng)2.安裝 MATLAB2007a軟件三、實(shí)驗(yàn)原理1. 離散傅里葉變換 (DFT)設(shè)序列為x(n),長(zhǎng)度為 N, 則N 1X(e j k)=DFTx(n)=x(n) e-j k n,n 0其中 k=2 -,M1),通常 M>N, 以便觀察頻譜的細(xì)節(jié)。|X(ej kk (k=

2、0,1,2,)|-x(n) 的幅頻譜。M2.譜分析參數(shù)選擇1) 設(shè)信號(hào) x(t) 最高頻率為fc,對(duì)其進(jìn)行取樣得x(n), 根據(jù)取樣定理，取樣頻率fs必須滿足：fs>=2fc。2)設(shè)譜分辨率為 F,則最小記錄時(shí)間 tpmin=1/F；取樣點(diǎn)數(shù)N2fc/F;為使用快速傅里葉變換(FFT) 進(jìn)行譜分析， N 還須滿足： N=2E（ E為整數(shù)）。3.用 FFT 計(jì)算信號(hào) x(n)的頻譜。 設(shè) x(n)為實(shí)信號(hào) 快速傅里葉變換（FFT ）是 DFT 的一種快速算法，其使得DFT的運(yùn)算速度大為加快。1)對(duì)信號(hào) x(n)作 N點(diǎn) FFT, 得頻譜 X(k)(k=0N-1)X(k)=X R (k)+j

3、X I (k) (k=0N/2-1)， X R(k) X(k) 的實(shí)部； XI (k) X(k) 的虛部。Matlab 語(yǔ)句： Y=fft(x,N)其中： x-x(n);Y-X(k)222) 幅頻譜 :|X(k)|=xR (K )xI ( K )，由于 x(n) 為實(shí)信號(hào)，因此 |X(k)|對(duì)稱，Matlab 語(yǔ)句： abs(Y)iii) 功率譜 :PSD(k)=|X(k)| 2/N=X(k)X * (k)/NMatlab 語(yǔ)句： PSD=Y.*conj(Y)/N其中： conj(Y)- X * (k)X(k) 的共軛 4.讀頻譜圖頻譜圖中任意頻率點(diǎn) k對(duì)應(yīng)實(shí)際頻率為 :f k =fs/N*k

4、。5.用 FFT 實(shí)現(xiàn)線性卷積運(yùn)算用FFT 實(shí)現(xiàn) y(n)=x(n)*h(n) 的步驟為：1) 設(shè)x(n)及 h(n) 的長(zhǎng)度分別為 N 1和 N2。為使循環(huán)卷積等于線性卷積，用補(bǔ)0的方法使 x(n),h(n) 長(zhǎng)度均為 N，則 N須滿足 NN1+N2-1；為用 FFT 計(jì)算 DFT,則 N還須滿足 N=2E (E為整數(shù) )。2) 用 FFT 計(jì)算 X(k),H(k) (N點(diǎn) )。3)Y(k)= ifft; y(n)=ifftY(K)。四、實(shí)驗(yàn)內(nèi)容1.根據(jù)公式設(shè)計(jì)DFT 原理程序，并計(jì)算：x(n)=1,1,1,1 的 4,16,64點(diǎn) DFT 并繪圖。%DFT/IDFT程序DFT.mclccl

5、earxn=input('x(n)= ');M=length(xn);% 輸入序列 x(n)=1 1 1 1%x(n) 的長(zhǎng)度 MN=input(' 變換區(qū)間 N=');% 變換區(qū)間 Nxn=xn zeros(1,N-M);% 補(bǔ) 0,使 xn長(zhǎng)度為 Nn=0:N-1;k=0:N-1;nk=n'*k;wn=exp(-j*2*pi/N);% 旋轉(zhuǎn)因子 wnwnK=wn.nk;xk=xn*wnK% 作 x(n)的 DFT=xksubplot(211);stem(k,abs(xk),'.');grid on;% 顯示 xk的幅頻譜 (離散曲線

6、) subplot(212);plot(k,abs(xk);grid on;% 顯示 xk的幅頻譜 (連續(xù)曲線 )運(yùn)行結(jié)果：4321000.511.522.534321000.511.522.5343210051015432100510154321001020304050607043210010203040506070問(wèn)：由此得出怎樣的結(jié)論？答： n越大越接近原來(lái)的dft2.理解 DIT-FFT 算法原理程序，并用它計(jì)算X(k)=FFTR 4(n), 分別取 N=4,8,16和 64,繪出幅頻譜|X(k)| 。% 程序DIT.mclearclcx=input('x= ');%輸

7、入序列N=input('N= ');%做 fft的點(diǎn)數(shù)x(length(x)+1:N)=zeros(1,N-length(x);% 補(bǔ) 0 x(1:N)l=log2(N);x1=zeros(1,N);for j1=1:N% 倒序x1(j1)=x(bin2dec(fliplr(dec2bin(j1-1,l)+1);end%FFT(DIT)%M=2;while(M<=N)W=exp(-2*j*pi/M);% 旋轉(zhuǎn)因子 WV=1;for k=0:1:M/2-1%k 為每級(jí)蝶形運(yùn)算旋轉(zhuǎn)因子的個(gè)數(shù)for i=0:M:N-1%i 為各群的首序號(hào)p=k+i;q=p+M/2;A=x1(p

8、+1);B=x1(q+1)*V;x1(p+1)=A+B;% 本級(jí)蝶形運(yùn)算，x1最終存放 X(k)x1(q+1)=A-B;endV=V*W;% 旋轉(zhuǎn)因子 W 的變化endM=2*M;%第M級(jí)end%subplot(211);stem(x,'.');grid on;%畫圖title('x(n)');%標(biāo)題subplot(212);stem(abs(x1),'.');grid on;%畫圖title('|X(k)|');%標(biāo)題x(n)10.5011.522.533.54|X(k)|4321011.522.533.54x(n)10.501

9、2345678|X(k)|4321012345678x(n)10.500246810121416|X(k)|432100246810121416x(n)10.50010203040506070|X(k)|432100102030405060703.FFT 譜分析設(shè)信號(hào)為 x(t)=sin(2 ff以取樣頻率1t)+sin(22t)+ 隨機(jī)噪聲， f1=50Hz, f 2=120Hz,fs=1kHz 對(duì) x(t) 進(jìn)行取樣 ,樣本長(zhǎng)度 tp=0.25s,得 x(n), 對(duì) x(n)作 256點(diǎn) FFT, 得頻譜 X(k), 畫原信號(hào)x(n),幅頻譜 |X(k)| 以及功率譜 PSD(k), 對(duì)信

10、號(hào)進(jìn)行譜分析。% 程序 pufenxi.mclearclcfs=1000;t=0:1/fs:0.25;%時(shí)間范圍N=256;%做 fft的點(diǎn)數(shù)f1=50;f2=120;%信號(hào)頻率s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%產(chǎn)生 x(n)x=s+randn(size(t);% 信號(hào) +噪聲 x(n)Y=fft(x,N);%對(duì) x做 N點(diǎn) fftPSD=Y.*conj(Y)/N;%做功率譜f=fs/N*(0:N/2-1);%將頻率點(diǎn)轉(zhuǎn)化為實(shí)際頻率subplot(311);plot(x);%畫原信號(hào)subplot(312);plot(f,abs(Y(1:N/2);%畫幅度譜（

11、N/2點(diǎn)）subplot(313);plot(f,PSD(1:N/2);%畫功率譜（ N/2 點(diǎn)）1) 畫出圖形窗口顯示的圖形，并注名每個(gè)圖形的含義。50-550100150200250300015010050050100150200250300350400450500060402005010015020025030035040045050002)回答下列問(wèn)題：i)觀察幅頻譜圖，可以發(fā)現(xiàn)，信號(hào)x(n)含有的兩個(gè)頻率分量分別是50.8Hz 和121.1Hz 。ii)在該程序中的“ f=fs/N*(0:N/2-1); ”下添加 “ k=0:N/2-1; ”，“ plot(f,abs(Y(1:N/2

12、); 改為 “”plot(k,abs(Y(1:N/2);重新運(yùn)行該程序并觀察幅頻譜圖，圖中兩峰值對(duì)應(yīng)的下標(biāo)分別是”13 和再將該程序中的31。它們的含義為主頻點(diǎn)。N 改為 512, 重新運(yùn)行該程序并觀察幅頻譜圖，這時(shí)圖中兩峰值對(duì)應(yīng)的下標(biāo)分別是26和61。結(jié)果是否和上面的相同？不同為什么？N不同。iii) 本例的頻譜分辨率F是3.9Hz,改變 f2=60Hz, 問(wèn)：在幅頻譜中，能否分辨f1和 f2對(duì)應(yīng)的頻率分量？不 能。為什么？間隔小于頻譜分辨率。再改變 f2 =52Hz, 問(wèn)：在幅頻譜中，能否分辨f1和f2對(duì)應(yīng)的頻率分量？不能。為什么？間隔小于頻譜分辨率。再改變 f2 =600Hz ，在幅頻譜

13、中， f2對(duì)應(yīng)的頻率分量出現(xiàn)在398.45Hz ；問(wèn)：在 fs=1000Hz 的情況下，能否正確檢測(cè)f2對(duì)應(yīng)的頻率分量？不能。為什么？不符合采樣定理。為了正確檢測(cè) f2對(duì)應(yīng)的頻率分量， 則fs至少取多少 Hz?1200Hz 。在該程序中改變 fs，驗(yàn)證你的結(jié)論。iv)比較幅頻譜和功率譜，可以發(fā)現(xiàn)功率譜具有突出主頻點(diǎn)的特性。4.FFT 實(shí)現(xiàn)任意兩個(gè)序列的快速卷積。% 程序 fftjuanji.mclearclcx1=input('x1=');x2=input('x2=');%輸入序列N1=length(x1);N2=length(x2);% 序列 x1(n),x2

14、(n) 的長(zhǎng)度E=ceil(log2(N1+N2-1);%ceil- 向 +方向取整N=2E;%做 FFT 的點(diǎn)數(shù)x1=x1,zeros(1,N-N1);%補(bǔ)零使 x1,x2長(zhǎng)度為 Nx2=x2,zeros(1,N-N2);X1=fft(x1,N);%對(duì) x1做 N點(diǎn)的 fftX2=fft(x2,N);Y=X1.*X2;%數(shù)列 X1 和 X2 的乘積y=ifft(Y,N)%對(duì) y做 N點(diǎn)的 ifft結(jié)果分析：1)回到 MATLAB 窗口，鍵入：x1=1 1 1, x2=1 2, 回車。結(jié)果 :y=13322)問(wèn)：可用 Matlab中的什么函數(shù)驗(yàn)算上述卷積結(jié)果？Y=conv(x1,x2)5. 利

15、用譜分析觀察太陽(yáng)黑子周期性。以 100年中記錄到的太陽(yáng)黑子出現(xiàn)次數(shù)為信號(hào) x(n),對(duì) x(n)作功率譜，從中觀察太陽(yáng)黑子周期性。% 程序 taiyangheizi.m clearclcx=101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 .132 131 118 90 67 60 47 41 21 16 6 4 7 14 34 45 43 48 .17 36 50 62 67 71 48 28 8 13 57 122 138 103 86 63 37 24 .11 15 40 62 98 124 96 66 64 54 39 21 7 4 23 55 94 96 .77 59 44 47 30 16 7 37 74;%100 年中太陽(yáng)黑子出現(xiàn)的次數(shù)subplot(211);plot(x)% 畫 x(n)N=128;fs=1;%fs=1Hz,N=128 點(diǎn)s=x-mean(x);% 對(duì) x作零均值化處理（去除直流分量）Y= fft(s,n);% 對(duì) s做 N 點(diǎn)fftPSD=Y.*conj(Y)/N;% 做功率譜 PSDf=fs/N*(0:N/2-1);% 將頻率定標(biāo)為實(shí)際頻率 fsubplot(212);plot(f,PSD(1:N/