兩種含義

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

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

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

從DFS角度來(lái)說(shuō)：FFT本質(zhì)上就DFT，而DFT和DFS沒(méi)有什么不同。

第一句很好理解，因?yàn)镕FT就是DFT的快速算法。

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

現(xiàn)在看DFS：

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

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

因?yàn)殡x散序列通過(guò)DFT（或DFS）變換后還是離散序列，這意味著變換前后都很方便用計(jì)算機(jī)處理，所以DFT在實(shí)踐中非常重要。

另外還可以從DTFT的角度來(lái)理解FFT。

先看DTFT的公式：

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

比如下圖：

圖1 一個(gè) [1 1 1 1] 的有限長(zhǎng)序列分別進(jìn)行長(zhǎng)度為8，16，128，1024的FFT

可見(jiàn)DFS選擇不同的周期N延拓非周期信號(hào)的時(shí)候，相當(dāng)于對(duì)DTFT的每個(gè)周期進(jìn)行采樣點(diǎn)為N的采樣。因?yàn)镈TFT以[0, 2Pi)為周期，則DFS中0對(duì)應(yīng)頻率0，N對(duì)應(yīng)頻率2Pi。

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

N的選擇

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

代碼：

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時(shí)，F(xiàn)FT的運(yùn)行時(shí)間大概是121000時(shí)的6倍，而N取值131072（2^17）時(shí)算得最快。

一些性質(zhì)

共軛對(duì)稱

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

相位信息

信號(hào)的相位信息，在經(jīng)過(guò)FFT之后，通過(guò) arctan(X) 反映出來(lái)，比如求這個(gè)信號(hào)的相位信息：

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));

圖二 信號(hào)的相位信息

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

所以我手動(dòng)去掉了某些極小值。

Y(intersect(find(real(Y)

加窗

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

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

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

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

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

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

線性卷積

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

線性卷積的公式是：

可見(jiàn)，加入x1的長(zhǎng)度是N1，x2的長(zhǎng)度是N2，那么x3的長(zhǎng)度就是N1+N2-1。

但是因?yàn)镕FT隱含的周期性，F(xiàn)FT的時(shí)域卷積是周期信號(hào)的線性卷積。而周期信號(hào)一旦卷積起來(lái)，就好像是按照周期N循環(huán)一樣，所以稱為循環(huán)卷積。

線性卷積的快速解法：

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

但還有一種情況是實(shí)時(shí)的線性卷積，這種情況很常見(jiàn)，比如實(shí)時(shí)錄制的音頻流卷積一個(gè)濾波器，再實(shí)時(shí)的輸出處理后的信號(hào)。我們也可以利用FFT來(lái)節(jié)省計(jì)算資源，但顯然音頻流的長(zhǎng)度是不可知的，也就無(wú)法通過(guò)上面的方式進(jìn)行計(jì)算。

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

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

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