兩種含義

學習信號處理經常會被各種變換搞得暈頭轉向，這也是很正常的事，大可不必驚慌。暈的原因有兩個，一是各種變換公式看上去差不多，中間有或多或少的聯系，但又很不一樣，可能適用于不同的情況，有不同的限制等；二是這門學科確實起名太混淆了點，比如DTFT（Discrete Time Fourier Transform，離散時間傅里葉變換），DFT（Discrete Fourier Transform，離散傅里葉變換），DFS（Discrete Fourier Series，離散傅里葉級數），FFT（Fast Fourier Transform，快速傅里葉變換）。這些變換名字差不多，公式也差不多，不太容易搞清楚各自確切的含義和互相的聯系。

本文不去很深入的研究每一個變換，只關注工程實踐中最常用的FFT，來簡單剖析一下它的含義，性質和一些用法。

FFT的含義可以從兩個角度去理解，首先是DFS角度，其次是DTFT角度。

從DFS角度來說：FFT本質上就DFT，而DFT和DFS沒有什么不同。

第一句很好理解，因為FFT就是DFT的快速算法。

快了多少呢？DFT的時間復雜度是N^2，FFT的時間復雜度是NlogN。所以，假設N=1024，那就快了N/logN=102.4倍，可以說天壤之別。

現在看DFS：

DFS是離散周期序列x[n]到離散周期序列X[k]的變換，x[n]和X[k]周期相同均為N。因為非周期信號工程中更常見，科學家們就截取了DFS對的各自一個周期，定義為DFT變換。這樣DFT本質上就是DFS，只是隱含了周期性的假設。

所以，對一個非周期信號x[n]進行長度為N的DFT變換，首先是對x[n]進行周期為N的周期延拓，再求這個周期信號的DFS變換X(k)，最后截取X(k)一個周期，就得到x[n]的DFT。

因為離散序列通過DFT（或DFS）變換后還是離散序列，這意味著變換前后都很方便用計算機處理，所以DFT在實踐中非常重要。

另外還可以從DTFT的角度來理解FFT。

先看DTFT的公式：

DTFT是離散非周期信號到連續周期信號的變換，頻域是周期連續信號，顯然是以2Pi為周期（所有頻率都是以2Pi為周期）。DFS是對這個非周期信號進行周期延拓，在頻域就是對DTFT進行采樣。我們可以認為DFS設定的信號的周期延拓的周期N（也就是FFT的長度N）為DTFT一個周期內的采樣個數，N設定的越長，采樣越頻繁。

比如下圖：

圖1 一個 [1 1 1 1] 的有限長序列分別進行長度為8，16，128，1024的FFT

可見DFS選擇不同的周期N延拓非周期信號的時候，相當于對DTFT的每個周期進行采樣點為N的采樣。因為DTFT以[0, 2Pi)為周期，則DFS中0對應頻率0，N對應頻率2Pi。

注意上面都是說的DFS，而FFT本質上就是DFS。

N的選擇

計算FFT的時候，最好選擇長度N為2的n次冪，即N=2^n。因為FFT是采用“分治”（divide-and-conquer）思想實現的算法，層層二分顯然是最好的分治，這樣長度就需要是2的N次冪。現代的FFT算法也支持長度N為任意值，都可以得到正確的結果，只是速度上的快慢差別而已。最慢的情況是N是一個質數，這樣算法無法進行任何的分解，好一點的情況是N是一個非質數，最好的情況是N是2的n次冪。比如我們用 Matlab（R2016a）做一個實驗：

代碼：

A=[1,2,3];

tic    

fft(A, 121000);

toc;

tic

fft(A, 121001);

toc;

tic

fft(A, 131072);

toc;

結果：

Elapsed time is 0.006575 seconds.

Elapsed time is 0.039950 seconds.

Elapsed time is 0.004747 seconds.

當N等于質數121001時，FFT的運行時間大概是121000時的6倍，而N取值131072（2^17）時算得最快。

一些性質

共軛對稱

從FFT的公式可輕松推出，FFT變換后，第k個頻點和第N-k個頻點是共軛關系。這樣，第k個頻點的信息量完全等于第N-k個頻點的信息量，因此在頻域我們可以只處理差不多一半的頻點，處理完畢后再利用共軛的性質把另一半恢復出來，減少了一半的計算量。但這個“差不多一半”不是精確的“N/2”，而是N/2+1。因為頻點的取值范圍是[0, N-1]，0頻點沒有頻點和它共軛，但必須把它算在信息里面。

相位信息

信號的相位信息，在經過FFT之后，通過 arctan(X) 反映出來，比如求這個信號的相位信息：

smallvalue = 1e-10;

fs = 1000;

N = 1000;

t = (0:N-1)/fs;

f=linspace(0,fs,N+1);

f1= 50;

y = cos(2*pi*f1*t+pi/4);

Y = fft(y);

Y(intersect(find(real(Y)< smallvalue), find(imag(Y)< smallvalue)))=0;

plot(f(1:N),angle(Y));

圖二 信號的相位信息

注意，如果FFT計算出的復數等于0的話，因為 Matlab 的計算精度問題，其實都是一些很小的值，比如實部 re = 1.05e-12，虛部 im = 3e-12。這些值對功率譜沒什么影響，但取 arctan(im/re) 時卻是一些很大的值，范圍在[-Pi, Pi]，這顯然是不合理的。

所以我手動去掉了某些極小值。

Y(intersect(find(real(Y)

加窗

在實際應用中，對一段信號做FFT通常要先加窗，根本原因也是FFT隱藏的周期性。

比如分別對下圖第一列中的兩個信號求FFT，首先是做周期延拓，再求DFS。第一行的信號，周期延拓后正好是一個完美的連續函數，DFT完全不受影響，仍然只有原始信號的頻率，其實就是一個單一頻率。而第二行的信號，周期延拓后變成了一個不連續的函數，簡單的說，也就是引入了其他頻率分量，因此DFT后產生了原始信號沒有的頻率，這就叫頻譜泄漏。

而加窗使得信號周期延拓后變得連續，減少了頻譜泄漏。

圖3 頻譜泄漏的產生

注意，如果正好是周期信號一個周期的延拓，那樣是沒有頻譜泄漏的，理論上也不需要加窗。但實際工程應用中的信號過于復雜，不可能做到這一點，所以加窗就變成了一個標準操作。

加窗導致信號的能量發生了變化，想當于每個點乘以 sum_of_window / len_of_window，那么IFFT之后需要乘以 len_of_window / sum_of_window 把能量補償回來。

線性卷積

一個線性時不變系統的輸出就是輸入信號和系統沖激響應的線性卷積。

線性卷積的公式是：

可見，加入x1的長度是N1，x2的長度是N2，那么x3的長度就是N1+N2-1。

但是因為FFT隱含的周期性，FFT的時域卷積是周期信號的線性卷積。而周期信號一旦卷積起來，就好像是按照周期N循環一樣，所以稱為循環卷積。

線性卷積的快速解法：

1. 選取 N>=(N1+N2-1)；
2. 計算兩個序列 x1[n] 和 x2[n] 的N點FFT X1[k] 和 X2[k]；
3. 計算乘積 X3[k]=X1[k]*X2[k]；
4. 計算 X3[k] 的IFFT逆變換得到 x1[n] 和 x2[n] 的循環卷積，由于N足夠大，其實等于線性卷積。

但還有一種情況是實時的線性卷積，這種情況很常見，比如實時錄制的音頻流卷積一個濾波器，再實時的輸出處理后的信號。我們也可以利用FFT來節省計算資源，但顯然音頻流的長度是不可知的，也就無法通過上面的方式進行計算。

但可以通過如下途徑來“分段”循環卷積。

簡單說明一下：h是濾波器系數，假如它的長度是L，則定義FFT的長度N=2L，在h的后面補零。in_data是輸入數據，out_data是輸出數據。長度為N的輸入數據需要兩步才能處理完成，每一步都需要做長度為N的FFT/IFFT，之后，丟棄生成的前N/2長度的數據，保留后面N/2長度的數據到輸出緩沖區中。所以若產生N長度的數據，需要做兩次N長度的FFT/IFFT。具體過程請參考圖示進行推導。

圖4 分段循環卷積