實(shí)驗(yàn)四 FFT的計(jì)算機(jī)實(shí)現(xiàn)及譜分析一、實(shí)驗(yàn)?zāi)康呐c要求(1) 了解DFT的性質(zhì)及應(yīng)用。(2) 熟悉MATLAB編程的特點(diǎn)。 (3)掌握有限離散傅里葉變換的計(jì)算機(jī)實(shí)現(xiàn)方法。()進(jìn)一步加深DFT算法原理和基本性質(zhì)的理解(因?yàn)镕FT只是DFT的一種快速算法，所以FPT的運(yùn)算結(jié)果必然滿足DFT的基本性質(zhì))。 ()熟悉FFT算法原理和FFT子程序的應(yīng)用。 ()學(xué)習(xí)用FFT對(duì)連續(xù)信號(hào)和時(shí)域離散信號(hào)進(jìn)行譜分析的方法，了解可能出現(xiàn)的分析誤差及其原因，以便在實(shí)際中正確應(yīng)用FFT。二、實(shí)驗(yàn)原理與方法.的計(jì)算原理設(shè)，有限離散信號(hào)，則它的有限離散頻譜為令則上式為把信號(hào)按下標(biāo)偶數(shù)項(xiàng)和奇數(shù)項(xiàng)分成兩部分，即令它們的有限離散頻譜分別為經(jīng)分析，我們可得如下的時(shí)域分解FFT算法若我們把按偶奇序號(hào)分成兩部分，則可得如下的頻域分解FFT算法、MATLAB提供的快速傅立葉變換函數(shù)（1）fftfft函數(shù)用于計(jì)算快速傅立葉變換，其語法格式為：B = fft(I)B = fft(I)返回信號(hào)I的fft變換矩陣，輸入信號(hào)I和輸出信號(hào)B大小相同。（2）fftshiftMATLAB提供的fftshift函數(shù)用于將變換后的信號(hào)頻譜中心從矩陣的原點(diǎn)移到矩陣的中心，其語法格式為：B = fftshift(I)對(duì)于矩陣I，B = fftshift(I)將I的一、三象限和二、四象限進(jìn)行互換。（3）ifftifft函數(shù)用于計(jì)算信號(hào)的傅立葉反變換，其語法格式為：B = ifft(I)B = ifft(A)返回信號(hào)I的傅立葉反變換矩陣，輸入信號(hào)I和輸出信號(hào)B大小相同。其語法格式含義與fft函數(shù)的語法格式相同，可以參考fft函數(shù)的說明。時(shí)域分解FFT算法的幾點(diǎn)說明：1、對(duì)一般的，如何求信號(hào)的一項(xiàng)頻譜排列？我們可通過二進(jìn)制的逆序得到：假設(shè)的二進(jìn)制為，則2、碟形運(yùn)算兩節(jié)點(diǎn)的距離：第一級(jí)（第一列）兩節(jié)點(diǎn)間的距離為1；第二級(jí)（第二列）兩節(jié)點(diǎn)間的距離為2；.第m級(jí)（第m列）兩節(jié)點(diǎn)間的距離為；3、的確定因子最后一列有種，順序?yàn)闉椋灰蜃拥箶?shù)第二列有種，為 ；因子倒數(shù)第三列有種，為；因子第一列有一種，為1頻域分解FFT算法的幾點(diǎn)說明：1、對(duì)一般的，最后得到的逆序頻譜如何排列才能得到原信號(hào)的頻譜？我們可通過二進(jìn)制的逆序得到：假設(shè)的二進(jìn)制為，則2、碟形運(yùn)算兩節(jié)點(diǎn)的距離：第一級(jí)（第一列，m1）兩節(jié)點(diǎn)間的距離為；.第m級(jí)（第m列）兩節(jié)點(diǎn)間的距離為；3、的確定因子第一列有種，順序?yàn)闉椋灰蜃拥诙杏蟹N，為 ；因子第三列有種，為；因子最后一列有一種，為1三 實(shí)驗(yàn)內(nèi)容及步驟 1、用三種不同的DFT程序計(jì)算的傅里葉變換,并比較三種程序計(jì)算機(jī)運(yùn)行時(shí)間。（1）用for 語句的M函數(shù)文件dft1.m，用循環(huán)變量逐點(diǎn)計(jì)算X(k)；（2）編寫時(shí)域分解FFT 算法或頻域分解FFT 算法的M函數(shù)文件dft2.m；（3）調(diào)用FFT庫(kù)函數(shù)，直接計(jì)算X(k)；（4）分別利用上述三種不同方式編寫的DFT程序計(jì)算序列x(n)的傅立葉變換，并畫出相應(yīng)的幅頻和相頻特性，再比較各個(gè)程序的計(jì)算機(jī)運(yùn)行時(shí)間。()復(fù)習(xí)DFT的定義、性質(zhì)和用DFT作譜分析的有關(guān)內(nèi)容。 ()復(fù)習(xí)FFT算法原理與編程思想，并對(duì)照DITFFT運(yùn)算流圖和程序框圖，讀懂本實(shí)驗(yàn)提供的FFT子程序。()編制信號(hào)產(chǎn)生子程序，產(chǎn)生以下典型信號(hào)供譜分析用： ()編寫主程序。 ()按實(shí)驗(yàn)內(nèi)容要求，上機(jī)實(shí)驗(yàn)，并寫出實(shí)驗(yàn)報(bào)告。M函數(shù)文件如下： dft1.m:functionAm,pha=dft1(x)N=length(x);w=exp(-j*2*pi/N);for k=1:N sum=0; for n=1:N sum=sum+x(n)*w(k-1)*(n-1); end Am(k)=abs(sum); pha(k)=angle(sum);enddft3.m:functionAm,pha=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk); 11源程序及運(yùn)行結(jié)果：（1）N=256;x=ones(1,8),zeros(1,N-8);t=cputime;Am1,pha1=dft1(x);t1=cputime-tn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(1)subplot(2,1,1), plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2), plot(w,pha1,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians); 運(yùn)行時(shí)間：t1=0.288s（3）N=256;x=ones(1,8),zeros(1,N-8);t=cputime;Am3,pha3=dft3(x);t3=cputime-t;n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(3)subplot(2,1,1), plot(w,Am3,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);subplot(2,1,2), plot(w,pha3,r); grid;title(Phase Part);xlabel(frequency in radians);ylabel(argXexp(jw)/radians);運(yùn)行時(shí)間：t3=0從以上運(yùn)行結(jié)果可以看出，調(diào)用FFT庫(kù)函數(shù)直接計(jì)算X(k)速度最快，用循環(huán)變量逐點(diǎn)計(jì)算運(yùn)行速度最慢。因此FFT算法大大提高了DFT的實(shí)用性。四、研究高密度頻譜和高分辨率頻譜設(shè)有連續(xù)信號(hào)Xa(t)=cos(26.5103t)+cos(27103t)+cos(29103t)以采樣頻率fs=32kHz對(duì)信號(hào)Xa(t)采樣,分析下列三種情況的幅頻特性。（1） 采集數(shù)據(jù)長(zhǎng)度N=16點(diǎn)，做N=16點(diǎn)的DFT，并畫出幅頻特性。（2） 采集數(shù)據(jù)長(zhǎng)度N=16點(diǎn)，補(bǔ)零到256點(diǎn)，做N=256點(diǎn)的DFT，并畫出幅頻特性。（3） 采集數(shù)據(jù)長(zhǎng)度N=256點(diǎn)，做N=256點(diǎn)的DFT，并畫出幅頻特性。觀察三幅不同頻率特性圖，分析和比較它們的特點(diǎn)以及形成的原因。源程序如下：fs=32000; %采樣頻率N=16; %采集16點(diǎn)n=0:N-1; %做16點(diǎn)DFTxa=cos(2*pi*6.5*103*n/fs)+cos(2*pi*7*103*n/fs)+cos(2*pi*9*103*n/fs);Am1,pha1=dft3(xa);n=0:(length(xa)-1);w=(2*pi/length(xa)*n;figure(1)plot(w,Am1,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);x=xa(1:16),zeros(1,240); %補(bǔ)零到256Am2,pha2=dft2(x); %做256點(diǎn)DFTn=0:(length(x)-1);w=(2*pi/length(x)*n;figure(2)plot(w,Am2,b); grid;title(Magnitude part);xlabel(frequency in radians);ylabel(|X(exp(jw)|);N0=256; %采集256點(diǎn)n0=0:N0-1; %做256點(diǎn)DFTxa0=cos(2*pi*6.5*103*n0/fs)+cos(2*pi*7*103*n0/fs)+cos(2*pi*9*103*n0/fs);Am3,pha3=dft3(xa0);n=0:(length(xa0)-1);w=(2*pi/length(xa0)*n;figure(3)plot(w,Am3,b); grid;title(Magnitude part);xlabel(frequency in radians); ylabel(|X(exp(jw)|); 運(yùn)行結(jié)果：N=16點(diǎn)DFT幅頻特性 高密度頻譜： 高分辨率頻譜：觀察運(yùn)行結(jié)果可知，N=16點(diǎn)時(shí)的DFT，由于點(diǎn)數(shù)太少，很難反映出頻譜的細(xì)節(jié)特征，只能分辨出兩個(gè)頻率分量。將16點(diǎn)DFT補(bǔ)零到256點(diǎn)，做256點(diǎn)DFT,由結(jié)果可以看出，頻率間隔縮小（由2p/16減小為2 p/256），得到了比較平滑的連續(xù)曲線，此為高密度頻譜。它是由16點(diǎn)經(jīng)過補(bǔ)零所得到的，并沒有增加任何信息，頻率分辨率并未提高，仍然只能看出兩個(gè)頻率分量。第三次做的是256點(diǎn)的DFT，采集點(diǎn)數(shù)大大增加，分辨率提高，可以清楚看到三根譜線的位置，是為高分辨率頻譜。五、上機(jī)實(shí)驗(yàn)內(nèi)容(1)對(duì)三中所給出的信號(hào)逐個(gè)進(jìn)行譜分析。 (2)令，用FFT計(jì)算8點(diǎn)和16點(diǎn)離散傅里葉變換。 (3)令，重復(fù)（2）。(4)對(duì)第四部分給以實(shí)現(xiàn)并進(jìn)行說明。六、思考題（1）在N=8時(shí)，和的幅頻