Matlab中fft與fftshift命令的小結與分析

http://wang-yg.diandian.com/post/2011-03-12/40028916801

註：轉載請註明出處——by author.

我們知道Fourier分析是訊號處理裡很重要的技術，matlab提供了強大的訊號處理能力，但是有一些細節部分需要我們注意。
記訊號f(t)的起始時間為t_start, 終止時間為t_end, 採樣周期為t_s, 可以計算訊號的期間Duration為 t_end – t_start, 訊號離散化造成的採樣點數 N = Duration/t_s + 1;
根據Fourier分析的相關結論，我們知道時域的採樣將會造成頻域的周期化，該周期為採樣頻率f_s（著名的香農採樣定理基於此）.
於是， 經過matlab的fft函數處理後，得到資料的橫座標為0:f_s/(N-1):f_s。相關代碼如下所示：
%matlab fft 測試代碼
t_s = 0.01;
t_start = 0.5; t_end = 5;
t = t_start:t_s:t_end;
y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y_f = fft(y);
subplot(3,1,1);
plot(t,y); title('original signal');
Duration = t_end - t_start;
Sampling_points = Duration/t_s + 1;
f_s = 1/t_s;
f_x = 0:f_s/(Sampling_points-1):f_s;
subplot(3,1,2);
plot(f_x,abs(y_f)); title('fft transform');
subplot(3,1,3);
plot(f_x-f_s/2,abs(fftshift(y_f))); title('shift fft transform');
也就是說，如果我們不使用fftshift，其變換後的橫座標為0:f_s/(N-1):f_s， 如果使用fftshift命令，0頻率分量將會移到座標中心，這也正是matlab中協助中心給出的意思：對fft的座標進行了處理。實際上由於頻譜的周期性，我們這樣做是合理的，可以接受的。

請讀者特別要注意橫座標的差別。另外，根據函數的特性，頻譜應當只有在15Hz，40Hz出現峰值，但是fft變換後在60Hz，及85Hz處同樣出現了峰值，應當可以從fft的計算過程中得到相應的解釋。

事實上，如果我們用15Hz，60Hz來測試fft變換，也即是 y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*60*t);映像如下所示，沒有任何變化。

這種現象提醒我們，頻率在f_s以內，即 0<f<f_s，f 以及 f_s – f 都有可能是測試訊號的頻率譜，這就給我們帶來了歧義。並且從第三個子圖也可以看出，這時候的fftshift會給我們帶來錯誤的引導，也就是說，如果我們試圖採樣fft或者fftshift來分析訊號的頻率譜顯得不那麼靠譜了，matlab的fft譜線與訊號的實際頻率並不是一一對應的映射關係。這當然不是我們所期望看到的結果，所以實際分析訊號時，有關這個問題需要額外的注意。
實際上，這也就間接地證明了Nyquist採樣定理的合理性：採樣頻率要高於截止頻率的兩倍，上面的處理中我們所使用的採樣頻率為100Hz，於是當截止頻率超過50Hz時，就會出現混疊效應，特殊情況就如所示：完全一樣。於是，這也就告訴我們若要正確的顯示頻譜，需要仔細地考量採樣頻率與截止頻率的關係，若太小，則有可能出現混疊，若太大，則計算代價過高。