直到現在才意識到，DFT或FFT無論是在訊號處理、還是在影象分析、又或是模式識別中的地位和重要性，以至於在大學竟然沒有自己寫過訊號的FFT的程式，果然出來混遲早是要還的，大學欠下的債現在該還了。花了一天看DFT、FFT，去證明蝶形運算、快速演算法，又花了一天去寫程式，去驗證Fourier的理論，去實現自己的想法。 1、Syntax and Description Y = fft(x) — returns the discrete Fourier Transform (DFT) of vector x,computed with a fast Fourier transform (FFT) algorithm. If the input x is a matrix, Y = fft(x) return the Fourier transform of each column of the matrix. Y = fft(x,n)

— returns the n-points DFT. fft(x) is equivalent to fft(x,n) where n is the size of x.If the length of x is less than n, x is padded with trailing zeros to length n. If the length of x is greater than n, the sequence x is truncated. When x is a matrix, the length of the columns are adjusted in the same manner. Y = fft(x,[ ],dim) and Y = fft(x,n,dim) — applies the FFT operation across the dimension dim. 2、程式示例

clc
clear all
close all

Fs=1000;						% Sampling Frequency
Ts=1/Fs;						% Sampling Time
f=50;							% Frequence of Signal
N=256;							% Length of Signal
t=(0:N-1)*Ts;					% Time Vector
x=sin(2*pi*f*t);				% Sine Signal
subplot(311)
plot(t,x);
ylabel('x(t)'),xlabel('t/s'),title('時域訊號')
grid on

NFFT=2^nextpow2(N);				% Next power of 2 from length of x.
y=fft(x,NFFT);            	    % 基-2快速傅立葉變換
X=abs(y);						% Compute Amplitude Spectrum of x
f=(0:N-1)/(N*Ts);         	    % Frency Vector
subplot(312)
plot(f,X);hold on
ylabel('X(f)'),xlabel('f/Hz'),title('頻域訊號')
grid on 
ind=find(X==max(X),1,'first');	%尋找最大值位置，返回1個linear inds：ind
x0=f(ind);						%根據位置得到橫座標（頻率）
y0=X(ind);						%根據位置得到縱座標（幅度）
plot(x0,y0,'ro');hold on
text(x0+1,y0-0.1,num2str(x0,'頻率=%f')); 
		%   text(x,y,'string')在圖形中指定的位置(x,y)上顯示字串string
ind=find(X==max(X),1,'last');	%尋找最大值位置，返回1個linear inds：ind
x0=f(ind);						%根據位置得到橫座標（頻率）
y0=X(ind);						%根據位置得到縱座標（幅度）
plot(x0,y0,'ro');hold off
text(x0-10,y0+10,num2str(x0,'頻率=%f')); 
		%   text(x,y,'string')在圖形中指定的位置(x,y)上顯示字串string

x1=ifft(y,NFFT);       			 % Inverse FFT
subplot(313)
plot(t,x1);
ylabel('x1(t)'),xlabel('t/s'),title('逆變換訊號')
grid on
 

3、程式解讀 先從理論分析sin（2Πft）的頻譜（幅度譜）長啥樣，根據理論有： 由上式可以得出結論：sin（2Πft）的FT是複數，我們接下來主要研究它的幅度譜。顯然，正弦函式的幅度譜就是在f0處和-f0處的兩條幅度都為（正）無窮大的衝激（不考慮衝激的相位，僅考慮模值）。 但是問題來了，程式始終不可能求連續的FT，計算機也只能求DFT或者FFT，所以訊號會被離散、截短，會將沒有訊號的位置認為是訊號的週期延拓。這會導致頻譜洩露，以至於頻譜幅度不會是無窮大，只能說是很大，中心頻率附近的頻率分得洩露的頻譜（能量）從而幅度不為0。 時域訊號以Fs進行取樣離散化，頻域會以Fs週期延拓；時域訊號以觀測時間NTs週期延拓，頻域將以1/（NTx）（頻域解析度）為間隔離散化。由於計算機執行的是DFT，我們只能得到主值序列【0到N-1的x（n）和X（k）】，所以在-f0處的譜線就會出現在Fs-f0。 接下來由兩個實驗加強理解傅立葉變換。

（1）改變取樣頻率Fs，取樣點數N=256，訊號頻率f0=50Hz。 a、Fs=1000Hz（訊號頻率的20倍），結果如下：

b、Fs=500Hz（訊號頻率的10倍），結果如下：

c、Fs=100Hz（Nyquist頻率），其實對於正弦訊號來說，在Nyquist取樣頻率下已經失真。結果如下：

（2）調整取樣點數N，f0=50Hz。 a、N=256，結果如下：

b、N=512，結果如下：

c、N=1024，結果如下：

雖然實驗不夠完善，但是也能粗略的說明問題。（1）中主要改變抽樣頻率，可以發現抽樣頻率必須要遠大於由取樣定理所確定的取樣頻率，那畢竟是“理論上”的極限，這樣才能有很高得保真度；（2）中主要改變訊號得觀測長度，取得得觀測點數N越多，取樣訊號與原始訊號（時間上無始無終）越接近，因為截短而造成得頻譜洩露程度就越輕。 完成！！！