數字信号處理3: 快速傅裡葉變換（FFT）（含代碼）

文章目錄

• 1. FFT推導
• 2. FFT為什麼快?
• 3. 一些加速措施
• • 3.1 查表法計算三角函數
  • 3.2 奇偶分解
• 4. FFT代碼

在之前的文章《傅裡葉變換》中，我們已經推導了連續傅裡葉變換和離散傅裡葉變換。由于計算機的發展，離散傅裡葉變換（DFT）可謂是信号處理的殺手锏。但是離散傅裡葉變換計算量巨大，通常在實時信号處理時是無法使用的，直到快速傅裡葉變換（FFT）算法被發現。

與DFT不同，FFT是一種算法而非理論，是以它無法隻靠一兩個公式就描述出來。接下來我們從DFT出發進行FFT的推導。

當然，大多數人可能來看這個部落格就隻是為了尋找可用的FFT代碼，但是肯定還有一些同學對這個過程會感到好奇，是以我還是先講推導過程，代碼在最後一部分。

1. FFT推導

首先大家應該還記得DFT的公式

X [ k w 0 ] = ∑ n = < N > x [ n ] e − j k w 0 n X[kw_0]=\sum_{n=<N>}x[n]e^{-jkw_0n} X[kw0​]=n=<N>∑​x[n]e−jkw0​n

x [ n ] = 1 N ∑ k = < N > X [ k w 0 ] e j k w 0 n x[n]=\frac{1}{N}\sum_{k=<N>}X[kw_0]e^{jkw_0n} x[n]=N1​k=<N>∑​X[kw0​]ejkw0​n

其中， N N N為DFT長度，也就是多少個資料點， w 0 = 2 π / N w_0=2\pi/N w0​=2π/N為基波頻率。首先從第一個公式開始。

為了友善表示，用 X [ k ] X[k] X[k]來表示這是DFT的第幾個點（畢竟是離散的），對于DFT可以寫成更具體的形式：

X [ k ] = ∑ n = 0 N − 1 x [ n ] e − j k w 0 n X[k]=\sum_{n=0}^{N-1}x[n]e^{-jkw_0n} X[k]=n=0∑N−1​x[n]e−jkw0​n

這種形式是更符合計算機表示的。

将上式按照 n n n的奇偶性分成兩部分：

X [ k ] = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k w 0 ( 2 n ) + ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k w 0 ( 2 n + 1 )   = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k w 0 ( 2 n ) + e − j k w 0 ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k w 0 ( 2 n )   = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k 2 π N / 2 n + e − j k w 0 ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k 2 π N / 2 n X[k]=\sum_{n=0}^{N/2-1} x[2n]e^{-jkw_0(2n)}+\sum_{n=0}^{N/2-1}x[2n+1]e^{-jkw_0(2n+1)}\\ \space\\ =\sum_{n=0}^{N/2-1} x[2n]e^{-jkw_0(2n)}+e^{-jkw_0}\sum_{n=0}^{N/2-1}x[2n+1]e^{-jkw_0(2n)}\\ \space\\ =\sum_{n=0}^{N/2-1} x[2n]e^{-jk\frac{2\pi}{N/2}n}+e^{-jkw_0}\sum_{n=0}^{N/2-1}x[2n+1]e^{-jk\frac{2\pi}{N/2}n} X[k]=n=0∑N/2−1​x[2n]e−jkw0​(2n)+n=0∑N/2−1​x[2n+1]e−jkw0​(2n+1) =n=0∑N/2−1​x[2n]e−jkw0​(2n)+e−jkw0​n=0∑N/2−1​x[2n+1]e−jkw0​(2n) =n=0∑N/2−1​x[2n]e−jkN/22π​n+e−jkw0​n=0∑N/2−1​x[2n+1]e−jkN/22π​n

注意 k k k的取值範圍仍然是 [ 0 , N − 1 ] [0, N-1] [0,N−1]，但是由于

e − j ( k + N ) 2 π N n = e − j k 2 π N n ⋅ e − j 2 n π = e − j k 2 π N n e^{-j(k + N)\frac{2\pi}{N}n}=e^{-jk\frac{2\pi}{N}n}\cdot e^{-j2n\pi}=e^{-jk\frac{2\pi}{N}n} e−j(k+N)N2π​n=e−jkN2π​n⋅e−j2nπ=e−jkN2π​n

同理

e − j ( k + N / 2 ) 2 π N / 2 n = e − j k 2 π N / 2 n ⋅ e − j 2 n π = e − j k 2 π N / 2 n e^{-j(k + N/2)\frac{2\pi}{N/2}n}=e^{-jk\frac{2\pi}{N/2}n}\cdot e^{-j2n\pi}=e^{-jk\frac{2\pi}{N/2}n} e−j(k+N/2)N/22π​n=e−jkN/22π​n⋅e−j2nπ=e−jkN/22π​n

是以不妨令 k 1 k_1 k1​的取值範圍為 [ 0 , N / 2 − 1 ] [0, N/2-1] [0,N/2−1]。令

X 1 [ k 1 ] = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k 1 2 π N / 2 n X_1[k_1]=\sum_{n=0}^{N/2-1} x[2n]e^{-jk_1\frac{2\pi}{N/2}n} X1​[k1​]=n=0∑N/2−1​x[2n]e−jk1​N/22π​n

X 2 [ k 1 ] = ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k 1 2 π N / 2 n X_2[k_1]=\sum_{n=0}^{N/2-1}x[2n+1]e^{-jk_1\frac{2\pi}{N/2}n} X2​[k1​]=n=0∑N/2−1​x[2n+1]e−jk1​N/22π​n

則

X [ k ] = X 1 [ k 1 ] + e − j k w 0 X 2 [ k 1 ] X[k]=X_1[k_1]+e^{-jkw_0}X_2[k_1] X[k]=X1​[k1​]+e−jkw0​X2​[k1​]

其中， k 1 = k % ( N / 2 ) k_1=k \% (N/2) k1​=k%(N/2)，百分号是取餘數。

很顯然， X 1 [ k 1 ] X_1[k_1] X1​[k1​]就是 x [ n ] x[n] x[n]的偶數項的DFT，而 X 2 [ k 1 ] X_2[k_1] X2​[k1​]就是 x [ n ] x[n] x[n]的奇數項的DFT，隻是基波頻率有所變化（畢竟樣點數變成了原來的一半）。

如果将上式繼續按照這種方式分解。

X 1 [ k 1 ] = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k 1 2 π N / 2 n   = ∑ n = 0 N / 4 − 1 x [ 4 n ] e − j k 1 2 π N / 2 ( 2 n ) + ∑ n = 0 N / 4 − 1 x [ 4 n + 2 ] e − j k 1 2 π N / 2 ( 2 n + 1 )   = ∑ n = 0 N / 4 − 1 x [ 4 n ] e − j k 1 2 π N / 4 n + e − j k 1 2 π N / 2 ∑ n = 0 N / 4 − 1 x [ 4 n + 2 ] e − j k 1 2 π N / 4 n X_1[k_1]=\sum_{n=0}^{N/2-1} x[2n]e^{-jk_1\frac{2\pi}{N/2}n}\\ \space\\ =\sum_{n=0}^{N/4-1}x[4n]e^{-jk_1\frac{2\pi}{N/2}(2n)}+\sum_{n=0}^{N/4-1}x[4n+2]e^{-jk_1\frac{2\pi}{N/2}(2n+1)}\\ \space\\ =\sum_{n=0}^{N/4-1}x[4n]e^{-jk_1\frac{2\pi}{N/4}n}+e^{-jk_1\frac{2\pi}{N/2}}\sum_{n=0}^{N/4-1}x[4n+2]e^{-jk_1\frac{2\pi}{N/4}n} X1​[k1​]=n=0∑N/2−1​x[2n]e−jk1​N/22π​n =n=0∑N/4−1​x[4n]e−jk1​N/22π​(2n)+n=0∑N/4−1​x[4n+2]e−jk1​N/22π​(2n+1) =n=0∑N/4−1​x[4n]e−jk1​N/42π​n+e−jk1​N/22π​n=0∑N/4−1​x[4n+2]e−jk1​N/42π​n

同樣，可以令 k 2 = k 1 % ( N / 4 ) k_2=k_1\%(N/4) k2​=k1​%(N/4)，即 k 2 k_2 k2​的取值範圍為 [ 0 , N / 4 − 1 ] [0, N/4-1] [0,N/4−1]。于是可得：

X 1 [ k 1 ] = X 3 [ k 2 ] + e − j k 1 w 0 2 X 4 [ k 2 ] X_1[k_1]=X_3[k_2]+e^{-jk_1w_02}X_4[k_2] X1​[k1​]=X3​[k2​]+e−jk1​w0​2X4​[k2​]

這樣就很顯然了，整個公式就是把長的DFT轉化為短的DFT，每次把一個長度為N的DFT按照序号的奇偶分為兩個長度為 N / 2 N/2 N/2的DFT。一直到最後分成一個數的DFT，一個數的DFT就是這個數本身。

X [ 0 ] = ∑ n = 0 0 x [ n ] e − j 0 = x [ 0 ] X[0]=\sum_{n=0}^{0}x[n]e^{-j0}=x[0] X[0]=n=0∑0​x[n]e−j0=x[0]

當然上面是以計算系統序号從0開始的情況來說的（絕大多數計算機語言的序号都是從0開始的），其實如果計算系統序号從1開始也是一樣的：

X [ 1 ] = ∑ n = 1 1 x [ n ] e − j w 0 = x [ 1 ] X[1]=\sum_{n=1}^{1}x[n]e^{-jw_0}=x[1] X[1]=n=1∑1​x[n]e−jw0​=x[1]

因為是一個數的DFT，是以 N = 1 N=1 N=1， w 0 = 2 π / N = 2 π w_0=2\pi/N=2\pi w0​=2π/N=2π。由歐拉公式 e − j w 0 = 1 e^{-jw_0}=1 e−jw0​=1。

稍微有些跑題，接着講。我們計算出來了一個數的DFT之後，根據上面分解的那些公式顯然可以用它們來組合出兩個數的DFT，進而得出四個數的DFT，一直到 N N N。其實這整個過程也隐含了一個很重要的條件，就是FFT的長度必須為2的幂次方。

好了，現在我們知道，長的DFT可以通過短的DFT組合計算。但是還有兩個頭疼的問題：

• 在我們分解的過程中， x [ n ] x[n] x[n]的順序已經被打亂了，那麼這個順序要怎麼排？
• 每次分解後的奇數項求和前面有個系數，這個系數如何确定？

接下來用一個具體的例子。

假設現在要求一個8個點的DFT。

表中綠色表示每次分解後的偶數項，黃色表示每次分解的奇數項。分解3次後到底。可以預見，對于長度為 N N N的序列(N必須為2的幂次方)，分解需要進行 l o g 2 N log_2N log2​N步。

乍一看似乎看不出有什麼規律，其實這種序号變化可以使用一種位反轉的方法來計算。如下表:

需要注意，位反轉并非二進制取反，而是“低位與高位交換”，并且反轉點需要根據具體的長度來定。比如如果 N = 16 N=16 N=16時，原序号為1，反轉後就是8。

然後看那個系數問題。

根據之前的分解過程，第一次分解後，系數是

e − j k w 0 e^{-jkw_0} e−jkw0​

第二次是

e − j k 1 w 0 2 e^{-jk_1w_02} e−jk1​w0​2

第三次是

e − j k 2 w 0 4 e^{-jk_2w_04} e−jk2​w0​4

對于取值範圍， k n k_{n} kn​總是 k n − 1 k_{n-1} kn−1​的一半。

第n次分解，就是

e − j k n w 0 2 n − 1 e^{-jk_nw_02^{n-1}} e−jkn​w0​2n−1

OK，以上基本就是所有FFT的推導了，結合上述過程，我們就可以得到FFT計算的蝶形圖：

圖中 W M k = e − j k 2 π M W_M^k=e^{-jk\frac{2\pi}{M}} WMk​=e−jkM2π​， k k k的取值範圍為 [ 0 , M − 1 ] 。 [0, M-1]。 [0,M−1]。這裡的 M M M并非原序列 x [ n ] x[n] x[n]的長度，而是每次被分解的那個序列的長度。比如在第二次分解時，我們得到 e − j k 1 w 0 2 = e − j k 1 2 π N / 2 e^{-jk_1w_02}=e^{-jk_1\frac{2\pi}{N/2}} e−jk1​w0​2=e−jk1​N/22π​， M = N / 2 = 4 M=N/2=4 M=N/2=4。

2. FFT為什麼快?

在未使用任何程式上的加速技巧時：

對于一個長度為N的序列，使用傳統DFT方法，我們需要計算 N 2 N^2 N2次複數乘法， N 2 N^2 N2次求和， N 2 N^2 N2次複數幂的計算（當然你也可以使用歐拉公式将其轉化為三角函數的計算，但這仍然是不小的開支）。

而使用FFT時，我們需要計算 N + N / 2 + N / 4 + ⋯ + 4 N + N/2 + N/4 + \cdots + 4 N+N/2+N/4+⋯+4（等比數列）次複數幂計算，總共有 l o g 2 N − 1 log_2N-1 log2​N−1項，也就是 2 N ( 1 − ( 1 2 ) l o g 2 N − 1 ) 2N(1-(\frac{1}{2})^{log_2N-1}) 2N(1−(21​)log2​N−1)。複數乘法以及加法都要進行 N l o g 2 N Nlog_2N Nlog2​N次。

3. 一些加速措施

3.1 查表法計算三角函數

對于計算機來說，無論是計算e的複數幂，或者求三角函數，都是非常耗時的。在arm加速庫中，對于三角函數是采用查表法，預先生成一個周期的sin值儲存在數組裡，需要精度高就儲存數組長一些，使用的時候根據角度在數組中查找對應值。當然，如果你不嫌煩的話，僅生成1/4個周期的sin值即可，不過查找起來比較麻煩。

3.2 奇偶分解

首先看歐拉公式，對于DFT公式中e的幂次方部分，：

e − j k 2 π N n = c o s ( k 2 π N n ) − j s i n ( k 2 π N n ) e^{-jk\frac{2\pi}{N}n}=cos(k\frac{2\pi}{N}n)-jsin(k\frac{2\pi}{N}n) e−jkN2π​n=cos(kN2π​n)−jsin(kN2π​n)

我們首先看一下不同頻率的複信号。由于是隻看頻率，是以可剔除 n n n，隻看 k k k。（這裡不了解的同學可以回顧一下上一篇講DFT的文章中通俗了解的那部分）。

e − j k 2 π N = c o s ( k 2 π N ) − j s i n ( k 2 π N ) e^{-jk\frac{2\pi}{N}}=cos(k\frac{2\pi}{N})-jsin(k\frac{2\pi}{N}) e−jkN2π​=cos(kN2π​)−jsin(kN2π​)

從這兩個三角函數中可以看到， c o s ( k 2 π N ) cos(k\frac{2\pi}{N}) cos(kN2π​)是關于 k = a N 2 k=a\frac{N}{2} k=a2N​偶對稱的，而 s i n ( k 2 π N ) sin(k\frac{2\pi}{N}) sin(kN2π​)是關于 k = a N 2 k=a\frac{N}{2} k=a2N​奇對稱的。其中a是整數。

是以我們可以得到一個重要結論：

對于純實信号，DFT後 X [ k ] X[k] X[k]的實部是關于 N 2 \frac{N}{2} 2N​偶對稱的， X [ k ] X[k] X[k]的虛部是關于 N 2 \frac{N}{2} 2N​奇對稱的。

對于純虛信号，DFT後 X [ k ] X[k] X[k]的實部是關于 N 2 \frac{N}{2} 2N​奇對稱的， X [ k ] X[k] X[k]的虛部是關于 N 2 \frac{N}{2} 2N​偶對稱的。

由于我們接觸的基本都是實信号，如此一來如果按照傳統FFT，就會有相當多的計算是浪費的。是以，我們不妨把純實序列 x [ n ] x[n] x[n]當作複數序列，它的偶數項是實部，奇數項是虛部。這樣一來我們僅需要做計算 N / 2 N/2 N/2長度的FFT即可。

待上一步計算完成後，我們需要把實部和虛部的FFT結果分離出來。這個時候就要用到之前對稱性的結論，可以知道 X [ k ] X[k] X[k]的實部和虛部都是由一個偶對稱序列和一個奇對稱序列疊加而成。如何把偶對稱序列和奇對稱序列分離出來？這就需要使用奇偶分解：

x E [ n ] = x [ n ] + x [ N − n ] 2   x O [ n ] = x [ n ] − x [ N − n ] 2 x_E[n]=\frac{x[n] + x[N - n]}{2}\\ \space\\ x_O[n]=\frac{x[n] - x[N - n]}{2} xE​[n]=2x[n]+x[N−n]​ xO​[n]=2x[n]−x[N−n]​

N為 x [ n ] x[n] x[n]的長度。同時我們令 x [ N ] = x [ 0 ] x[N]=x[0] x[N]=x[0]。這對于 X [ k ] X[k] X[k]也是适用的，因為如果不限制k的取值範圍， X [ k ] X[k] X[k]是周期為N的周期函數。這一點可以從DFT公式得出。

接下來我們就可以分解出偶對稱序列和奇對稱序列。 X [ k ] X[k] X[k]由實部序列和虛部序列組合而成：

X [ k ] = X r [ k ] + X i [ k ] X[k]=X_r[k]+X_i[k] X[k]=Xr​[k]+Xi​[k]

由于我們把長度為 N N N的實數FFT轉化為長度為 N / 2 N/2 N/2的複數FFT，是以上式中k的範圍是 [ 0 , N / 2 − 1 ] [0, N/2-1] [0,N/2−1]。 X r [ k ] X_r[k] Xr​[k]是實部序列， X i [ k ] X_i[k] Xi​[k]是虛部序列。

設 X e X_{e} Xe​是 x [ n ] x[n] x[n]的偶數項的FFT結果， X o X_{o} Xo​是奇數項的FFT結果（注意此時奇數項被我們視作純虛數序列）。 X e r X_{er} Xer​是 X e X_{e} Xe​的實部序列， X e i X_{ei} Xei​是 X e X_{e} Xe​的虛部序列。 X o r X_{or} Xor​是 X o X_{o} Xo​的實部序列， X o i X_{oi} Xoi​是 X o X_{o} Xo​的虛部序列。

X e [ k ] = X e r [ k ] + X e i [ k ] X o [ k ] = X o r [ k ] + X o i [ k ] X_{e}[k] = X_{er}[k] + X_{ei}[k]\\ X_{o}[k] = X_{or}[k] + X_{oi}[k] Xe​[k]=Xer​[k]+Xei​[k]Xo​[k]=Xor​[k]+Xoi​[k]

那麼

X r [ k ] = X e r [ k ] + X o r [ k ] X_r[k]=X_{er}[k] + X_{or}[k] Xr​[k]=Xer​[k]+Xor​[k]

由前面的結論， X e r X_{er} Xer​是偶對稱的， X o r X_{or} Xor​是奇對稱的。故

X e r [ k ] = X r [ k ] + X r [ N 2 − k ] 2   X o r [ k ] = X r [ k ] − X r [ N 2 − k ] 2 X_{er}[k]=\frac{X_r[k] + X_r[\frac{N}{2}-k]}{2}\\ \space\\ X_{or}[k]=\frac{X_r[k] - X_r[\frac{N}{2}-k]}{2} Xer​[k]=2Xr​[k]+Xr​[2N​−k]​ Xor​[k]=2Xr​[k]−Xr​[2N​−k]​

同理

X e i [ k ] = X i [ k ] − X i [ N 2 − k ] 2   X o i [ k ] = X i [ k ] + X i [ N 2 − k ] 2 X_{ei}[k]=\frac{X_i[k] - X_i[\frac{N}{2}-k]}{2}\\ \space\\ X_{oi}[k]=\frac{X_i[k] + X_i[\frac{N}{2}-k]}{2} Xei​[k]=2Xi​[k]−Xi​[2N​−k]​ Xoi​[k]=2Xi​[k]+Xi​[2N​−k]​

這樣一來，我們就得到了 X e [ k ] X_{e}[k] Xe​[k]和 X o [ k ] X_{o}[k] Xo​[k]。回顧我們将實數FFT轉化為複數DFT的那步，我們将偶數項和奇數項分開做FFT，這和蝶形運算非常相似，不過在蝶形運算中，我們仍然将奇數項看作實數，是以，我們将 X o [ k ] X_o[k] Xo​[k]每一項乘 − j -j −j，即可得到奇數項的FFT結果，此時它是正常的實數FFT結果。接下來，我們就可以進行蝶形運算的最後一步合并，得到的就是最終的FFT結果。

4. FFT代碼

代碼請通路我的github：FFT

在代碼中，我使用了查表法。對于純實數的FFT，我将其轉換為了複數FFT。在代碼中并沒有使用複數，而是使用float數組代替複數數組，格式為{real, imag, real ,imag …}。