來源：學步園
FFT（Fast Fourier Transform，快速傅立葉變換）是離散傅立葉變換的快速算法，也是我們在數字信號處理技術中經常會提到的一個概念。在大學的理工科課程中，在完成高等數學的課程后，數字信號處理一般會作為通信電子類專業的專業基礎課程進行學習，原因是其中涉及了大量的高等數學的理論推導，同時又是各類應用技術的理論基礎。
關於傅立葉變換的經典著作和文章非常多，但是看到滿篇的復雜公式推導和羅列，我們還是很難從直觀上去理解這一復雜的概念，我想對於普通的測試工程師來說，掌握FFT的概念首先應該搞清楚這樣幾個問題：(1) 為什么需要FFT (2) 變換究竟是如何進行的 (3) 變換前后信號有何種對應關系。
在這篇文章中我嘗試用更加淺顯的講解，盡量不使用公式推導來說一說FFT的那些事兒。
為什么需要FFT？
FFT(快速傅立葉變換)是離散傅立葉變換的快速算法.
傅立葉變換的物理意義在哪里？
傅立葉原理表明：任何連續測量的時序或信號，都可以表示為不同頻率的正弦波信號的無限疊加。而根據該原理創立的傅立葉變換算法利用直接測量到的原始信號，以累加方式來計算該信號中不同正弦波信號的頻率、振幅和相位。當然這是從數學的角度去看傅立葉變換。
那么從物理的角度去看待傅立葉變換，它其實是幫助我們改變傳統的時間域分析信號的方法轉到從頻率域分析問題的思維，下面的一幅立體圖形可以幫助我們更好得理解這種角度的轉換：
所以，最前面的時域信號在經過傅立葉變換的分解之后，變為了不同正弦波信號的疊加，我們再去分析這些正弦波的頻率，可以將一個信號變換到頻域。有些信號在時域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征了。這就是很多信號分析采用FFT變換的原因。另外，FFT可以將一個信號的頻譜提取出來，這在頻譜分析方面也是經常用的。
傅立葉變換提供給我們這種換一個角度看問題的工具，看問題的角度不同了，問題也許就迎刃而解！
變換是如何進行的？
首先，按照被變換的輸入信號類型不同，傅立葉變換可以分為4種類型：
1、 非周期性連續信號傅立葉變換（Fourier Transform）
2、 周期性連續信號傅立葉級數(Fourier Series)
3、 非周期性離散信號離散時域傅立葉變換（Discrete Time Fourier Transform）
4、 周期性離散信號離散傅立葉變換(Discrete Fourier Transform)
下面是四種原信號圖例：
這里我們要討論是離散信號，對於連續信號我們不作討論，因為計算機只能處理離散的數值信號，我們的最終目的是運用計算機來處理信號的。所以對於離散信號的變換只有離散傅立葉變換（DFT）才能被適用，對於計算機來說只有離散的和有限長度的數據才能被處理，對於其它的變換類型只有在數學演算中才能用到，在計算機面前我們只能用DFT方法，我們要討論的FFT也只不過是DFT的一種快速的算法。
DFT的運算過程是這樣的：
其中，
X(k)—頻域值
X(n)—時域采樣點
n—時域采樣點的序列索引
k—頻域值的索引
N—進行轉換的采樣點數量
可見，在計算機或者示波器上進行的DFT，使用的輸入值是數字示波器經過ADC后采集到的采樣值，也就是時域的信號值，輸入采樣點的數量決定了轉換的計算規模。變換后的頻譜輸出包含同樣數量的采樣點，但是其中有一半的值是冗余的，通常不會顯示在頻譜中，所以真正有用的信息是N/2+1個點。
FFT的過程大大簡化了在計算機中進行DFT的過程，簡單來說，如果原來計算DFT的復雜度是N2次運算（N代表輸入采樣點的數量），進行FFT的運算復雜度是Nlg10（N），因此，計算一個1,000采樣點的DFT，使用FFT算法只需要計算3,000次，而常規的DFT算法需要計算1,000,000次！
我們以一個4個點的DFT變換為例來簡單說明FFT是怎樣實現快速算法的：
計算得出：
其中的紅色部分在FFT中是必須計算的分量，其他藍色部分不需要直接計算，可以由紅色的分量直接推導得到，比如：
x(1)e-j0 = -1*x(1)e-jπ
x(2)e-j0 = x(2)e-j2π
… …
這樣，已經計算出的紅色分量只需要計算機將結果保存下來用於之后計算時調用即可，因此大大減少了DFT的計算量。
變換前后信號有何種對應關系？
我們以一個實際的信號為例來說明：
示波器采樣得到的數字信號，就可以做FFT變換了。N個采樣點，經過FFT之后，就可以得到N個點的FFT結果。為了方便進行FFT運算，通常N取2的整數次方。
假設采樣頻率為Fs，信號頻率F，采樣點數為N。那么FFT之后結果就是一個為N點的復數。每一個點就對應着一個頻率點。這個點的模值，就是該頻率值下的幅度特性。具體跟原始信號的幅度有什么關系呢？假設原始信號的峰值為A，那么FFT的結果的每個點（除了第一個點直流分量之外）的模值就是A的N/2倍。而第一個點就是直流分量，它的模值就是直流分量的N倍。而每個點的相位呢，就是在該頻率下的信號的相位。第一個點表示直流分量（即0Hz），而最后一個點N的再下一個點（實際上這個點是不存在的，這里是假設的第N+1個點，也可以看做是將第一個點分做兩半分，另一半移到最后）則表示采樣頻率Fs，這中間被N-1個點平均分成N等份，每個點的頻率依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到頻率為為Fs/N，如果采樣頻率Fs為1024Hz，采樣點數為1024點，則可以分辨到1Hz。1024Hz的采樣率采樣1024點，剛好是1秒，也就是說，采樣1秒時間的信號並做FFT，則結果可以分析精確到1Hz，如果采樣2秒時間的信號並做FFT，則結果可以分析精確到0.5Hz。如果要提高頻率分辨率，則必須增加采樣點數，也即采樣時間。頻率分辨率和采樣時間是倒數關系。
下面這幅圖更能夠清晰地表示這種對應關系：
變換之后的頻譜的寬度（Frequency Span）與原始信號也存在一定的對應關系。根據Nyquist采樣定理，FFT之后的頻譜寬度（Frequency Span）最大只能是原始信號采樣率的1/2，如果原始信號采樣率是4GS/s，那么FFT之后的頻寬最多只能是2GHz。時域信號采樣周期（Sample Period）的倒數，即采樣率（Sample Rate）乘上一個固定的系數即是變換之后頻譜的寬度，即 Frequency Span = K*（1/ΔT），其中ΔT為采樣周期，K值取決於我們在進行FFT之前是否對原始信號進行降采樣（抽點），因為這樣可以降低FFT的運算量。如下圖所示：
可見，更高的頻譜分辨率要求有更長的采樣時間，更寬的頻譜分布需要提高對於原始信號的采樣率，當然我們希望頻譜更寬，分辨率更精確，那么示波器的長存儲就是必要的！它能提供在高采樣率下采集更長時間信號的能力。
頻譜泄露
所謂頻譜泄露，就是信號頻譜中各譜線之間相互干擾，使測量的結果偏離實際值，同時在真實譜線的兩側的其它頻率點上出現一些幅值較小的假譜。產生頻譜泄露的主要原因是采樣頻率和原始信號頻率不同步，造成周期的采樣信號的相位在始端和終端不連續。簡單來說就是因為計算機的FFT運算能力有限，只能處理有限點數的FFT，所以在截取時域的周期信號時，沒有能夠截取整數倍的周期。信號分析時不可能取無限大的樣本。只要有截斷不同步就會有泄露。如下圖所示：
圖中被測信號的開始端相位和截止端相位相同，表示在采集時間內有整數倍周期的信號被采集到，所以此時經行FFT運算后得出的頻譜不會出現泄露。
上圖的信號頻率為2.1MHz，采集時間內沒有截取整數倍周期的信號，FFT運算之后譜線的泄露現象嚴重，可以看到能量較低的譜線很容易被臨近的能量較高的譜線的泄露給淹沒住。
因此，避免頻譜泄露的方法除了盡量使采集速率與信號頻率同步之外，還可以采用適當的窗函數。
另外一個方法是采集信號時間足夠長，基本上可以覆蓋到整個有效信號的時間跨度。這種方法經常在瞬態捕捉中被使用到，比如說沖擊試驗，如果捕捉的時間夠長，捕捉到的信號可以一直包括了振動衰減為零的時刻。在這種情況下，可以不加窗函數。
窗函數其實就是一個加權函數，它在截取的信號時間段內有值，時間段之外值為0:，記為：
w(t)=g(t) -T/2<t<T/2
w(t)=0 其它
加窗在時域上表現的是點乘，因此在頻域上則表現為卷積。卷積可以被看成是一個平滑的過程。這個平滑過程可以被看出是由一組具有特定函數形狀的濾波器，因此，原始信號中在某一頻率點上的能量會結合濾波器的形狀表現出來，從而減小泄漏。基於這個原理，人們通常在時域上直接加窗。
大多數的信號分析儀一般使用矩形窗（rectangular），漢寧（hann），flattop和其它的一些窗函數。
不同的窗函數對頻譜譜線的影響不同，基本形狀可以參看下圖：
可以看到，不同的窗函數的主瓣寬度和旁瓣的衰減速度都不一樣，所以對於不同信號的頻譜應該使用適當的窗函數進行處理。
矩形窗(Rectangular)：加矩形窗等於不加窗，因為在截取時域信號時本身就是采用矩形截取，所以矩形窗適用於瞬態變化的信號，只要采集的時間足夠長，信號寬度基本可以覆蓋整個有效的瞬態部分。
漢寧窗(Von Hann)：如果測試信號有多個頻率分量，頻譜表現的十分復雜，且測試的目的更多關注頻率點而非能量的大小。在這種情況下，需要選擇一個主瓣夠窄的窗函數，漢寧窗是一個很好的選擇。
flattop窗：如果測試的目的更多的關注某周期信號頻率點的能量值，比如，更關心其EUpeak,EUpeak-peak,EUrms，那么其幅度的准確性則更加的重要，可以選擇一個主瓣稍寬的窗，flattop窗在這樣的情況下經常被使用。
好了，說了半天，看着公式也暈，下面以一個實際的信號來做說明。
假設我們有一個信號，它含有2V的直流分量，頻率為50Hz、相位為-30度、幅度為3V的交流信號，以及一個頻率為75Hz、相位為90度、幅度為1.5V的交流信號。用數學表達式就是如下：
S=2+3cos(2pi50t-pi30/180)+1.5cos(2pi75t+pi90/180)
式中cos參數為弧度，所以-30度和90度要分別換算成弧度。我們以256Hz的采樣率對這個信號進行采樣，總共采樣256點。按照我們上面的分析，Fn=(n-1)*Fs/N，我們可以知道，每兩個點之間的間距就是1Hz，第n個點的頻率就是n-1。我們的信號有3個頻率：0Hz、50Hz、75Hz，應該分別在第1個點、第51個點、第76個點上出現峰值，其它各點應該接近0。實際情況如何呢？我們來看看FFT的結果的模值如圖所示：
從圖中我們可以看到，在第1點、第51點、和第76點附近有比較大的值。我們分別將這三個點附近的數據拿上來細看：
1點： 512+0i
2點： -2.6195E-14 - 1.4162E-13i
3點： -2.8586E-14 - 1.1898E-13i
50點：-6.2076E-13 - 2.1713E-12i
51點：332.55 - 192i
52點：-1.6707E-12 - 1.5241E-12i
75點：-2.2199E-13 -1.0076E-12i
76點：3.4315E-12 + 192i
77點：-3.0263E-14 +7.5609E-13i
很明顯，1點、51點、76點的值都比較大，它附近的點值都很小，可以認為是0，即在那些頻率點上的信號幅度為0。接着，我們來計算各點的幅度值。分別計算這三個點的模值，結果如下：
1點： 512
51點：384
76點：192
按照公式，可以計算出直流分量為：512/N=512/256=2；50Hz信號的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號的幅度為192/(N/2)=192/(256/2)=1.5。可見，從頻譜分析出來的幅度是正確的。
然后再來計算相位信息。直流信號沒有相位可言，不用管它。先計算50Hz信號的相位，atan2(-192, 332.55)=-0.5236,結果是弧度，換算為角度就是180(-0.5236)/pi=-30.0001。再計算75Hz信號的相位，atan2(192, 3.4315E-12)=1.5708弧度，換算成角度就是1801.5708/pi=90.0002。可見，相位也是對的。
根據FFT結果以及上面的分析計算，我們就可以寫出信號的表達式了，它就是我們開始提供的信號。
總結
假設采樣頻率為Fs，采樣點數為N，做FFT之后，某一點n（n從1開始）表示的頻率為：Fn=(n-1)*Fs/N；該點的模值除以N/2就是對應該頻率下的信號的幅度（對於直流信號是除以N）；該點的相位即是對應該頻率下的信號的相位。相位的計算可用函數atan2(b,a)計算。atan2(b,a)是求坐標為(a,b)點的角度值，范圍從-pi到pi。要精確到xHz，則需要采樣長度為1/x秒的信號，並做FFT。要提高頻率分辨率，就需要增加采樣點數，這在一些實際的應用中是不現實的，需要在較短的時間內完成分析。解決這個問題的方法有頻率細分法，比較簡單的方法是采樣比較短時間的信號，然后在后面補充一定數量的0，使其長度達到需要的點數，再做FFT，這在一定程度上能夠提高頻率分辨力。具體的頻率細分法可參考相關文獻。
% [附錄：本測試數據使用的matlab程序]
close all; % 先關閉所有圖片
Adc=2; % 直流分量幅度
A1=3; % 頻率F1信號的幅度
A2=1.5; % 頻率F2信號的幅度
F1=50; % 信號1頻率(Hz)
F2=75; % 信號2頻率(Hz)
Fs=256; % 采樣頻率(Hz)
P1=-30; % 信號1相位(度)
P2=90; % 信號相位(度)
N=256; % 采樣點數
t=[0:1/Fs:N/Fs]; % 采樣時刻
% 信號
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
% 顯示原始信號
plot(S);
title('原始信號');
figure;
Y = fft(S,N); % 做FFT變換
Ayy = (abs(Y)); % 取模
plot(Ayy(1:N)); % 顯示原始的FFT模值結果
title('FFT 模值');
figure;
Ayy=Ayy/(N/2); % 換算成實際的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; % 換算成實際的頻率值
plot(F(1:N/2),Ayy(1:N/2)); % 顯示換算后的FFT模值結果
title('幅度-頻率曲線圖');
figure;
Pyy=[1:N/2];
for i=1:N/2
Pyy(i)=phase(Y(i)); % 計算相位
Pyy(i)=Pyy(i)*180/pi; % 換算為角度
end;
plot(F(1:N/2),Pyy(1:N/2)); % 顯示相位圖
title('相位-頻率曲線圖');