離散傅里葉變換和快速傅里葉變換

1、實(shí)驗(yàn)報(bào)告課程名稱：信號(hào)分析與處理 指導(dǎo)老師 成績(jī)： 實(shí)驗(yàn)名稱：離散傅里葉變換和快速傅里葉變換 實(shí)驗(yàn)類型： 同組學(xué)生姓名： 一、實(shí)驗(yàn)?zāi)康暮鸵螅ū靥睿┒?、?shí)驗(yàn)內(nèi)容和原理（必填）三、主要儀器設(shè)備（必填）四、操作方法和實(shí)驗(yàn)步驟五、實(shí)驗(yàn)數(shù)據(jù)記錄和處理六、實(shí)驗(yàn)結(jié)果與分析（必填）七、討論、心得一、實(shí)驗(yàn)?zāi)康暮鸵?. 掌握DFT的原理和實(shí)現(xiàn)2. 掌握FFT的原理和實(shí)現(xiàn)，掌握用FFT對(duì)連續(xù)信號(hào)和離散信號(hào)進(jìn)行譜分析的方法。二、實(shí)驗(yàn)內(nèi)容和原理2.1 DTFT和DFT序列x(n)的離散事件傅里葉變換（DTFT）表示為：，如果x(n)為因果有限長(zhǎng)序列，n=0,1,.,N-1，則x(n)的DTFT表示為：，x(n)的離

2、散傅里葉變換（DFT）表達(dá)式為：，序列的N點(diǎn)DFT是DTFT在 0,2上的N點(diǎn)等間隔采樣，采樣間隔為2/N。通過(guò)DFT，可以完成由一組有限個(gè)信號(hào)采樣值x(n)直接計(jì)算得到一組有限個(gè)頻譜采樣值X(k)。X(k)的幅度譜為，XR(k)和XI(k)分別為X(k)的實(shí)部和虛部。X(k)的相位譜為。離散傅里葉反變換（IDFT）定義為。2.2 FFT快速傅里葉變換（FFT）是DFT的快速算法，它減少了DFT的運(yùn)算量，使數(shù)字信號(hào)的處理速度大大提高。三、主要儀器設(shè)備PC一臺(tái)，matlab軟件四、實(shí)驗(yàn)內(nèi)容4.1第一題求有限長(zhǎng)離散時(shí)間信號(hào)x(n)的離散時(shí)間傅里葉變換（DTFT）X(ej)并繪圖。（1）已知；（2）

3、已知。4.1.1理論分析1) 由DTFT計(jì)算式， X（）是實(shí)數(shù)，可以直接作出它的圖像。Figure 1 X（）曲線2) 由DTFT計(jì)算式：可以發(fā)現(xiàn)X（）周期為2；而X（）的相位在2周期內(nèi)有約十次振蕩。4.1.2編程計(jì)算作圖編寫(xiě)一個(gè)計(jì)算DTFT的函數(shù)。function DTFT(x,n1,n2)w=-2*pi:2*pi/1000:2*pi; %表示X=zeros(size(w); for i=n1:n2 %DTFT計(jì)算式 X=X+x(i-n1+1)*exp(-1)*j*w*i);endangle(X);subplot(2,1,1);plot(w, abs(X),r);xlabel(Omega);

4、ylabel(|X(Omega)|);hold on; %作幅頻圖subplot(2,1,2);plot(w,angle(X),b);xlabel(Omega);ylabel(angle(Omega); %作相頻圖end輸入序列x，和n的取值范圍，即可計(jì)算其DTFT。1) 輸入：x=1 1 1 1 1;DTFT(x,-2,2);（因?yàn)閄（）是實(shí)數(shù)，所以實(shí)際計(jì)算過(guò)程中對(duì)相頻曲線取了絕對(duì)值）結(jié)果：Figure 2 X（）的頻譜可以看出，X（）的相位只有0和兩種取值，X（）是實(shí)函數(shù)，而且其幅度頻譜與理論計(jì)算得到的相同。2) 輸入：n=0:10;x=2.n;DTFT(x,0,10)結(jié)果：Figure

5、3 第1題（2）中X（）的頻譜4.2第二題已知有限長(zhǎng)序列x(n)=8,7,9,5,1,7,9,5,試分別采用DFT和FFT求其離散傅里葉變換X(k)的幅度、相位圖。4.2.1理論分析由FFT蝶形運(yùn)算得到，X(k)= 51,7,-9-j4,7,3,7,-9+j4,74.2.2編程計(jì)算作圖1. DFT編寫(xiě)一個(gè)計(jì)算DTFT的函數(shù)。DFT（序列x，長(zhǎng)度N）function DFT(x,N)k=0:N-1;X=zeros(size(k);for n=0:N-1X=X+x(n+1)*exp(-1)*j*2*pi/N*n*k);%DFT計(jì)算式endsubplot(2,1,1);stem(k,abs(X),.

6、); xlabel(k);ylabel(|X(k)|);hold on; %幅頻圖subplot(2,1,2);stem(k,angle(X),*);xlabel(k);ylabel(Angle(k);%相頻圖end輸入：x=8 7 9 5 1 7 9 5;DFT(x,8);結(jié)果：Figure 4 第2題DFT結(jié)果2. FFT編寫(xiě)一個(gè)利用matlab自帶函數(shù)計(jì)算FFT并繪圖的函數(shù)FFT1（序列x，長(zhǎng)度N）function FFT1(x,N)X=fft(x,8); %用自帶的fft函數(shù)計(jì)算k=0:N-1;subplot(2,1,1);stem(k,abs(X),.); xlabel(k);yla

7、bel(|X(k)|);hold on; %幅頻圖subplot(2,1,2);stem(k,angle(X),*);xlabel(k);ylabel(Angle(k);%相頻圖end輸入：x=8 7 9 5 1 7 9 5;FFT1(x,8);結(jié)果：Figure 5 第2題FFT結(jié)果因?yàn)镕FT只是DFT的一種快速算法，所以FFT的結(jié)果與DFT結(jié)果相同。DFT和FFT的結(jié)果，符合理論計(jì)算得到的，X(k)= 51,7,-9-j4,7,3,7,-9+j4,7 。4.3第三題已知連續(xù)時(shí)間信號(hào)x(t)=3cos8t, X()= ，該信號(hào)從t=0開(kāi)始以采樣周期Ts=0.1 s進(jìn)行采樣得到序列x(n)，試

8、選擇合適的采樣點(diǎn)數(shù)，分別采用DFT和 FFT求其離散傅里葉變換X(k)的幅度、相位圖，并將結(jié)果與X(k)的幅度、相位圖，并將結(jié)果與X()相比較。4.3.1理論分析1. 原信號(hào)的頻譜：X()=，只在8不為0. 且在8處相位為0 。2. 采樣角頻率s=2028，滿足采樣定理。3. 采樣后的信號(hào)，為X()以20為周期的延拓。所以只在（8+20k） （k為任意整數(shù)）處不為0. 如取區(qū)間0,20內(nèi)，只有8和16處不為0 。進(jìn)行N點(diǎn)DFT后，將20的區(qū)間映射為0,N區(qū)間。理想情況下僅在 n= 和兩處不為0 。4. ，周期為5，所以取采樣點(diǎn)數(shù)為5的倍數(shù)時(shí)，不會(huì)發(fā)生泄漏；而采樣點(diǎn)數(shù)不是5的倍數(shù)時(shí)，則會(huì)發(fā)生泄漏

9、。Figure 6 原始信號(hào)的頻譜X（）4.3.2編程計(jì)算作圖編寫(xiě)一個(gè)獲得信號(hào)的N點(diǎn)樣本的函數(shù) sample（點(diǎn)數(shù)N）function x=sample(N)t=0:0.1:(N-1)*0.1; %0.1s為間隔x=3*cos(8*pi*t); %x即采樣結(jié)果。End輸入:X=sample(N);FFT1(X,N);即可獲得采樣N點(diǎn)的頻譜圖。因?yàn)镈FT結(jié)果與FFT是完全一樣的，所以這里只使用FFT作圖。取采樣點(diǎn)數(shù)N=5 16 20 104 獲得以下頻譜圖：Figure 7 N=5Figure 8 N=16Figure 9 N=20Figure 10 N=104可以看出，N=5和20時(shí)，由于是周

10、期的整數(shù)倍，頻譜只有兩條譜線，且滿足前面理論計(jì)算得出的公式n= 和，沒(méi)有發(fā)生泄漏，且這兩條譜線對(duì)應(yīng)的相位是0.所以頻譜與原信號(hào)頻譜在形式上時(shí)相同的。而N=16和N=54時(shí)，則都發(fā)生了頻譜泄漏，頻譜與原信號(hào)頻譜就很不同了。但相比之下N=54時(shí)譜線更加接近原譜線。驗(yàn)證了“為減小泄漏誤差，如果待分析的信號(hào)實(shí)現(xiàn)不知道確切周期，則截取較長(zhǎng)時(shí)間長(zhǎng)度的樣點(diǎn)進(jìn)行分析”這個(gè)說(shuō)法。同時(shí)也可以發(fā)現(xiàn)，雖然幅頻圖中顯示幅值為0，但相頻圖中相應(yīng)的位置仍有譜線。這可能是matlab浮點(diǎn)運(yùn)算造成的誤差，即本來(lái)為0處其實(shí)是一個(gè)非常小的復(fù)數(shù)，所以仍有一定相位。4.4第四題4.4.1理論分析若噪聲信號(hào)較小，則采樣后的頻譜仍能較準(zhǔn)確

11、地反映原信號(hào)的特征。4.4.2編程計(jì)算對(duì)原采樣程序稍加改編，加入一個(gè)噪聲信號(hào)p*randn(1,N)。p表示噪聲信號(hào)的強(qiáng)度。function x=samplenoise(N,p)t=0:0.1:(N-1)*0.1;x=3*cos(8*pi*t)+p*randn(1,N);end取采樣點(diǎn)數(shù)N=20進(jìn)行分析。輸入：X=samplenoise(20,p); %取P=1和10兩種情況。FFT1(X,20);Figure 11 N=20 噪聲較?。╬=1）Figure 12 N=20 噪聲較大（p=10）可見(jiàn)，較小的噪聲對(duì)信號(hào)的頻譜的影響不大，仍能較精確地獲得頻譜圖。而噪聲較大時(shí)則很難準(zhǔn)確獲得原信號(hào)的頻

12、譜。4.5第五題3.5已知序列，X(k)是x(n)的6點(diǎn)DFT，設(shè)。（1） 若有限長(zhǎng)序列y(n)的6點(diǎn)DFT是，求y(n)。（2） 若有限長(zhǎng)序列w(n)的6點(diǎn)DFT W(k)是的實(shí)部，求w(n)。（3） 若有限長(zhǎng)序列q(n)的3點(diǎn)DFT是，k=0,1,2，求q(n)。由題意得到：x(n)=4,3, 2, 1, 0, 0 n=(0,1,2,3,4,5)4.5.1理論分析1) 由DFT的性質(zhì)可以得到，如果y（n）的DFT為，那么y(n)就是x(n)圓周左移4位得到的。所以y（n）=0,0,4,3,2,12) 由題意，X(k)=10 3.5-j4.3301 2.5-8.66j 2 2.5+j0.86

13、6 3.5+j4.3301取實(shí)部，則W(k)=10 3.5 2.5 2.5 2 3.5 按照IDFT計(jì)算式計(jì)算得到w(n)=4 1.5 1 1 1.5 由DFT的性質(zhì)，因?yàn)閃(k)是實(shí)數(shù)，所以對(duì)應(yīng)的w(n)也是實(shí)數(shù)。3) 由題意，Q(k)= 10 2.5-8.66j 2.5+j0.866 按照IDFT計(jì)算式計(jì)算得到q(n)=5 3 24.5.2編程計(jì)算1)x=4 3 2 1 0 0;X=fft(x,6); %求DFTk=0:5;ifft(exp(j*4*pi/3*k).*X) %求IFFT 結(jié)果：y（n）=0,0,4,3,2,1 即x(n)圓周左移4位，與理論值相同。2)%接第（1）題的程序W

14、=real(X); %取實(shí)部ifft(W)結(jié)果：w(n)=4 1.5 1 1 1.5 與理論計(jì)算值相同。3) %接第（1）題的程序 Q=X(1) X(3) X(5); %Q是X（2k），由于matlab的矩陣是從X（1）開(kāi)始，ifft(Q) %所以對(duì)應(yīng)的應(yīng)該是第1、3、5個(gè)元素結(jié)果：q(n)=5 3 2 ，與理論值相同。4.6第六題已知信號(hào)，其中f1=4 Hz、f2=4.02 Hz、f3=5 Hz，采用采樣頻率為20 Hz進(jìn)行采樣，求:（1） 當(dāng)采樣長(zhǎng)度N分別為512和2048情況下x(t)的幅度頻譜；（2） 當(dāng)采樣長(zhǎng)度N為32，且增補(bǔ)N個(gè)零點(diǎn)、4N個(gè)零點(diǎn)、8N個(gè)零點(diǎn)、16N個(gè)零點(diǎn)情況下x(t

15、)的幅度頻譜。4.6.1理論分析1. 首先20Hz的采樣頻率是滿足采樣定理的。2. 頻率分辨率是DFT中譜線間的最小間隔，單位是Hz。對(duì)于長(zhǎng)度為N的序列，頻率分辨率為fs/N, 為采樣頻率。3. 因?yàn)椴蓸狱c(diǎn)數(shù)N不是周期的整數(shù)倍，所以一定會(huì)存在頻譜泄露情況。4. fs=20Hz時(shí),N=512,則分辨率 20/5120.039Hz 0.02Hz 所以不能區(qū)分開(kāi)信號(hào)中頻率為4Hz和4.02Hz的兩個(gè)分量。5. N=2048,則分辨率20/20480.01Hz0.02Hz 可以區(qū)分開(kāi)4Hz和4.02Hz的這兩個(gè)分量。4.6.2編程作圖編寫(xiě)一個(gè)取N個(gè)點(diǎn)并補(bǔ)充t*N個(gè)0的函數(shù)sample2(取樣點(diǎn)數(shù)N,補(bǔ)

16、零t)function x=sample2(N,t)n=0:0.02:0.02*(N-1); %N個(gè)點(diǎn)x=sin(2*pi*4*n)+sin(2*pi*4.02*n)+sin(2*pi*5*n); %取樣過(guò)程x(N+1:N*(t+1)=0; %補(bǔ)零end1) 輸入：x=sample2(N,0); %N=512 , 4096FFT1(x,N);結(jié)果：Figure 13 N=512幅頻圖N=512時(shí)，可以看到對(duì)4Hz和5Hz的分量是明顯區(qū)分開(kāi)來(lái)的。對(duì)峰值附近放大來(lái)看：Figure 14 N=512 幅頻圖放大看到k=103和104處有峰值。理論計(jì)算中， 由于頻譜泄露，在k=103和104處出現(xiàn)峰值

17、，實(shí)際上并沒(méi)有將4Hz和4.02Hz的兩個(gè)分量區(qū)分開(kāi)來(lái)。Figure 15 N=2048 幅頻圖N=2048時(shí)，顯然區(qū)分開(kāi)了4Hz和5Hz的分量。峰值附近放大看：在k=409、410和k=412有峰值。理論上 由于頻譜泄露，在409.6附近的409和410處同時(shí)出現(xiàn)峰值，在411.6附近的412出現(xiàn)峰值。所以采樣點(diǎn)數(shù)N=2048時(shí)一定程度上區(qū)分開(kāi)了4Hz和4.02Hz的分量。2) 輸入x=sample2(32,t); %t=取1,4,8,16FFT1(x,N); %N=32*(t+1)結(jié)果：Figure 16 補(bǔ)32個(gè)0 Figure 17 補(bǔ)128個(gè)0Figure 18 補(bǔ)256個(gè)0Figure 19 補(bǔ)512個(gè)0可以發(fā)現(xiàn)，采樣點(diǎn)數(shù)相同，都是32。20/32=0.61Hz，頻譜圖可以區(qū)分4Hz和5Hz分量，但不可以區(qū)分4Hz和4.02Hz的分量。補(bǔ)不同個(gè)數(shù)的0，幅頻圖的