實(shí)驗(yàn)二用FFT作譜分析

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

2、寫(xiě)主程序。/* DIT-FFT函數(shù)（C語(yǔ)言）*/fft基2DITFFT函數(shù)要求：指向復(fù)數(shù)數(shù)組指針X，F(xiàn)FT長(zhǎng)度為2m，m為正整數(shù)FFT輸出結(jié)果放在輸入復(fù)數(shù)數(shù)組中。/* 計(jì)算N點(diǎn)FFT子程序*/* xr：=信號(hào)序列實(shí)部，xi：=信號(hào)序列虛部，N：=FFT變換區(qū)間長(zhǎng)度N=2M */* 如果信號(hào)長(zhǎng)度小于N，應(yīng)該給xr，xi后面補(bǔ)0 */* 計(jì)算如果X（K）的實(shí)部和虛部分別儲(chǔ)存在數(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/* 計(jì)算旋轉(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 );/* 計(jì)算各級(jí)蝶形 */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;iPar

4、tKB = 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) 按實(shí)驗(yàn)內(nèi)容要求，上機(jī)實(shí)驗(yàn)并寫(xiě)出實(shí)驗(yàn)報(bào)告

5、。本實(shí)驗(yàn)采用的是MATLAB語(yǔ)言，因此FFT子程序直接調(diào)用MATLAB語(yǔ)言中的FFT函數(shù)就可以實(shí)現(xiàn)。下面給出完整的MATLAB程序%實(shí)驗(yàn)二，用FFT做譜分析b=menu(請(qǐng)選擇信號(hào)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(請(qǐng)選擇FFT變換區(qū)間長(zhǎng)度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(請(qǐng)選擇FFT變換區(qū)間長(zhǎng)度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(

7、pi*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)的共軛對(duì)稱分量 fo=(x-k)/2;%求X8(k)的共軛反對(duì)稱分量 xr=ifft(fe,N);%xr=x4(n) b=4; f

9、igure(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(請(qǐng)選擇信號(hào)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) 按實(shí)驗(yàn)內(nèi)容要求，上機(jī)實(shí)驗(yàn)，并寫(xiě)出實(shí)驗(yàn)報(bào)告。 3. 上機(jī)實(shí)驗(yàn)內(nèi)容(1) 對(duì) 2 中所給出的信號(hào)逐個(gè)進(jìn)行譜分析。下面給出針

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

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

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

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