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

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

2、潤厲釤瘞睞櫪廡賴賃軔。矚慫潤厲釤瘞睞櫪廡賴賃。2. 譜分析參數(shù)選擇1) 設(shè)信號x(t)最高頻率為fc,對其進(jìn)行取樣得x(n),根據(jù)取樣定理，取樣頻率fs必須滿足:fs=2fc。聞創(chuàng)溝燴鐺險(xiǎn)愛氌譴凈禍測。聞創(chuàng)溝燴鐺險(xiǎn)愛氌譴凈禍。2) 設(shè)譜分辨率為F,則最小記錄時(shí)間tpmin=1/F ;取樣點(diǎn)數(shù)N 2fc/F為使用快速傅里葉變換(FFT)進(jìn)行譜分析，N還須滿足：N=2E (E為整數(shù))。3. 用FFT計(jì)算信號x(n)的頻譜。設(shè)x(n)為實(shí)信號快速傅里葉變換(FFT)是DFT的一種快速算法，其使得 DFT的運(yùn)算速度大為加快。1)對信號 x(n) 作 N 點(diǎn) FFT,得頻譜 X(k)(k=0N-1)X(

3、k)=XR(k)+jXI(k) (k=0N/2-1), XR(k) X(k)的實(shí)部；Xl(k) X(k)的虛部。殘騖樓諍錈瀨濟(jì)溆塹籟婭驟。殘騖樓諍錈瀨濟(jì)溆塹籟婭。Matlab 語句：Y=fft(x,N)其中：x-x(n);Y-X(k):222)幅頻譜:|X(k)|= xR(K) xI (K)，由于x(n)為實(shí)信號，因此|X(k)|對稱，Matlab 語句：abs(Y)iii)功率譜:PSD(k)=|X(k)|2/N=X(k)X*(k)/NMatlab 語句：PSD=Y .*conj(Y)/N其中：conj(Y)- X*(k)X(k) 的共軛4讀頻譜圖頻譜圖中任意頻率點(diǎn)k對應(yīng)實(shí)際頻率為:fk=f

4、s/N*k 。5用FFT實(shí)現(xiàn)線性卷積運(yùn)算用FFT實(shí)現(xiàn)y(n)=x(n)*h(n)的步驟為：1) 設(shè)x(n)及h(n)的長度分別為 N1和N2。為使循環(huán)卷積等于線性卷積，用補(bǔ)0的方法使x(n),h(n)長度均為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.mclc%輸入序列x(n)=1 1 1 1%x

5、(n)的長度Mclearxn=in put(x (n)=);M=le ngth(x n);xn=xn zeros(1,N-M); n=0:N-1;k=0:N-1;nk= n *k;wn=exp(-j*2*pi/N);wn K=w n.nk;xk=x n*wnK%補(bǔ)0,使xn長度為N%旋轉(zhuǎn)因子wn%作 x(n)的 DFT=xksubplot(211);stem(k,abs(xk),.);grid on;%顯示xk的幅頻譜(離散曲線)subplot(212);plot(k,abs(xk);grid on;%顯示xk的幅頻譜(連續(xù)曲線)運(yùn)行結(jié)果：問：由此得出怎樣的結(jié)論?答：n越大越接近原來的 dft

6、2. 理解 DIT-FFT 算法原理程序，并用它計(jì)算 X(k)=FFTR4(n), 分別取 N=4,8,16和 64, 繪出幅頻譜 |X(k)| 。%程序 DIT.mclearclcx=input(x= );%輸入序列N=input(N= );%做 fft 的點(diǎn)數(shù)x(length(x)+1:N)=zeros(1,N-length(x); 錐顧。%補(bǔ) 0 x(1:N) 釅錒極額閉鎮(zhèn)檜豬訣錐顧葒。 釅錒極額閉鎮(zhè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)

7、%M=2;while(M=N)W=exp(-2*j*pi/M);%旋轉(zhuǎn)因子 WV=1;for k=0:1:M/2-1%k 為每級蝶形運(yùn)算旋轉(zhuǎn)因子的個(gè)數(shù)for i=0:M:N-1%i 為各群的首序號p=k+i;q=p+M/2;A=x1(p+1);B=x1(q+1)*V;x1(p+1)=A+B;%本級蝶形運(yùn)算，x1最終存放X(k)x1(q+1)=A-B;endV=V*W;%旋轉(zhuǎn)因子W的變化endM=2*M;%第M級end%subplot(211);stem(x,.);grid on;%畫圖title(x( n);%標(biāo)題subplot(212);stem(abs(x1),.);grid on;%畫圖

8、title(|X(k)|);%標(biāo)題x(n)4-2111I|X(k)|3201.52.53.5x(n)|X(k)|x(n)|X(k)|0.510203040506070Hx(n)|X(k)|3.FFT譜分析設(shè)信號為 x(t)=sin(2n f1t)+sin(2鹵廡詒爾。彈貿(mào)攝爾霽斃攬磚鹵廡詒。%隨2機(jī)+噪聲，f1=50Hz, f2=120Hz,以取樣頻率 彈貿(mào)攝爾霽斃攬磚fs=1kHz對x(t)進(jìn)行取樣，樣本長度tp=0.25s，得x(n),對x(n)作256點(diǎn)FFT,得頻譜X(k),畫原信 號x(n),幅頻譜|X(k)|以及功率譜 PSD(k),對信號進(jìn)行譜分析。謀蕎摶篋飆鐸懟類蔣薔點(diǎn)鉍。謀蕎

9、摶篋飆鐸懟類蔣薔點(diǎn)。%程序 pufenxi.mclearclcfs=1000;t=0:1/fs:0.25;%時(shí)間范圍N=256;%做fft的點(diǎn)數(shù)f1=50;f2=120;%信號頻率s=si n( 2*pi*f1*t)+si n(2*pi*f2*t);%產(chǎn)生x(n)x=s+randn( size(t);時(shí)盡繼價(jià)騷巹。Y=fft(x,N);PSD=Y.*conj(Y)/N;f=fs/N*(0:N/2-1);齊。煢楨廣鰳鯡選塊網(wǎng)羈淚鍍。subplot(311);plot(x);subplot(312);plot(f,abs(Y(1:N/2);盡損鵪慘歷蘢鴛賴縈。subplot(313);plot(f

10、,PSD(1:N/2);媽羥為贍債蟶練淨(jìng)櫧。% 信號 +噪聲x（n）廈礴懇蹣駢時(shí)盡繼價(jià)騷巹癩。廈礴懇蹣駢%對x做N點(diǎn)fft%做功率譜%將頻率點(diǎn)轉(zhuǎn)化為實(shí)際頻率煢楨廣鰳鯡選塊網(wǎng)羈淚鍍%畫原信號%畫幅度譜（N/2點(diǎn)）鵝婭盡損鵪慘歷蘢鴛賴縈詰。鵝婭%畫功率譜（N/2點(diǎn)）籟叢媽羥為贍債蟶練淨(jìng)櫧撻。籟叢50-5畫出圖形窗口顯示的圖形，并注名每個(gè)圖形的含義。05010015020025030011.IE1I1i：1J 1la .ri_Ac/ - 11rrrrrrr1150100500050100150200250300350400450500604020005010015020025030035040045

11、05002）回答下列問題：i）觀察幅頻譜圖，可以發(fā)現(xiàn)，信號 x（n）含有的兩個(gè)頻率分量分別是 50.8Hz 和 121.1Hz。“plot(f,abs(Y(1:N/2); 改為 “p”lot(k,abs(Y(1:N/2);”重新運(yùn)行該程序并觀察幅頻譜圖，圖中兩峰值對應(yīng)的下標(biāo)分別是 13和 31。它們的含義為主頻點(diǎn)。再將該程序中的 N 改為 512, 重新運(yùn)行該程序并觀察幅頻譜圖，這時(shí)圖中 兩峰值對應(yīng)的下標(biāo)分別是 26和 61 。結(jié)果是否和上面的相同？ 不同 為什么？ N 不同。iii) 本例的頻譜分辨率 F是3.9Hz,改變f2=60Hz,問：在幅頻譜中，能否分辨fl和f2對應(yīng)的頻率分量？不能

12、。為什么？間隔小于頻譜分辨率。再改變 f2=52Hz, 問：在幅頻譜中，能否分辨 f1 和 f2 對應(yīng)的頻率分量？不能。 為什么？間隔小于頻譜分辨率。再改變 f2=600Hz ，在幅頻譜中， f2 對應(yīng)的頻率分量出現(xiàn)在398.45Hz；問：在 fs=1000Hz 的情況下，能否正確檢測 f2 對應(yīng)的頻率分量？不能。 為什么？不符合采樣定理。為了正確檢測f2對應(yīng)的頻率分量，則 fs至少取多少Hz?1200Hz。在該程序中改變 fs，驗(yàn)證 你的結(jié)論。 預(yù)頌圣鉉儐歲齦訝驊糴買闥。預(yù)頌圣鉉儐歲齦訝驊糴買。iv) 比較幅頻譜和功率譜，可以發(fā)現(xiàn)功率譜具有突出主頻點(diǎn)的特性。4. FFT 實(shí)現(xiàn)任意兩個(gè)序列的快

13、速卷積。%程序 fftjuanji.mclearclc x2=x2,zeros(1,N-N2);x1=input(x1=);x2=input(x2=);N1=length(x1);N2=length(x2);E=ceil(log2(N1+N2-1);N=2AE;%輸入序列%序列 x1(n),x2(n) 的長度%ceil-向+2方向取整%做 FFT 的點(diǎn)數(shù)X1=fft(x1,N);%對 x1 做 N 點(diǎn)的 fftX2=fft(x2,N);Y=X1.*X2; %數(shù)列 X1 和 X2 的乘積y=ifft(Y ,N)%對 y 做 N 點(diǎn)的 ifft結(jié)果分析：1)回到 MATLAB 窗口，鍵入：x1=1

14、 1 1, x2=1 2, 回車。結(jié)果 :y= 1 3 3 22)問：可用 Matlab 中的什么函數(shù)驗(yàn)算上述卷積結(jié)果？ Y=conv(x1,x2)5. 利用譜分析觀察太陽黑子周期性。以100年中記錄到的太陽黑子出現(xiàn)次數(shù)為信號x(n),對x(n)作功率譜，從中觀察太陽黑子周期性。%程序 taiyangheizi.mclearclcx=101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 . 滲釤嗆儼勻諤鱉調(diào)硯錦鋇絨。滲釤嗆儼勻諤鱉 調(diào)硯錦鋇。132 131 118 90 67 60 47 41 21 16 6 4 7 14 34 45 43

15、48 . 鐃誅臥瀉噦圣騁貺頂廡縫勵(lì)。鐃誅臥瀉噦圣騁 貺頂廡縫。42 28 10 8 2 0 1 5 12 14 35 46 41 30 24 16 7 4 2 8 . 擁締鳳襪備訊顎輪爛薔報(bào)贏。 擁締鳳襪備訊顎輪爛薔 報(bào)。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 . 壇摶鄉(xiāng)囂懺蔞鍥鈴氈淚躋馱。 壇摶鄉(xiāng)囂懺蔞鍥鈴 氈淚躋。77 59 44 47 30 16 7 37 74;%100 年中太陽黑子出現(xiàn)的次數(shù)subplot(211);plot(x)N=128;fs=1;s=x-mea n( x);Y=fft(s,n);PSD=Y.*conj(Y)/N;f=fs/N*(0:N/2-1);% 畫 x(n)%fs=1Hz,N=128 點(diǎn)%對x作零均值化處理(去除直流分量)%對s做N點(diǎn)fft%做功率譜PSD%將頻率定標(biāo)為實(shí)際頻率fsubplot(212);plot(f,PSD(1:N/2);%畫功率譜(N/2 點(diǎn))4x 10填寫空格中的畫圖語句并繪出結(jié)果圖形。1tblL/i f