數(shù)字信號處理實(shí)驗(yàn)：實(shí)驗(yàn)三快速傅立葉變換

1、實(shí)驗(yàn)三 快速傅立葉變換,一、實(shí)驗(yàn)?zāi)康?加深對離散傅立葉變換(DFT)的理解。 掌握利用MATLAB語言進(jìn)行離散傅立葉變換和逆變換的方法。 加深對離散傅立葉變換基本性質(zhì)的理解。 掌握離散傅立葉變換快速算法的應(yīng)用,二、實(shí)驗(yàn)原理及方法,有限長序列通過離散傅里葉變換(DFT)將其頻域離散化成有限長序列.但其計(jì)算量太大(與N的平方成正比）, 很難實(shí)時(shí)地處理問題, 因此引出了快速傅里葉變換(FFT)。 FFT并不是一種新的變換形式,它只是DFT的一種快速算法.并且根據(jù)對序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法,DFT的快速算法FFT是數(shù)字信號處理的基本方法和基本技術(shù)，是必須牢牢掌握的。 時(shí)間抽選F

2、FT算法的理論推導(dǎo)和流圖詳見數(shù)字信號處理教材。該算法遵循兩條準(zhǔn)則： （1）對時(shí)間奇偶分；（2）對頻率前后分。 這種算法的流圖特點(diǎn)是： （1）基本運(yùn)算單元都是蝶形 任何一個(gè)長度為N=2M的序列，總可通過M次分解最后成為2點(diǎn)的DFT計(jì)算。如圖所示,WNk稱為旋轉(zhuǎn)因子 計(jì)算方程如下： Xm+1(p)=Xm(p)+WNkXm(q) Xm+1(q)=Xm(p)-WNkXm(q,2）同址（原位）計(jì)算 這是由蝶形運(yùn)算帶來的好處，每一級蝶形運(yùn)算的結(jié)果 Xm+1(p)無須另外存儲，只要再存入Xm(p)中即可，Xm+1(q) 亦然。這樣將大大節(jié)省存儲單元。 （3）變址計(jì)算 輸入為“混序”（碼位倒置）排列，輸出按自

3、然序排 列，因而對輸入要進(jìn)行“變址”計(jì)算（即碼位倒置計(jì)算）。 “變址”實(shí)際上是一種“整序”的行為，目的是保證“同址,FFT的應(yīng)用,凡是利用付里葉變換來進(jìn)行分析、綜合、變換的地方，都可以利用FFT算法來減少其計(jì)算量。 FFT主要應(yīng)用在 1、快速卷積 2、快速相關(guān) 3、頻譜分析,快速傅立葉變換的MATLAB實(shí)現(xiàn),提供fft函數(shù)計(jì)算DFT 格式 X=fft(x) X=fft(x,N) 如果x的長度小于N，則在其后填零使其成為N點(diǎn)序列，反之對x進(jìn)行截?cái)?，若省略變量N，則DFT的長度即為x的長度。 如果N為2的冪，則得到高速的基-2FFT算法；若N不是2的乘方，則為較慢的混合算法。 如果x是矩陣，則X是

4、對矩陣的每一列向量作FFT,快速傅立葉逆變換（IFFT,函數(shù)調(diào)用格式 y=ifft(x) y=ifft(x,N) 當(dāng)N小于x長度時(shí)，對x進(jìn)行截?cái)?，?dāng)N大于x長度時(shí)，對x進(jìn)行補(bǔ)零,fftshift函數(shù),功能： 對fft的輸出進(jìn)行重新排列，將零頻分量移到頻譜的中心。 調(diào)用格式 y=fftshift(x) 當(dāng)x為向量時(shí)，fftshift(x)直接將x中左右兩半交換而產(chǎn)生y。 當(dāng)x為矩陣時(shí)，fftshift(x)直接將x中左右、上下進(jìn)行交換而產(chǎn)生y,由題目可得 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t) fs=100 N=128/1024,例：已知信號由15Hz幅值0.

5、5的正弦信號和40Hz幅值2的正弦信號組成，數(shù)據(jù)采樣頻率為100Hz,試?yán)L制N=128點(diǎn)DFT的幅頻圖,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),利用FFT進(jìn)行功率譜的噪聲分析,已知帶有測量噪聲信號 其中f1=50Hz,f2=120Hz, 為均值為零、方差為1的隨機(jī)信號，采樣頻率為1000Hz，數(shù)據(jù)點(diǎn)數(shù)N=512。試?yán)L制信號的功率譜圖

6、,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,序列長度和FFT的長度對信號頻譜的影響,已知信號 其中f1=15Hz,f2=40Hz,采樣頻率為100Hz. 在下列情況下繪制其幅頻譜。 Ndata=32,Nfft=32; Ndata=32,Nfft=128,fs=100;

7、 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,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

8、(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,線性卷積的FFT算法,在MATLAB實(shí)現(xiàn)卷積的函數(shù)為CONV，對于N值較小的向量，這是十分有效的。對于N值較大的向量卷積可用FFT加快計(jì)算速度。 由DFT性質(zhì)可知，若DFTx1(n)=X1(k),DFTx2(n)=X2(k)則 若DFT和IDFT均采用FFT和IFFT算法，可提高卷積速度,計(jì)算x1(n)和x2(n)的線性卷積的FFT算法可由下面步驟實(shí)現(xiàn),計(jì)算X1(k)=FFTx1(n); 計(jì)算X2(k)=FFTx2(n); 計(jì)算Y(k)=X1(k) X2(k); 計(jì)算x1(n)*x2(n)=IFFTY(k,用函數(shù)conv和FFT計(jì)算同一序列的卷積，比較其計(jì)算時(shí)間,clock函數(shù)讀取瞬時(shí)時(shí)鐘 etime(t1,t2)函數(shù)計(jì)算時(shí)刻t1,t2間所經(jīng)歷的時(shí)