【转】用C语言实现FFT算法-白红宇

【转】用C语言实现FFT算法

阅读量：5031 次

发布时间：2019-06-12

本文共 6009 字，大约阅读时间需要 20 分钟。

傅里叶变换

快速傅里叶变换（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_i\cdot 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_j\cdot 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 \right|$ 定义为 $z$ 在复平面的距离到原点的距离， $\left| z \right|=\sqrt{a^2+b^2}$ 。幅角 $\theta$ 为实轴的正半轴正方向（逆时针）旋转到 $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^{'}=a-bi$ , $z^{'}$ 称为 $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$ 其实是由定义来的...）

将单位圆等分成 $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^{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$ 项（ $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_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$ 与 $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_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})$ 作为系数向量，做一次离散傅里叶变换得到的向量每一项都除以 $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 ;

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