實驗二用FFT作譜分析

1、實驗二 用FFT作譜分析 1. 實驗?zāi)康?1) 進一步加深DFT算法原理和基本性質(zhì)的理解(因為FFT只是DFT的一種快速算法，所以FFT的運算結(jié)果必然滿足DFT的基本性質(zhì))。(2) 熟悉FFT算法原理和FFT子程序的應(yīng)用。(3) 學(xué)習(xí)用FFT對連續(xù)信號和時域離散信號進行譜分析的方法，了解可能出現(xiàn)的分析誤差及其原因，以便在實際中正確應(yīng)用FFT。2. 實驗步驟(1) 復(fù)習(xí)DFT的定義、性質(zhì)和用DFT作譜分析的有關(guān)內(nèi)容。(2) 復(fù)習(xí)FFT算法原理與編程思想，并對照DIT-FFT運算流圖和程序框圖， 讀懂本實驗提供的FFT子程序。(3) 編制信號產(chǎn)生子程序， 產(chǎn)生以下典型信號供譜分析用： (4) 編寫

2、主程序。/* DIT-FFT函數(shù)（C語言）*/fft基2DITFFT函數(shù)要求：指向復(fù)數(shù)數(shù)組指針X，F(xiàn)FT長度為2m，m為正整數(shù)FFT輸出結(jié)果放在輸入復(fù)數(shù)數(shù)組中。/* 計算N點FFT子程序*/* xr：=信號序列實部，xi：=信號序列虛部，N：=FFT變換區(qū)間長度N=2M */* 如果信號長度小于N，應(yīng)該給xr，xi后面補0 */* 計算如果X（K）的實部和虛部分別儲存在數(shù)組xr和xi中 */Void Fft(double xr, double xi, int N, int M)int L, B, J, P, k, i;double rPartKB, iPartKB;double rCf128,

3、 iCf128/* 計算旋轉(zhuǎn)因子 */double PI2= 8.0 * atan(1.0);for(i=0; iN; i+)rCfi=cos(i*PI2/N);iCfi=sin(i*PI2/N);ChangeOrder( xr, xi, N );/* 計算各級蝶形 */for(L=1; L=M; L+)B=(int)(pow(2, (L-1)+0.5);for(J=0; J=B-1; J+)P=J*(int)(pow(2,(M-L)+0.5);for(k=J; k=N-1; k+=(int)(pow(2,L)+0.5)rPartKB = xrk+B*rCfP-xik+B*iCfP;iPart

4、KB = xik+B*rCfP+xrk+B*iCfPxrk+B = xrk- rPartKB;xik+B = xik- iPartKB;xrk = xrk + rPartKB;xik = xik + iPartKB;/* 倒序子程序*/void ChangeOrdor(double xr, double xi, int N )int LH, N1, I, J, K;double T;LH = N/2; J = LH; N1 = N 2;for(I=1;I=N1;I+)if(I=K)J = J-K;K = (int)(K/2+0.5);J=J+K;(1) 按實驗內(nèi)容要求，上機實驗并寫出實驗報告。

5、本實驗采用的是MATLAB語言，因此FFT子程序直接調(diào)用MATLAB語言中的FFT函數(shù)就可以實現(xiàn)。下面給出完整的MATLAB程序%實驗二，用FFT做譜分析b=menu(請選擇信號x1(n)-x8(n),x1(n),x2(n),x3(n),x4(n),x5(n),x6(n),x7=x4+x5,x8=x4+jx5,Exit);if b=9 b=0;endi=0;close all;while(b)if b=6temp=menu(請選擇FFT變換區(qū)間長度N,N=16,N=32,N=64);if temp=1 N=16; elseif temp=2 N=32; else N=64; endfs=64;

6、n=0:N-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);elsetemp=menu(請選擇FFT變換區(qū)間長度N,N=8,N=16,N=32); if temp=1 N=8; elseif temp=2 N=16; else N=32; endif b=1x=1 1 1 1 0 0 0 0;else if b=2 x=1 2 3 4 4 3 2 1; else if b=3x=4 3 2 1 1 2 3 4;else if b=4 n=0:N-1; x=cos(0.25*pi*n); else if b=5 n=0:N-1;x=sin(p

7、i*n)/8);else if b=7 n=0:N-1; x=cos(n*pi/4)+sin(n*pi/8); else if b=8n=0:N-1;x=cos(n*pi/4)+j*sin(n*pi/8);end endend end end end end end%TO Calculate FFTf=fft(x,N);i=i+1;figure(i);printf(x,abs(f),abs(N),abs(b);if N=16 if b=7 k=conj(f); x4=(f+k)/2;%ReX7(k)=x4(k) figure(i+2); subplot(2,2,1); stem(abs(x4)

8、,.); xlabel(k); ylabel(|X4(k)|); title(恢復(fù)后的X4(k); x5=(f-k)/2;%jImX7(k)=X5(k) subplot(2,2,3); Stem(abs(x5),.); xlabel(k); ylabel(|X5(k)|); title(恢復(fù)后的X5(k); end if b=8 k(1)=conj(f(1); for m=2:N k(m)=conj(f(N-m+2); end fe=(x+k)/2;%求X8(k)的共軛對稱分量 fo=(x-k)/2;%求X8(k)的共軛反對稱分量 xr=ifft(fe,N);%xr=x4(n) b=4; fi

9、gure(i+1) printf(xr,abs(fe),abs(N),abs(b); xi=ifft(fo,N)/j;%xi=x5(n) b=5; figure(i+2) printf(xi,abs(f),abs(N),abs(b); endendb=menu(請選擇信號x1(n)-x8(n),x1(n),x2(n),x3(n),x4(n),x5(n),x6(n),x7=x4+x5,x8=x4+jx5,Exit);if b=9 b=0;endclose all;end(5) 按實驗內(nèi)容要求，上機實驗，并寫出實驗報告。 3. 上機實驗內(nèi)容(1) 對 2 中所給出的信號逐個進行譜分析。下面給出針對

10、各信號的FFT變換區(qū)間N以及對連續(xù)信號x6(t)的采樣頻率Fs，供實驗時參考。 x1(n)，x2(n)，x3(n)，x4(n)，x5(n)：N=8，16x6(t)：Fs=64Hz，N=16，32，64x1(n)N=8 N=16 由以上兩個圖分析如下： 離散傅里葉變換的N點變換在頻域范圍內(nèi)表現(xiàn)為對傅里葉變換即Z變換在單位圓上的抽樣。所以N取8點時，k=0,1,2,3,4,5,6,7與N取16點時，k=0,2,4,6,8,10,12,14的離散傅里葉變換值對應(yīng)相等，即它們都等于原信號在w=0、/8、/4、3/8、4/8、5/8、6/8、6/8的傅里葉變換，這在上面兩圖可以明顯看出。所以，離散傅里葉

11、變換實際上是對該序列在頻域范圍內(nèi)以2/N的間隔進行抽樣。x2(n) N=8 N=16x3(n) N=8 N=16X4(n) N=8 N=16從以上兩組圖可以看出，若按X1(n)進行分析，則明顯不對。現(xiàn)分析如下：原信號周期為16，所以當(dāng)N=8時，未能取完一個周期的值，N=16則取完了一個周期的值，所以這是兩個不同的序列，所以按照X1(n)的分析方式是不對的，因為本身它們的傅里葉變換就是不一樣的。由于離散傅里葉變換是該序列周期延拓后所對應(yīng)的傅里葉級數(shù)變換的主值序列，所以，當(dāng)N=16時，所得的DFT值與X5(n)的傅里葉級數(shù)變換的主值序列是一致的，而N=8時是X5(n)的部分序列的周期延拓后的傅里葉

12、級數(shù)變換的主值序列，因此兩者的值是不同的。X6(n) N=16 N=32N=64原連續(xù)信號的周期為0.5，當(dāng)采樣頻率Fs=64Hz時，所形成的序列周期為0.5*64=32。所以只有N32,才能取完一個周期的序列。這一點，從上面三個圖可以清晰看出。其中N=32和N=16的圖形分析，可以參考x4(n)的分析。N=32和N=64的圖形分析，可以參考x5(n)的分析。(2) 令x(n)=x4(n)+x5(n)，用FFT計算8點和16點離散傅里葉變換。 X(k)=DFT x(n) X7 N=8 N=16兩者都不能顯示一個周期內(nèi)的所有序列，盡管N=8時的序列是N=16時序列的一部分，但是它們確屬兩個不同的

13、序列。所以它們的傅里葉變換不同，即不能按照x1(n)的進行分析而且周期延拓后所取得傅里葉級數(shù)的主值序列不同，即DFT變換值不同。(3) 令x(n)=x4(n)+jx5(n)， 重復(fù)(2)。X8 N=8N=16通過結(jié)合x4(n)和x5(n)的頻譜分析，從以上的圖可以看出，將原信號的DFT變換分為共軛對稱部分合共個反對稱部分，則可以得出原信號的實部對應(yīng)離散傅里葉變換的共軛對稱部分，原信號的虛部對應(yīng)離散信號的共軛反對稱部分。4. 思考題(1) 在N=8時，x2(n)和x3(n)的幅頻特性會相同嗎？為什么？N=16呢？時一樣，時不一樣。因為DFT變換可以看成是將該序列進行周期延拓后的傅里葉級數(shù)變換的主值序列。當(dāng)N=8時，兩序列進行周期延拓后序列相同，所以其傅里葉級數(shù)變換的主值序列也相同，進而DFT變換也相同。而當(dāng)N=16時，兩序列進行周期延拓后序列不相同，所