兩種含義

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

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

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

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

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

快了多少呢？DFT的時間復(fù)雜度是N^2，F(xiàn)FT的時間復(fù)雜度是NlogN。所以，假設(shè)N=1024，那就快了N/logN=102.4倍，可以說天壤之別。

現(xiàn)在看DFS：

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

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

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

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

先看DTFT的公式：

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

比如下圖：

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

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

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

N的選擇

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

代碼：

A=[1,2,3];

tic    

fft(A, 121000);

toc;

tic

fft(A, 121001);

toc;

tic

fft(A, 131072);

toc;

結(jié)果：

Elapsed time is 0.006575 seconds.

Elapsed time is 0.039950 seconds.

Elapsed time is 0.004747 seconds.

當(dāng)N等于質(zhì)數(shù)121001時，F(xiàn)FT的運行時間大概是121000時的6倍，而N取值131072（2^17）時算得最快。

一些性質(zhì)

共軛對稱

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

相位信息

信號的相位信息，在經(jīng)過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計算出的復(fù)數(shù)等于0的話，因為 Matlab 的計算精度問題，其實都是一些很小的值，比如實部 re = 1.05e-12，虛部 im = 3e-12。這些值對功率譜沒什么影響，但取 arctan(im/re) 時卻是一些很大的值，范圍在[-Pi, Pi]，這顯然是不合理的。

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

Y(intersect(find(real(Y)

加窗

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

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

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

圖3 頻譜泄漏的產(chǎn)生

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

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

線性卷積

一個線性時不變系統(tǒng)的輸出就是輸入信號和系統(tǒng)沖激響應(yīng)的線性卷積。

線性卷積的公式是：

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

但是因為FFT隱含的周期性，F(xiàn)FT的時域卷積是周期信號的線性卷積。而周期信號一旦卷積起來，就好像是按照周期N循環(huán)一樣，所以稱為循環(huá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] 的循環(huán)卷積，由于N足夠大，其實等于線性卷積。

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

但可以通過如下途徑來“分段”循環(huán)卷積。

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

圖4 分段循環(huán)卷積