【轉】用C語言實現FFT算法

本文轉載自查看原文 2018-10-15 23:39 11874 數據結構/算法

傅里葉變換

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

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

考慮到兩個多項式 A(x),B(x) 的乘積 C(x) ，假設 A(x) 的項數為，其系數構成的維向量為 $(a_0,a_1,a_2,...,a_{n-1})$ ， B(x) 的項數為，其系數構成的維向量為 $(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_i\cdot x^i$

多項式的點值表示法

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

其中

$y_i=\sum_{j=0}^{n-1}a_j\cdot x_i^j$

定理:

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

證明：

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

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

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

關於多項式的乘法

已知在一組插值節點 (x_0,x_1,x_2,...,x_n) 中 A(x),B(x) (假設個多項式的項數相同，沒有的視為 )的點值向量分別為 $(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 這樣的數叫做復數，復數集合用來表示。其中稱為實部 $(real\;part)$ ，稱為虛部 $(imaginary\;part)$ , 為虛數單位，指滿足 x^2=-1 的一個解 $\sqrt{-1}$ ；此外，對於這樣對復數開偶次冪的數叫做虛數 $( imaginary\;number)$ .

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

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

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

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

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

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

(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^{'}=a-bi$ , $z^{'}$ 稱為的復共軛 $(complex\;conjugate)$ .

共軛復數有一些性質

$z\cdot z^{'}=a^2+b^2$

$\left| z \right| = \left | z^{'} \right|$

0x03 復數中的單位根

復平面中的單位圓

其中 $\bar{OA}$ 單位根，表示為 $e^{i\theta}$ ，可知 $e^{i\theta}=cos\theta+i\cdot sin\theta$

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

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

其中幅角為正且最小的向量稱為次單位向量，記為 $\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^{2\pi\cdot \frac{k}{n}i}$ 。

所以 $\omega_{n}^{k}=cos(2\pi\cdot\frac{k}{n})+i\cdot sin(2\pi\cdot\frac{k}{n})$

關於單位根有兩個性質

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

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

證明一：

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

證明二：

由計算的公式：

$\omega_{2n}^{2k}=cos(2\pi\frac{2k}{2n})+i\cdot sin(2\pi\frac{2k}{2n})=cos(2\pi\frac{k}{n})+i\cdot sin(2\pi\frac{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(2\pi\frac{k+\frac{n}{2}}{n})+i\cdot sin(2\pi\frac{k+\frac{n}{2}}{n})=cos(2\pi\frac{k}{n}+\pi)+i\cdot sin(2\pi\frac{k}{n}+\pi)=-cos(2\pi\frac{k}{n})-i\cdot sin(2\pi\frac{k}{n})=-\omega_{n}^{k}$

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

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

首先我們單獨考慮一個項（ n=2^x ）的多項式 A(x) ，其系數向量為 $(a_0,a_1,a_2,...,a_{n-1})$ 。我們將次單位根的 ~ 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_1\cdot x^1+a_2\cdot x^2+a_3\cdot x^3+...+a_{n-1}\cdot x^{n-1}$

考慮將其按照奇偶分組

$A(x)=(a_0+a_2\cdot x^2+a_{4}\cdot x^{4}...+a_{n-2}\cdot x^{n-2})+(a_1\cdot x^1+a_3\cdot x^3+a_{5}\cdot x^{5}+...+a_{n-1}\cdot x^{n-1})$

$A(x)=(a_0+a_2\cdot x^2+a_{4}\cdot x^{4}...+a_{n-2}\cdot x^{n-2})+x\cdot(a_1+a_3\cdot x^2+a_{5}\cdot x^{4}+...+a_{n-1}\cdot x^{n-2})$

令

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

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

則可得到

$A(x)=A1(x^2)+x\cdot A2(x^2)$

分類討論

設 $0\leq k\leq \frac{n}{2}-1$ ， $k\in 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+\frac{n}{2}$ 取遍了 [0,n-1] 中的個整數，保證了可以由這個點值反推解出系數（上文已證明）。

於是我們可以知道

如果已知了 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）。

即由維點值向量 $(A(x_0),A(x_1),...,A(x_{n-1}))$ 推出維系數向量 $(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_1\cdot x+d_2\cdot 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_i\cdot(\omega_{n}^{-k})^i$ ，

我們考慮對 d_i 進行還原

於是

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

由和式的性質

$c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{i})^j\cdot(\omega_{n}^{-k})^i$

$=\sum_{j=0}^{n-1}a_j\sum_{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}^{2\delta}+...+\omega_{n}^{(n-1)\delta}$

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

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

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

當 $\omega_{n}^{\delta}\ne1$ 即 $\delta\ne0$ 時

由等比數列求和公式

$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\ne k$ .

所以

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

將 S(j,k) 帶入原式

$c_k=\sum_{j=0}^{n-1}a_j\cdot S(j,k)=\sum_{j=0}^{n-1}a_j\cdot [j=k]\cdot n=a_k\cdot 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})$ 作為系數向量，做一次離散傅里葉變換得到的向量每一項都除以之后得到的 $(\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 > 依照精度要求一般為 $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 算法中，合並操作被稱作是蝴蝶操作。

慮合並兩個子問題的過程，這一層有項需要處理。假設 $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 ;

注意代碼中的 omega[k] 為 $\omega_{n}^{k}$ ,而在代碼中需要得到的是 $\omega_{l}^{k}$ 。