快速傅立葉變換及其應(yīng)用

關(guān)于快速傅立葉變換及其應(yīng)用第1頁，講稿共28頁，2023年5月2日，星期三一、實(shí)驗(yàn)?zāi)康?/p>

了解計(jì)算DFT算法存在的問題及改進(jìn)途徑。掌握幾種DFT算法（時(shí)間抽取算法DIT算法，頻率抽取算法DIF算法，線性調(diào)頻Z變換即CZT法）。學(xué)習(xí)并掌握FFT的應(yīng)用。第2頁，講稿共28頁，2023年5月2日，星期三二、實(shí)驗(yàn)原理有限長(zhǎng)序列通過離散傅里葉變換(DFT)將其頻域離散化成有限長(zhǎng)序列.但其計(jì)算量太大(與N的平方成正比）,很難實(shí)時(shí)地處理問題,因此引出了快速傅里葉變換(FFT)。FFT并不是一種新的變換形式,它只是DFT的一種快速算法.并且根據(jù)對(duì)序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法.第3頁，講稿共28頁，2023年5月2日，星期三DFT的快速算法—FFT是數(shù)字信號(hào)處理的基本方法和基本技術(shù)，是必須牢牢掌握的。時(shí)間抽選FFT算法的理論推導(dǎo)和流圖詳見《數(shù)字信號(hào)處理》教材。該算法遵循兩條準(zhǔn)則：（1）對(duì)時(shí)間奇偶分；（2）對(duì)頻率前后分。這種算法的流圖特點(diǎn)是：（1）基本運(yùn)算單元都是蝶形

任何一個(gè)長(zhǎng)度為N=2M的序列，總可通過M次分解最后成為2點(diǎn)的DFT計(jì)算。如圖所示：第4頁，講稿共28頁，2023年5月2日，星期三WNk稱為旋轉(zhuǎn)因子計(jì)算方程如下：Xm+1(p)=Xm(p)+WNkXm(q)Xm+1(q)=Xm(p)-WNkXm(q)第5頁，講稿共28頁，2023年5月2日，星期三（2）同址（原位）計(jì)算這是由蝶形運(yùn)算帶來的好處，每一級(jí)蝶形運(yùn)算的結(jié)果Xm+1(p)無須另外存儲(chǔ)，只要再存入Xm(p)中即可，Xm+1(q)亦然。這樣將大大節(jié)省存儲(chǔ)單元。（3）變址計(jì)算輸入為“混序”（碼位倒置）排列，輸出按自然序排列，因而對(duì)輸入要進(jìn)行“變址”計(jì)算（即碼位倒置計(jì)算）?！白冎贰睂?shí)際上是一種“整序”的行為，目的是保證“同址”。第6頁，講稿共28頁，2023年5月2日，星期三FFT的應(yīng)用凡是利用付里葉變換來進(jìn)行分析、綜合、變換的地方，都可以利用FFT算法來減少其計(jì)算量。FFT主要應(yīng)用在1、快速卷積2、快速相關(guān)3、頻譜分析第7頁，講稿共28頁，2023年5月2日，星期三快速傅立葉變換的MATLAB實(shí)現(xiàn)提供fft函數(shù)計(jì)算DFT格式

X=fft(x)

X=fft(x,N)如果x的長(zhǎng)度小于N，則在其后填零使其成為N點(diǎn)序列，若省略變量N，則DFT的長(zhǎng)度即為x的長(zhǎng)度。如果N為2的冪，則得到高速的基-2FFT算法；若N不是2的乘方，則為較慢的混合算法。如果x是矩陣，則X是對(duì)矩陣的每一列向量作FFT。第8頁，講稿共28頁，2023年5月2日，星期三由題目可得x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)fs=100N=128/1024例：已知信號(hào)由15Hz幅值0.5的正弦信號(hào)和40Hz幅值2的正弦信號(hào)組成，數(shù)據(jù)采樣頻率為100Hz,試?yán)L制N=128點(diǎn)DFT的幅頻圖。第9頁，講稿共28頁，2023年5月2日，星期三fs=100;N=128;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);f=(0:length(y)-1)'*fs/length(y);mag=abs(y);stem(f,mag);title(‘N=128點(diǎn)’)第10頁，講稿共28頁，2023年5月2日，星期三第11頁，講稿共28頁，2023年5月2日，星期三利用FFT進(jìn)行功率譜的噪聲分析已知帶有測(cè)量噪聲信號(hào)其中f1=50Hz,f2=120Hz,為均值為零、方差為1的隨機(jī)信號(hào)，采樣頻率為1000Hz，數(shù)據(jù)點(diǎn)數(shù)N=512。試?yán)L制信號(hào)的頻譜圖和功率譜圖。第12頁，講稿共28頁，2023年5月2日，星期三t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(1,length(t));Y=fft(y,512);P=Y.*conj(Y)/512;%求功率f=1000*(0:255)/512;subplot(2,1,1);plot(y);subplot(2,1,2);plot(f,P(1:256));第13頁，講稿共28頁，2023年5月2日，星期三第14頁，講稿共28頁，2023年5月2日，星期三序列長(zhǎng)度和FFT的長(zhǎng)度對(duì)信號(hào)頻譜的影響。已知信號(hào)其中f1=15Hz,f2=40Hz,采樣頻率為100Hz.

在下列情況下繪制其幅頻譜。

Ndata=32,Nfft=32;Ndata=32,Nfft=128;第15頁，講稿共28頁，2023年5月2日，星期三fs=100;Ndata=32;Nfft=32;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)’*fs/length(y);subplot(2,1,1)plot(f(1:Nfft/2),mag(1:Nfft/2))title(‘Ndata=32,Nfft=32’)第16頁，講稿共28頁，2023年5月2日，星期三Nfft=128;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)’*fs/length(y);subplot(2,1,2)plot(f(1:Nfft/2),mag(1:Nfft/2))title(‘Ndata=32,Nfft=128’)第17頁，講稿共28頁，2023年5月2日，星期三第18頁，講稿共28頁，2023年5月2日，星期三快速傅立葉逆變換（IFFT）函數(shù)調(diào)用格式

y=ifft(x)y=ifft(x,N)當(dāng)N小于x長(zhǎng)度時(shí)，對(duì)x進(jìn)行截?cái)?，?dāng)N大于x長(zhǎng)度時(shí)，對(duì)x進(jìn)行補(bǔ)零。第19頁，講稿共28頁，2023年5月2日，星期三對(duì)信號(hào)進(jìn)行DFT，對(duì)其結(jié)果進(jìn)行IDFT，并將IDFT的結(jié)果和原信號(hào)進(jìn)行比較。f1=40Hzf2=15HzFs=100Hz第20頁，講稿共28頁，2023年5月2日，星期三fs=100;N=128;n=0:N-1;t=n/fs;x=sin(2*pi*40*t)+sin(2*pi*15*t);subplot(2,2,1)plot(t,x)title('originalsignal')y=fft(x,N);mag=abs(y);f=(0:length(y)-1)'*fs/length(y);subplot(2,2,2)plot(f,mag)title('FFTtooriginalsignal')第21頁，講稿共28頁，2023年5月2日，星期三xifft=ifft(y);magx=real(xifft);ti=[0:length(xifft)-1]/fs;subplot(2,2,3)plot(ti,magx);title('signalfromIFFT')yif=fft(xifft,N);mag=abs(yif);subplot(2,2,4)plot(f,mag)title('FFTtosignalfromIFFT')第22頁，講稿共28頁，2023年5月2日，星期三第23頁，講稿共28頁，2023年5月2日，星期三線性卷積的FFT算法在MATLAB實(shí)現(xiàn)卷積的函數(shù)為CONV，對(duì)于N值較小的向量，這是十分有效的。對(duì)于N值較大的向量卷積可用FFT加快計(jì)算速度。由DFT性質(zhì)可知，若DFT[x1(n)]=X1(k),DFT[x2(n)]=X2(n)則

若DFT和IDFT均采用FFT和IFFT算法，可提高卷積速度。第24頁，講稿共28頁，2023年5月2日，星期三計(jì)算x1(n)和x2(n)的線性卷積的FFT算法可由下面步驟實(shí)現(xiàn)計(jì)算X1(k)=FFT[x1(n)];計(jì)算X2(k)=FFT[x2(n)];計(jì)算Y(k)=X1(k)X2(k);計(jì)算x1(n)*x2(n)=IFFT[Y(k)].第25頁，講稿共28頁，2023年5月2日，星期三用函數(shù)conv和FFT計(jì)算同一序列的卷積，比較其計(jì)算時(shí)間。L=5000;N=L*2-1;n=1:L;x1=0.5*n;x2=2*n;t0=clock;yc=conv(x1,x2);conv_time=etime(clock,t0)t0=clock;yf=ifft(fft(x1,N).*fft(x2,N));fft_time