一小時學會快速傅里葉變換（Fast Fourier Transform）

01-24

0x00 寫在前面

為了讓更多人能夠看到這個教程，希望大家收藏之前，也要點贊哦！！！蟹蟹大家的認可和鼓勵。??????

傅里葉變換

快速傅里葉變換（Fast Fourier Transform，FFT）是一種可在 $O(nlogn)$ 時間內完成的離散傅里葉變換（Discrete Fourier transform，DFT）演算法。

在演算法競賽中的運用主要是用來加速多項式的乘法。

考慮到兩個多項式 $A(x),B(x)$ 的乘積 $C(x)$ ，假設 $A(x)$ 的項數為 $n$ ，其係數構成的 $n$ 維向量為 $(a_0,a_1,a_2,...,a_{n-1})$ ， $B(x)$ 的項數為 $m$ ，其係數構成的 $m$ 維向量為 $(b_0,b_1,b_2,...,b_{m-1})$ 。

我們要求 $C(x)$ 的係數構成的 $n+m-1$ 維的向量，先考慮樸素做法。

可以用這段代碼表示：

for ( int i = 0 ; i < n ; ++ i ) for ( int j = 0 ; j < m ; ++ j ) { c [i + j] += a [i] * b [j] ; }

思路非常清晰，其時間複雜度是 $O(n^2)$ 的。

所以我們來學習快速傅里葉變換。

0x01 關於多項式

多項式有兩種表示方法，係數表達法與點值表達法

多項式的係數表示法

設多項式 $A(x)$ 為一個 $n-1$ 次的多項式，顯然，所有項的係數組成的係數向量 $(a_0,a_1,a_2,...,a_{n-1})$ 唯一確定了這個多項式。

$A(x)=sum_{i=0}^{n-1}a_icdot x^i$

多項式的點值表示法

將一組互不相同的 $(x_0,x_1,x_2,...,x_n)$ （叫插值節點）分別帶入 $A(x)$ ，得到 $n$ 個取值 $(y_0,y_1,y_2,...,y_n)$ .

其中

$y_i=sum_{j=0}^{n-1}a_jcdot x_i^j$

定理:

一個 $n-1$ 次多項式在 $n$ 個不同點的取值唯一確定了該多項式。

證明：

假設命題不成立，存在兩個不同的 $n-1$ 次多項式 $A(x),B(x)$ ，滿足對於任何 $i in [0, n - 1]$ ，有 $A(x_i) = B(x_i)$ 。
令 $C(x) = A(x) - B(x)$ ，則 $C(x)$ 也是一個 $n-1$ 次多項式。對於任何 $i in [0, n - 1]$ ，都有 $C(x_i) = 0$ 。
即 $C(x)$ 有 $n$ 個根，這與代數基本定理（一個 $n-1$ 次多項式在複數域上有且僅有 $n-1$ 個根）相矛盾，故 $C(x)$ 並不是一個 $n-1$ 次多項式，推到矛盾。
原命題成立，證畢。

如果我們按照定義求一個多項式的點值表示，時間複雜度為 $O(n^2)$

已知多項式的點值表示，求其係數表示，可以使用插值。樸素的插值演算法時間複雜度為 $O(n^2)$ 。

關於多項式的乘法

已知在一組插值節點 $(x_0,x_1,x_2,...,x_n)$ 中 $A(x),B(x)$ (假設個多項式的項數相同，沒有的視為 $0$ )的點值向量分別為 $(y_{a0},y_{a1},y_{a_2},...,y_{an}),(y_{b0},y_{b1},y_{b_2},...,y_{bn})$ ，那麼 $C(x)=A(x)cdot B(x)$ 的點值表達式可以在 $O(n)$ 的時間內求出，為 $(y_{a0}cdot y_{b0},y_{a1}cdot y_{b1},y_{a_2}cdot y_{b2},...,y_{an}cdot y_{bn})$ 。

因為 $C(x)$ 的項數為 $A(x),B(x)$ 的項數之和。

設 $A(x),B(x)$ 分別有 $n,m$ 項所以我們帶入的插值節點有至少有 $n+m$ 個。

如果我們能快速通過點值表式求出係數表示，那麼就搭起了它們之間的一座橋了。

這也是快速傅里葉變換的基本思路，由係數表達式到點值表達式到結果的點值表達式再到結果的係數表達式。

0x02 關於複數的基本了解

我們把形如 $a+bi$ 這樣的數叫做複數，複數集合用 $C$ 來表示。其中 $a$ 稱為實部 $(real;part)$ ， $b$ 稱為虛部 $(imaginary;part)$ , $i$ 為虛數單位，指滿足 $x^2=-1$ 的一個解 $sqrt{-1}$ ；此外，對於這樣對複數開偶次冪的數叫做虛數 $( imaginary;number)$ .

每一個複數 $a+bi$ 都對應了一個平面上的向量 $(a,b)$ 我們把這樣的平面稱為複平面 $(complex;plane)$ ，它是由水平的實軸與垂直的虛軸建立起來的複數的幾何表示。

故每一個複數唯一對應了一個複平面上的向量，每一個複平面上的向量也唯一對應了一個複數。其中 $0$ 既被認為是實數，也被認為是虛數。

其中複數 $z=a+bi$ 的模長 $left| z ight|$ 定義為 $z$ 在複平面的距離到原點的距離， $left| z ight|=sqrt{a^2+b^2}$ 。幅角 $heta$ 為實軸的正半軸正方向（逆時針）旋轉到 $z$ 的有向角度。

由於虛數無法比較大小。複數之間的大小關係只存在等於與不等於兩種關係，兩個複數相等當且僅當實部虛部對應相等。對於虛部為 $0$ 的複數之間是可以比較大小的，相當於實數之間的比較。

複數之間的運算滿足結合律，交換律和分配律。

由此定義複數之間的運演算法則：

$(a+bi)+(c+di)=(a+c)+(b+d)i$

$(a+bi)-(c+di)=(a-c)+(b-d)i$ $(a+bi)cdot(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+cb)i$

$frac{a+bi}{c+di}=frac{(a+bi)cdot(c-di)}{(c+di)cdot(c-di)}=frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=frac{(ac+bd)}{c^2+d^2}+frac{(bc-ad)i}{c^2+d^2}$

複數運算的加法滿足平行四邊形法則，乘法滿足幅角相加，模長相乘。

對於一個複數 $z=a+bi$ ,它的共軛複數是 $z^{$ , $z^{$ 稱為 $z$ 的復共軛 $(complex;conjugate)$ .

共軛複數有一些性質

$zcdot z^{$

$left| z ight| = left | z^{$

0x03 複數中的單位根

複平面中的單位圓

其中 $ar{OA}$ 單位根，表示為 $e^{i heta}$ ，可知 $e^{i heta}=cos heta+icdot sin heta$

（順便一提著名的歐拉幅角公式 $e^{ipi}+1=0$ 其實是由定義來的...）

將單位圓等分成 $n$ 個部分（以單位圓與實軸正半軸的交點一個等分點），以原點為起點，圓的這 $n$ 個 $n$ 等分點為終點，作出 $n$ 個向量。

其中幅角為正且最小的向量稱為 $n$ 次單位向量，記為 $omega_{n}^{1}$ 。

（有沒有大佬幫我補張圖啊，畫不來）

其餘的 $n-1$ 個向量分別為 $omega_{n}^{2},omega_{n}^{3},......,omega_{n}^{n}$ ，它們可以由複數之間的乘法得來 $w_{n}^{k}=w_{n}^{k-1}cdot w_{n}^{1} (2 leq k leq n)$ 。

容易看出 $w_{n}^{n}=w_{n}^{0}=1$ 。

對於 $w_{n}^{k}$ ,它事實上就是 $e^{2picdot frac{k}{n}i}$ 。

所以 $omega_{n}^{k}=cos(2picdotfrac{k}{n})+icdot sin(2picdotfrac{k}{n})$

關於單位根有兩個性質

性質一（又稱為折半引理）：

$omega_{2n}^{2k}=omega_{n}^{k}$

證明一：

由幾何意義，這兩者表示的向量終點是一樣的。

證明二：

由計算的公式：

$omega_{2n}^{2k}=cos(2pifrac{2k}{2n})+icdot sin(2pifrac{2k}{2n})=cos(2pifrac{k}{n})+icdot sin(2pifrac{k}{n})=omega_{n}^{k}$

其實由此我們可以引申出

$omega_{mn}^{mk}=omega_{n}^{k}$

性質二（又稱為消去引理）

$omega_{n}^{k+frac{n}{2}}=-omega_{n}^{k}$

證明一：

由幾何意義，這兩者表示的向量終點是相反的，左邊較右邊在單位圓上多轉了半圈。

證明二：

由計算的公式：

$omega_{n}^{k+frac{n}{2}}=cos(2pifrac{k+frac{n}{2}}{n})+icdot sin(2pifrac{k+frac{n}{2}}{n})=cos(2pifrac{k}{n}+pi)+icdot sin(2pifrac{k}{n}+pi)=-cos(2pifrac{k}{n})-icdot sin(2pifrac{k}{n})=-omega_{n}^{k}$

最後一步由三角恆等變換得到。

0x04 離散傅里葉變換（Discrete Fourier Transform）

首先我們單獨考慮一個 $n$ 項（ $n=2^x$ ）的多項式 $A(x)$ ，其係數向量為 $(a_0,a_1,a_2,...,a_{n-1})$ 。我們將 $n$ 次單位根的 $0$ ~ $n-1$ 次冪分別帶入 $A(x)$ 得到其點值向量 $(A(w_n^{0}),A(w_n^{1}),A(w_n^{2}),...,A(w_n^{n-1}))$ 。

這個過程稱為離散傅里葉變換（Discrete Fourier Transform）。

如果樸素帶入，時間複雜度也是 $O(n^2)$ 的。

所以我們必須要利用到單位根 $omega$ 的特殊性質。

對於 $A(x)=a_0+a_1cdot x^1+a_2cdot x^2+a_3cdot x^3+...+a_{n-1}cdot x^{n-1}$

考慮將其按照奇偶分組

$A(x)=(a_0+a_2cdot x^2+a_{4}cdot x^{4}...+a_{n-2}cdot x^{n-2})+(a_1cdot x^1+a_3cdot x^3+a_{5}cdot x^{5}+...+a_{n-1}cdot x^{n-1})$

$A(x)=(a_0+a_2cdot x^2+a_{4}cdot x^{4}...+a_{n-2}cdot x^{n-2})+xcdot(a_1+a_3cdot x^2+a_{5}cdot x^{4}+...+a_{n-1}cdot x^{n-2})$

令

$A1(x)=(a_0+a_2cdot x+a_{4}cdot x^{2}...+a_{n-2}cdot x^{frac{n-2}{2}})$

$A2(x)=(a_1+a_3cdot x+a_{5}cdot x^{2}...+a_{n-1}cdot x^{frac{n-2}{2}})$

則可得到

$A(x)=A1(x^2)+xcdot A2(x^2)$

分類討論

設 $0leq kleq frac{n}{2}-1$ ， $kin Z$

$A(omega_{n}^{k})=A1(omega_{n}^{2k})+omega_{n}^{k}cdot A2(omega_{n}^{2k})$

由上文提到的折半引理

$A(omega_{n}^{k})=A1(omega_{frac{n}{2}}^{k})+omega_{n}^{k}cdot A2(omega_{frac{n}{2}}^{k})$

對於 $frac{n}{2}leq k+frac{n}{2}leq n-1$

$A(omega_{n}^{k+frac{n}{2}})=A1(omega_{n}^{2k+n})+omega_{n}^{k+frac{n}{2}}cdot A2(omega_{n}^{2k+n})$

其中 $omega_{n}^{2k+n}=omega_{n}^{2k}cdot omega_{n}^{n}=omega_{n}^{2k}=omega_{frac{n}{2}}^{k}$

由消去引理 $omega_{n}^{k+frac{n}{2}}=-omega_{n}^{k}$

故 $A(omega_{n}^{k+frac{n}{2}})=A1(omega_{frac{n}{2}}^{k})-omega_{n}^{k}cdot A2(omega_{frac{n}{2}}^{k})$

注意， $k$ 與 $k+frac{n}{2}$ 取遍了 $[0,n-1]$ 中的 $n$ 個整數，保證了可以由這 $n$ 個點值反推解出係數（上文已證明）。

於是我們可以知道

如果已知了 $A1(x),A2(x)$ 分別在 $omega_{frac{n}{2}}^{0},omega_{frac{n}{2}}^{1},...,omega_{frac{n}{2}}^{frac{n}{2}-1},$ 的取值，可以在 $O(n)$ 的時間內求出 $A(x)$ 的取值。

而 $A1(x),A2(x)$ 都是 $A(x)$ 一半的規模，顯然可以轉化為子問題遞歸求解。

時間複雜度：

$T(n)=2T(frac{n}{2})+O(n)=O(nlogn)$

0x05 離散傅里葉反變換（Inverse Discrete Fourier Transform）

使用快速傅里葉變換將點值表示的多項式轉化為係數表示，這個過程叫做離散傅里葉反變換（Inverse Discrete Fourier Transform）。

即由 $n$ 維點值向量 $(A(x_0),A(x_1),...,A(x_{n-1}))$ 推出 $n$ 維繫數向量 $(a_0,a_1,...,a_{n-1})$ 。

設 $(d_0,d_1,...,d_{n-1})$ 為 $(a_0,a_1,...,a_{n-1})$ 得到的離散傅里葉變換的結果。

我們構造一個多項式 $F(x)=d_0+d_1cdot x+d_2cdot x^2+...+d_{n-1}cdot x^{n-1}$

設向量 $(c_0,c_1,...,c_{n-1})$ 中

$c_k$ 為 $F(x)$ 在 $x=omega_{n}^{-k}$ 的點值表示

即 $c_k=sum_{i=0}^{n-1}d_icdot(omega_{n}^{-k})^i$ ，

我們考慮對 $d_i$ 進行還原

於是

$c_k=sum_{i=0}^{n-1}[sum_{j=0}^{n-1}a_jcdot(omega_{n}^{i})^j]cdot(omega_{n}^{-k})^i$

由和式的性質

$c_k=sum_{j=0}^{n-1}a_jsum_{i=0}^{n-1}(omega_{n}^{i})^jcdot(omega_{n}^{-k})^i$

$=sum_{j=0}^{n-1}a_jsum_{i=0}^{n-1}(omega_{n}^{i})^{j-k}$

令 $S(j,k)=sum_{i=0}^{n-1}(omega_{n}^{i})^{j-k}$

對其進行化簡

設 $j-k=delta$

則 $S(j,k)=omega_{n}^{0}+omega_{n}^{delta}+omega_{n}^{2delta}+...+omega_{n}^{(n-1)delta}$

其公比為 $omega_{n}^{delta}$

當 $omega_{n}^{delta}=1$ 即 $delta=0$ 時

$S(j,k)=n$ 此時 $delta=0Rightarrow j-k=0 Rightarrow j=k$

當 $omega_{n}^{delta} e1$ 即 $delta e0$ 時

由等比數列求和公式

$S(j,k)=frac{(omega_{n}^{delta})^{n}-1}{omega_{n}^{delta}}=frac{(omega_{n}^{n})^{delta}-1}{omega_{n}^{delta}}=frac{(1)^{delta}-1}{omega_{n}^{delta}}=frac{0}{omega_{n}^{delta}}=0$ ,此時 $j e k$ .

所以

$S(j,k)=[j=k]cdot n$

將 $S(j,k)$ 帶入原式

$c_k=sum_{j=0}^{n-1}a_jcdot S(j,k)=sum_{j=0}^{n-1}a_jcdot [j=k]cdot n=a_kcdot n$

所以 $a_k=frac{c_k}{n}$ .

其中 $a_k$ 為原多項式 $A(x)$ 的係數向量 $(a_0,a_1,...,a_n)$ 中的 $a_k$ .

由此得到：

對於多項式 $A(x)$ 由插值節點 $(omega_{n}^{0},omega_{n}^{1},omega_{n}^{2},...,omega_{n}^{n-1})$ 做離散傅里葉變換得到的點值向量 $(d_0,d_1,...,d_{n-1})$ 。我們將 $(omega_{n}^{0},omega_{n}^{-1},omega_{n}^{-2},...,omega_{n}^{-(n-1)})$ 作為插值節點， $(d_0,d_1,...,d_{n-1})$ 作為係數向量，做一次離散傅里葉變換得到的向量每一項都除以 $n$ 之後得到的 $(frac{c_0}{n},frac{c_1}{n},...,frac{c_{n-1}}{n})$ 就是多項式的係數向量 $(a_0,a_1,...,a_{n-1})$ 。

注意到 $omega_{n}^{-k}$ 是 $omega_{n}^{k}$ 的共軛複數。

這個過程稱為離散傅里葉反變換。

0x06 關於FFT在C++的實現

首先要解決複數運算的問題，我們可以使用C++STL自帶的 $std :: complex < T >$ 依照精度要求 $T$ 一般為 $double,long;double$ 。

也可以自己封裝，下面是我封裝的複數類。

struct Complex { double r, i ; Complex ( ) { } Complex ( double r, double i ) : r ( r ), i ( i ) { } inline void real ( const double& x ) { r = x ; } inline double real ( ) { return r ; } inline Complex operator + ( const Complex& rhs ) const { return Complex ( r + rhs.r, i + rhs.i ) ; } inline Complex operator - ( const Complex& rhs ) const { return Complex ( r - rhs.r, i - rhs.i ) ; } inline Complex operator * ( const Complex& rhs ) const { return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ; } inline void operator /= ( const double& x ) { r /= x, i /= x ; } inline void operator *= ( const Complex& rhs ) { *this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ; } inline void operator += ( const Complex& rhs ) { r += rhs.r, i += rhs.i ; } inline Complex conj ( ) { return Complex ( r, -i ) ; }} ;

我們由上面的分析可以得到這個遞歸的寫法。

bool inverse = false ;inline Complex omega ( const int& n, const int& k ) { if ( ! inverse ) return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ) ; return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ).conj ( ) ;}inline void fft ( Complex *a, const int& n ) { if ( n == 1 ) return ; static Complex buf [N] ; const int m = n >> 1 ; for ( int i = 0 ; i < m ; ++ i ) { buf [i] = a [i << 1] ; buf [i + m] = a [i << 1 | 1] ; } memcpy ( a, buf, sizeof ( int ) * ( n + 1 ) ) ; Complex *a1 = a, *a2 = a + m; fft ( a1, m ) ; fft ( a2, m ) ; for ( int i = 0 ; i < m ; ++ i ) { Complex t = omega ( n, i ) ; buf [i] = a1 [i] + t * a2 [i] ; buf [i + m] = a1 [i] - t * a2 [i] ; } memcpy ( a, buf, sizeof ( int ) * ( n + 1 ) ) ;}

但是這樣的 $FFT$ 要用到輔助數組，並且常數比較大。

能不能優化呢？

我們把每一次分組的情況推演出來

遞歸分類的每一層

觀察到每一個位置的數其實都是原來位置上的數的二進位後 $log_{2}n$ 位 $reverse$ 了一下。

於是我們可以想，先將原數組調整成最底層的位置（很好調整吧）。

然後從倒數第二層由底向上計算。

這就是我們一般用來實現 $FFT$ 的 $Cooley-Tukey$ 演算法。

考慮怎麼合併？

在 $Cooley-Tukey$ 演算法中，合併操作被稱作是蝴蝶操作。

慮合併兩個子問題的過程，這一層有 $n$ 項需要處理。假設 $A_1(omega_{ frac{n}{2} } ^ k)$ 和 $A_2(omega_{ frac{n}{2} } ^ k)$ 分別存在 $a(k)$ 和 $a(frac{n}{2} + k)$ 中， $A(omega_n ^ {k})$ 和 $A(omega_n ^ {k + frac{n}{2} })$ 將要被存放在 $buf(k)$ 和 $buf(frac{n}{2} + k)$ 中，合併的單位操作可表示為

$buf(k):=a(k)+omega_{n}^{k}cdot a(k+frac{n}{2})$

$buf(k+frac{n}{2}):=a(k)-omega_{n}^{k}cdot a(k+frac{n}{2})$

只要將合併順序換一下，再加入一個臨時變數，合併過程就可以在原數組中進行。

令 $t:=omega_{n}^{k}cdot a(k+frac{n}{2})$

合併過程如下：

$a(k+frac{n}{2}):=a(k)-t$

$a(k):=a(k)+t$ 。

至此，我們可以給出 $Cooley-Tukey$ 演算法的實現。

struct FastFourierTransform { Complex omega [N], omegaInverse [N] ; void init ( const int& n ) { for ( int i = 0 ; i < n ; ++ i ) { omega [i] = Complex ( cos ( 2 * PI / n * i), sin ( 2 * PI / n * i ) ) ; omegaInverse [i] = omega [i].conj ( ) ; } } void transform ( Complex *a, const int& n, const Complex* omega ) { for ( int i = 0, j = 0 ; i < n ; ++ i ) { if ( i > j ) std :: swap ( a [i], a [j] ) ; for( int l = n >> 1 ; ( j ^= l ) < l ; l >>= 1 ) ; } for ( int l = 2 ; l <= n ; l <<= 1 ) { int m = l / 2; for ( Complex *p = a ; p != a + n ; p += l ) { for ( int i = 0 ; i < m ; ++ i ) { Complex t = omega [n / l * i] * p [m + i] ; p [m + i] = p [i] - t ; p [i] += t ; } } } } void dft ( Complex *a, const int& n ) { transform ( a, n, omega ) ; } void idft ( Complex *a, const int& n ) { transform ( a, n, omegaInverse ) ; for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ; }} fft ;