MATLAB中fft函數(shù)用法、性質(zhì)、特性、缺陷全面深入解析(含程序)

MATLAB中fft函數(shù)用法、性質(zhì)、特性、缺陷全面深入解析〔含程序〕前言目前網(wǎng)絡(luò)上在信號(hào)處理方面對(duì)fft函數(shù)使用泛濫，但大多數(shù)人根本就不了解它的性能，使用方法也是千奇百怪，各種程序代碼魚(yú)龍混雜，亂七八糟，誤導(dǎo)了很多人，所以我在這里進(jìn)行一下全面的演示說(shuō)明，希望后來(lái)人能夠少走彎路。fft信號(hào)處理中獲得時(shí)域信號(hào)的頻域的方法叫做傅里葉變換，但教材上的變換都是理論推導(dǎo)的理想、連續(xù)、無(wú)限情況，在實(shí)際的matlab計(jì)算中，必須是實(shí)際、離散、有限的計(jì)算，所以實(shí)現(xiàn)起來(lái)比擬麻煩，需要各種轉(zhuǎn)換。fft函數(shù)就是為了實(shí)現(xiàn)快速傅里葉變換而設(shè)計(jì)的計(jì)算機(jī)算法，與傳統(tǒng)傅里葉變換的步驟不完全一樣，但通過(guò)適當(dāng)?shù)淖儞Q也能夠獲得時(shí)域信號(hào)的頻域圖形。 fft只是數(shù)值計(jì)算，不是理論推導(dǎo)，所以頻域圖形肯定有誤差，但我們能夠通過(guò)適當(dāng)調(diào)整變換參數(shù)來(lái)提高頻域圖形的精確性。建議大家能夠理論推導(dǎo)的傅里葉變換就直接在紙上推導(dǎo)好了再用matlab直接畫(huà)出來(lái)就是了。所以不是因?yàn)閒ft精確，而是因?yàn)閒ft能夠?qū)Υ髷?shù)據(jù)復(fù)雜信號(hào)變換計(jì)算快，自由度可控制性高，所以使用的人才多。fft(x)圖示fft的一些定性問(wèn)題解析時(shí)域信號(hào)x有多少個(gè)點(diǎn)，fft(x)就有多少個(gè)點(diǎn)，這時(shí)的fft(x)之所以是對(duì)稱的，是因?yàn)閒ft的方法是進(jìn)行周期計(jì)算從-2π,2π，所以后面的一半要翻轉(zhuǎn)到負(fù)半軸，是負(fù)頻率。所以有用的是fft(x)的前一半的點(diǎn)。我是經(jīng)過(guò)各項(xiàng)測(cè)試實(shí)驗(yàn)，才完全明白fft的性能。利用程序解釋如下：%這里的fft變換有補(bǔ)零clearall;Fs=1e3;%這是原始時(shí)域信號(hào)x的采樣頻率multiple=0.6;%x進(jìn)行fft前的補(bǔ)零倍數(shù)，可以增加頻率分辨率time=2;t=0:1/Fs:time-1/Fs;%這個(gè)會(huì)影響x的采樣點(diǎn)數(shù)，當(dāng)然點(diǎn)數(shù)越多，頻域越精確越平滑%x=1.5*cos(2*pi*100*t)+3*sin(2*pi*202.5*t)+2*sin(2*pi*1250*t);x=1.5*cos(2*pi*100*t)+3*sin(2*pi*201.25*t);%頻率最小分辨率————————————————-x=[x,zeros(1,length(x)*multiple)];%{注意可以通過(guò)對(duì)x進(jìn)行補(bǔ)零來(lái)增加頻率分辨率，只是因?yàn)閒ft的周期性計(jì)算方法，需要在后面對(duì)fft(x)后的數(shù)據(jù)進(jìn)行幅度值重新計(jì)算時(shí)，乘以多的零點(diǎn)的倍數(shù)，見(jiàn)后面的幅度值計(jì)算。%}%freqres=Fs/length(x);freqres=1/time;%{x的頻率必須是freqres的整數(shù)倍時(shí)，x的頻率才能被正確分辨，也就是相應(yīng)頻域值才能夠正常，比方freqres=1,在頻率202.5處的幅度就不正常，當(dāng)freqres=2時(shí)，在頻率1001處的幅度就不正常。%}%——————————————————%{x有多少個(gè)點(diǎn)，fft(x)就有多少個(gè)點(diǎn)，這時(shí)的fft(x)之所以是對(duì)稱的，-是因?yàn)閒ft的方法是進(jìn)行周期計(jì)算從[-2π,2π]，所以后面的一半要翻轉(zhuǎn)到負(fù)半軸，是負(fù)頻率。所以有用的是fft(x)的前一半的點(diǎn)再加上0點(diǎn)。%}xdft=fft(x);xdft=xdft(1:length(x)/2+1);%————這與fft方法的周期計(jì)算方法相關(guān)，也就是幅度值要重新計(jì)算一下——————xdft=1/length(x).*xdft;xdft(2:end-1)=(1+multiple)*2*xdft(2:end-1);%{幅度值首先由于周期對(duì)稱性，需要乘以2，然后由于補(bǔ)零的倍數(shù)是multiple，所以是(1+multiple)*2%}%重新計(jì)算頻率橫坐標(biāo)—————————————————————————————freq=0:Fs/length(x):Fs/2;%{Fs/2就是最高能夠顯示出來(lái)的頻率，但要注意實(shí)際上，fft(x)實(shí)際是計(jì)算出了任何頻率，但是在[Fs/2,Fs]這個(gè)區(qū)間的頻率，是關(guān)于Fs/2這里對(duì)稱的，也就是正負(fù)的頻率開(kāi)始混合或者混疊在一起。而超出了Fs的頻率那么會(huì)由于周期性，它的值先減去Fs后，重新從0點(diǎn)開(kāi)始。所以這個(gè)程序并不是說(shuō)畫(huà)出了在[0,Fs/2]的一定正確的頻率，而是可能會(huì)出現(xiàn)頻率的混合和混疊，一定要先預(yù)判一下頻率的區(qū)間范圍后，再使用fft這個(gè)函數(shù)。%}%——————————————————————————plot(freq,abs(xdft)); %畫(huà)頻域幅度譜title('頻域圖');fft的一些定量詳細(xì)解析與測(cè)試?yán)贸绦蚪忉屓缦拢?為了更好的展示，這里信號(hào)沒(méi)有補(bǔ)零clearall;clc;Fs=1e3;%這是原始時(shí)域信號(hào)x的采樣頻率multiple=0;%x進(jìn)行fft前的補(bǔ)零倍數(shù)，可以增加頻率分辨率，也會(huì)影響信號(hào)頻域值time=10;t=0:1/Fs:time-1/Fs;%這個(gè)會(huì)影響x的采樣點(diǎn)數(shù)，當(dāng)然點(diǎn)數(shù)越多，頻域越精確越平滑%x=1.5*cos(2*pi*100*t)+3*sin(2*pi*202.5*t)+2*sin(2*pi*1450*t);%x=1.5*cos(2*pi*100*t)+3*sin(2*pi*200*t)+2*cos(2*pi*100000.1*t);x=1.5*cos(2*pi*100*t)+3*sin(2*pi*200*t)+2*cos(2*pi*0.1*t);%頻率最小分辨率————————————————-x=[x,zeros(1,length(x)*multiple)];%{注意可以通過(guò)對(duì)x進(jìn)行補(bǔ)零來(lái)增加頻率分辨率，只是因?yàn)閒ft的周期性計(jì)算方法，需要在后面對(duì)fft(x)后的數(shù)據(jù)進(jìn)行幅度值重新計(jì)算時(shí)，乘以多的零點(diǎn)的倍數(shù)，見(jiàn)后面的幅度值計(jì)算。%}freqres=Fs/length(x);%{x的頻率必須是freqres的整數(shù)倍時(shí)，x的頻率才能被正確分辨，也就是幅度值才能夠正常，比方freqres=1,在頻率202.5處的幅度就不正常，當(dāng)freqres=2時(shí)，在頻率1001處的幅度就不正常。%}%——————————————————%{有多少個(gè)點(diǎn)，fft(x)就有多少個(gè)點(diǎn)，這時(shí)的fft(x)之所以是對(duì)稱的，-是因?yàn)閒ft的方法是進(jìn)行周期計(jì)算從[-2π,2π]，所以后面的一半要翻轉(zhuǎn)到負(fù)半軸，是負(fù)頻率。所以有用的是fft(x)的前一半的點(diǎn)再加上0點(diǎn)。%}xdft=fft(x);subplot(3,1,1);plot(0:length(xdft)-1,abs(xdft));title('abs');title('頻域圖');subplot(3,1,2);plot(0:length(xdft)-1,real(xdft));title('real');subplot(3,1,3);plot(0:length(xdft)-1,imag(xdft));title('imag');back=ifft(xdft);error=x-back;fori=2:length(xdft)/2store(i)=abs(xdft(i))-abs(xdft(length(xdft)-i+2));%這時(shí)的誤差為0,%store(i)=real(xdft(i+1))-real(xdft(length(xdft)-i+1));%這時(shí)的誤差為0,%store(i)=imag(xdft(i+1))+imag(xdft(length(xdft)-i+1));%這時(shí)的誤差為0,注意這時(shí)是加號(hào)endfigure;plot(store);title('正負(fù)頻率誤差');首先如上面程序所示，xdft=fft(x)，頻率分辨率是freqres，設(shè)L=length(xdft)=length(x)，顯示的頻率范圍是0:freqres:Fs-freqres，點(diǎn)數(shù)范圍是1，L。xdft的幅度譜在頻率范圍1*freqres，L2-1*freqres，也就是點(diǎn)數(shù)范圍2，L2內(nèi)，是與右邊點(diǎn)數(shù)范圍L，L2+2內(nèi)的頻率的實(shí)部是完全對(duì)稱的，虛部是正負(fù)對(duì)稱的，所以幅度譜是完全對(duì)稱的。xdft的幅度譜在頻率范圍L2+1*freqres，L-1*freqres，也就是點(diǎn)數(shù)范圍L2+2，L，是與左邊邊點(diǎn)數(shù)范圍L2，2內(nèi)的頻率的實(shí)部是完全對(duì)稱的，虛部是正負(fù)對(duì)稱的，所以幅度譜是完全對(duì)稱的，但是這時(shí)候左右的正負(fù)頻率是混疊的，這就是奈奎斯特準(zhǔn)那么fft變換存在的缺陷與問(wèn)題有人問(wèn)，為什么我畫(huà)的頻域圖與書(shū)上的不一樣啊，這就是你在對(duì)參數(shù)進(jìn)行設(shè)置時(shí)，沒(méi)有考慮到它們對(duì)于fft變換后的結(jié)果的影響，或者根本就沒(méi)設(shè)置正確。不要直接就進(jìn)行fft變換，變換前必須先估計(jì)信號(hào)的頻率范圍，和你想要獲得的圖形精確程度，然后適當(dāng)改變參數(shù)，確定不會(huì)發(fā)生下面這些問(wèn)題。1.混疊效應(yīng)如果x(t)的頻譜是帶限的，X(f)=0,|f|>fm，那么由抽樣定理，抽樣間隔最大滿足Ts=1/2fm如果f(t)的頻譜不是帶限的，那么抽樣后頻譜總要發(fā)生混疊，減小抽樣間隔Ts，fs增大，可減小混疊，但工作量增加。解決方法：預(yù)濾波，再抽樣，或者一般選擇Ts<1/(3~5)fm。2.泄漏(leakage)假設(shè)X(f)為有限帶寬頻譜，那么x(t)為時(shí)間無(wú)限的。為利用FFT分析x(t)的頻譜，必須截取x(t)有限范圍，即加窗。頻域卷積后，使原頻帶受限的頻譜擴(kuò)展開(kāi)來(lái)(有限帶寬拖了尾巴)，這種現(xiàn)象稱為泄漏。解決方法：改善窗的形狀或長(zhǎng)度。3.柵欄效應(yīng)〔PicketFenceEffect〕N點(diǎn)序列的DFT只能在有限的N個(gè)頻點(diǎn)上觀察頻譜，這相當(dāng)于從柵欄的縫隙中觀察景色，對(duì)于了解信號(hào)在整個(gè)頻域上的特性是不夠的。為了觀察到其他頻率的信息，需要對(duì)原信號(hào)x[n]做一些處理，以便在不同的頻點(diǎn)上采樣。解決方法：增加采樣點(diǎn)數(shù)，信號(hào)補(bǔ)零。4.補(bǔ)零與fft分辨率通過(guò)補(bǔ)零觀察到更多的頻點(diǎn)，但是這并不意味著補(bǔ)零能夠提高真正的頻譜分辨率。這是因?yàn)閤[n]實(shí)際上是x(t)采樣的主值序列，而將x[n]補(bǔ)零得到的x'[n]周期延拓之后與原來(lái)的序列并不相同，也不是x(t)的采樣。因此與是不同離散信號(hào)的頻譜。對(duì)于補(bǔ)零至M點(diǎn)的x'的DFT，只能說(shuō)它的分辨率2π/M僅具有計(jì)算上的意義，并不是真正的、物理意義上的頻譜。物理頻譜分辨率的提高只能在滿足采樣定理的條件下增加時(shí)域采樣長(zhǎng)度來(lái)實(shí)現(xiàn)。5.物理分辨率與計(jì)算分辨率頻率分辨率是指所用的算法能將信號(hào)中兩個(gè)考得很近的譜峰保持分開(kāi)的能力。