快速傅立葉變換是離散傅立葉變換的一種快速實現方式，快速傅立葉變換可用於多項式乘法、大數乘法、卷積等操作，把原本的O(n^2)計算量優化到了O(nlogn)，這是質的飛躍。我們現在能這麼快的網上衝浪，這個演算法居功至偉，讓我們為它鼓掌！
O(n^2)和O(nlogn)的差距在哪裡，這裡的log底數是2。
如果要處理一萬個資料，O(n^2)需要計算一億次，O(nlogn)則僅需要計算十三萬次。
接下來，下面我們會看到快速傅立葉變換是如何一步步地將O(n^2)優化到O(nlogn)的。
暴力破解法
根據上一篇文章所述，離散傅立葉公式如下:
舉個例子，當N為4的時候，離散傅立葉變換的運算相當於如下矩陣乘積:
按照這個操作，我們也可以實現N*N操作的離散傅立葉變換程式碼。
#include<cmath>
#include<iostream>
#include<complex>
#include<vector>
using namespace std;
const double pi = acos(-1.0);
vector<complex<double>> fft(vector<complex<double>> input)
{
vector<complex<double>> output;
int N = input.size(), k = 0, n = 0;
for (k = 0; k < N; k++)
{
complex<double> sum(0,0);
for (n = 0; n < N; n++)
{
complex<double> wn(polar<double>(1, -2*pi/N * n * k));
sum += input[n] * wn;
}
output.push_back(sum);
}
return output;
}
int main()
{
vector<complex<double>> in{1, 2, 3};
auto re = fft(in);
for (auto i: re)
{
cout << "real:" << i.real() << " imag:" << i.imag() << endl;
}
return 0;
}
離散傅立葉變換運算有兩層for迴圈，時間複雜度是O(n^2)。
從上面那張矩陣圖可以看到，公式的
這個$$ e^{-i2\pi } $$ 是固定不變的，變化的只有N, n, k。
再此多說一句:$$ e^{-i2\pi } $$是順時針旋轉一週的意思。同理如下圖的$$ e^{i\pi } $$ 也就是逆時針旋轉半周的意思。
至於旋轉就是這麼轉，這麼地轉動:
我們可以將$$ e^{-i2\pi } $$簡化表示成$$ \omega $$，並將N作為下標，且nk作為上標如：$$ \omega_{N}^{nk} $$
如果$$ \omega $$代表順時針旋轉一週，那麼$$ \omega_{N}^{nk} $$就代表，其步進值是順時針旋轉N分之一週，且步進次數是n*k次後的結果。
每一次旋轉的步進值就是單位複數根：
還是當N為4的時候舉個例子，運算矩陣如下：
命運輪迴，上述的ω再怎麼變，也跳不出下圖這四個紅點：
而當N為2時，ω再怎麼變，也跳不出下圖這兩個紅點：
旋轉定理
此時我們可以得到這樣的定理。
定理一
這是旋轉的週期性決定，旋轉一週就等於沒有旋轉。
定理二
當N為偶數的時候，圓的每個點都即為X軸對稱，也為Y軸對稱。
定理三
以上，當N、k均為偶數時，當N為之前的一半的時候，單位複數根也增大了一倍，所以是成立的。
根據定理一二三，我們可以把矩陣從這樣：
轉換成這樣:
看這裡的時候不必過分追究細節，具體推導在下面。
我們好像看到了一些東西，似乎我們只需要求得 $$ \omega _{2}^{1} $$ 和 $$ \omega _{4}^{1} $$ 即可？
進一步可以轉換成這樣:
嗯？好像看到了一點規律的感覺，當我們轉換一個N=4的序列的時候，我們似乎可以將它分成兩個N=2的序列，各自進行轉換。但是為什麼會這樣？
推導過程如下（也可以直接拉到最後看結論），以下參考羅博士的推導，我重新寫了公式。考慮一個4點DFT，其原始公式是：
奇偶下標項拆分一下:
提取公共項：
統一一下求和符號的下標n:
注意到N=4為偶數，根據單位復根的定理三，可以做一個比例變換，得到：
再對求和符號的下標做一個變換:
拆分成0,1和2,3:
其中k=2,3項可以轉換成以下公式:
根據定理一和定理二可以轉換成以下公式：
最終得到前後半部分的公式:
可以看到，如果把N=4的序列分成，偶數下標項和奇數下標項兩部分，分別做FFT得到A和B，那麼N=4的序列的傅立葉變換的結果：
這個時候再看看之前的這個矩陣，好像有點意思:
這裡為了加深理解，加上X0和X2，當N=2時的FFT如下：
遞迴法（自頂向下）
如果是這樣的話，原本之前我們的時間複雜度是O(N^2)，現在變成了O(N*N/2)。
進一步，如果N的長度是2次冪的話，那麼這個不斷拆分偶數和奇數下標項的過程可以一直遞迴下去。
例如，如果要計算一個N=8的序列，那麼拆分成兩個N=4的序列計算，再拆分成四個N=2的序列計算，再拆分八個N=2的序列計算，是不是有一點歸併排序的感覺了？
那麼最終計算量可以降低到O(N*logN)
void FFT_recursion(vector<complex<double>> &input)
{
int n = input.size();
if(n == 1)
return;
int mid= input.size()/2;
vector<complex<double>> A0, A1;
for(int i = 0; i < n; i += 2){//拆分奇偶下標項
A0.push_back(input[i]);
A1.push_back(input[i + 1]);
}
FFT_recursion(A0);
FFT_recursion(A1);
complex<double> w0(1,0);
complex<double> wn(polar<double>(1, -2*pi/n));//單位根
for(int i=0; i < mid; i++, w0*=wn){//合併多項式
input[i] = A0[i] + w0*A1[i];
input[i + mid] = A0[i] - w0*A1[i];
}
}
以上我們明白了優化的思路，上述程式碼雖然時間複雜度降到了O(nlogn)，但是空間複雜度也是O(nlogn)。下一篇將介紹時間複雜度O(nlogn)，同時空間複雜度為O(1)的操作……