FFT——快速傅里叶变换

本文主要是介绍FFT——快速傅里叶变换，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

FFT——快速傅里叶变换

卷积

一般来说在计算机上处理卷积通常是离散的，所以这里只介绍离散卷积

有两个序列\(\{a_n\},\{b_n\}\)​​，若将这两个序列按以下方式生成一个新序列\(\{c_n\}\)​

\[c_k=\sum\limits_{i=-\infty}^{+\infty} a_i\cdot b_{k-i} \]

则新序列\(\{c_n\}\)​​称为\(\{a_n\}\)​​和\(\{b_n\}\)​​​​​的卷积。用文字描述即为：新序列的第\(n\)项为原序列下标之和为\(n\)的乘积之和。

计算机能处理的是有限数据，上面的形式可以类推到有限序列：

有两个序列，\(a_0,a_1,...,a_{n-1}\)和\(b_0,b_1,...,b_{m-1}\)，新序列为\(c_0,c_1,...,c_{n+m-2}\)，生成方式如下式

\[c_k=\sum\limits_{\substack{0\leq i\leq n-1\\0\leq j\leq m-1\\ i+j=k}}a_i b_j \]

我们主要讨论的就是这种形式的卷积。

多项式乘法

一个含有\(n\)项的多项式一般表示为

\[f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \]

假如你还有另外一个含有\(m\)项的多项式\(g(x)\)​

\[g(x)=b_0+b_1x+b_2x^2+\cdots+b_{m-1}x^{m-1} \]

那么两个多项式相乘会得到一个新多项式\(h(x)\)​，它有\(n+m-1\)项

\[h(x)=f(x)\cdot g(x)=c_0+c_1x+c_2x^2+\cdots+c_{n+m-2}x^{n+m-2} \]

我们来观察一下这三个多项式的系数序列，发现它们满足

\[c_k=\sum\limits_{\substack{0\leq i\leq n-1\\0\leq j\leq m-1\\ i+j=k}}a_i b_j \]

发现形式和卷积一样，所以多项式相乘相当于对系数做卷积。

多项式的表示方式

系数表示

顾名思义，就是通过系数序列来表示一个多项式，很好理解。例如上面的\(f(x)\)​​可以用\(a_0,a_1,a_2\cdots,a_{n-1}\)​​​这个序列来表示。这也是通常用来表示多项式的方式。

点值表示

假如你知道\(n\)​个不同的点，这\(n\)​个点可以确定一个\(n-1\)​​次的多项式\(f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}\)​。

即已知\((x_1,f(x_1)),(x_2,f(x_2)),\cdots,(x_n,f(x_n))\)​，你可以确定\(a_0,a_1,a_2,\cdots,a_{n-1}\)​的值（证明可用线代知识）。

所以足够多的点也可以表示一个多项式。

点值表示有什么用呢？

假如你有一个\(n-1\)次多项式\(f(x)\)和\(m-1\)次多项式\(g(x)\)。你可以在这两条多项式的曲线上找出无数多个点。比如你可以找出\(2(n+m-1)\)个点，\((x_1,f(x_1)),(x_2,f(x_2)),\cdots,(x_{n+m-1},f(x_{n+m-1}))\)以及\((x_1,g(x_1)),(x_2,g(x_2)),\cdots,(x_{n+m-1},g(x_{n+m-1}))\)

现在问你根据这些点，是否能很轻易的找出\(f(x)\cdot g(x)\)​曲线上的\(n+m-1\)​个点呢？

显然是很容易的，\((x_1,f(x_1)\cdot g(x_1)),(x_2,f(x_2)\cdot g(x_2)),\cdots,(x_{n+m-1},f(x_{n+m-1})\cdot g(x_{n+m-1}))\)​

而根据上面的理论，这\(n+m-1\)个点就唯一确定了\(f(x)\cdot g(x)\)​​​，我们就轻松求出了多项式乘积的点值表示，这个过程是线性复杂度。

所以我们可以发现点值表示的优点就是乘积非常快速便捷，缺点就是没有系数表示直观。而多项式乘积相当于对系数做卷积，所以如果我们能快速进行系数表示和点值表示的转化，就可以快速地计算卷积。

离散傅里叶变换及其逆变换

离散傅里叶变换

我们很容易想到一个朴素的将多项式的系数表示转化成点值表示的方法，那就是随便找若干个点，然后一一代入式子计算，但是很遗憾这种方法是\(O(n^2)\)​的，和暴力方法计算卷积的复杂度是一致的，所以这个方法不够高效。

我们能否通过某些对称性来减少计算量？最容易想到的就是\(1\)​​和\(-1\)​​​，我们可以通过计算奇数项系数和以及偶数项系数和，同时得出\(1\)​​和\(-1\)​​​​处的点值表示，但是这只有两个点，显然不够。于是我们想到了模为\(1\)​​的复数，想到了共轭复数，想到了单位根，设单位根\(\omega_n=\mathrm{e}^\frac{2\pi \mathrm{i}}{n}\)​​​​​​。我们不妨让将\(\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}\)​​​​​代入，计算点值表示，这就是离散傅里叶变换（当然可以选择朴素做法暴力代入单位根，这也是离散傅里叶变换，但是如果用下面的做法，就叫快速傅里叶变换）

显然，单位根有如下性质。

\[\omega_{2n}^{k+n}=-\omega_{2n}^{k}\\ \omega_{2n}^{2k}=\omega_{n}^{k} \]

为了简单起见，我们认为\(n\)​为\(2\)​的幂次，那些多出来的项，系数认为是\(0\)​​。

我们将多项式拆成奇数项和偶数项

设

\[f_{012\cdots n-1}(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\\ f_{024\cdots n-2}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n-2}{2}}\\ f_{135\cdots n-1}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n-2}{2}} \]

显然有

\[f_{012\cdots n-1}(x)=f_{024\cdots n-2}(x^2)+xf_{135\cdots n-1}(x^2) \]

根据性质，容易发现，当\(k<\frac{n}{2}\)​时，有

\[f_{012\cdots n-1}(\omega_n^k)=f_{024\cdots n-2}(\omega_{\frac{n}{2}}^{k})+\omega_n^kf_1(\omega_{\frac{n}{2}}^{k})\\ f_{012\cdots n-1}(\omega_n^{k+\frac{n}{2}})=f_{024\cdots n-2}(\omega_{\frac{n}{2}}^{k})-\omega_n^k f_{135\cdots n-1}(\omega_{\frac{n}{2}}^{k})\\ \]

什么意思呢，就是我只要知道\(f_{024\cdots n-2}(\omega_{\frac{n}{2}}^k)\)和\(f_{135\cdots n-1}(\omega_{\frac{n}{2}}^{k})\)，我就能同时知道\(f_{012\cdots n-1}(\omega_n^k)\)和\(f_{012\cdots n-1}(\omega_n^{k+\frac{n}{2}})\)​

可以采用分治的方式计算，递归出口也很显然\(f_i(\omega_1^0)=a_i\)​。​

然而递归常数太大了，我们可以尝试画一下递归树，看它每一个节点分别对哪些系数进行处理

我们寻找这些系数下标的规律，由于第一次是偶数放前面，奇数放后面，所以观察它们下标的二进制

第一层，把最低位为0的那些提出来，这部分放在前面，最低位为1的那些提出来，这部分放在后面。

可以发现，若忽略最低位的话，两部分依旧各自从0递增（相当于少一位二进制），前半部分最低位是0，后半部分最低位是1

和原来比，原来是忽略最高位，两部分依旧各自从0递增（相当于少一位二进制），前半部分最高位是0，后半部分最高位是1

这样的话，如果一直递归下去，仔细想一下，发现会将每个下标的二进制倒置（发现不了的话就按照这个思路多写几层，应该就能发现了）

对于\(i\in[0,2^n)\)​​，二进制倒置记为\(rev_i\)​。​​

则显然有这个递推式，

rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (n - 1))

这样可以在\(O(n)\)​时间内预处理出下标反转，

根据预处理交换之后，即可原地进行变换，复杂度是\(O(n\log n)\)​的。（这部分可以配合下面的代码进行理解）

离散傅里叶逆变换

上面解决了从系数表示到点值表示的问题，那反过来呢？似乎就是解一个方程组，为了方便表示，设\(y_i=f(\omega_n^i)\)​。

我们可以把方程组写成矩阵形式。

\[\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{bmatrix}= \begin{bmatrix} (\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^0)^{n-1}\\ (\omega_n^1)^0&(\omega_n^1)^1&\cdots&(\omega_n^1)^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ (\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\cdots&(\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix} \]

设

\[A=\begin{bmatrix} (\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^0)^{n-1}\\ (\omega_n^1)^0&(\omega_n^1)^1&\cdots&(\omega_n^1)^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ (\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\cdots&(\omega_n^{n-1})^{n-1} \end{bmatrix}\\ Y=\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{bmatrix}, X=\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix} \]

容易发现\(X=A^{-1}Y\)​​，我们只要想办法求出\(A\)​的逆矩阵，然后再用离散傅里叶变换即可。

设\(A\)​的第\(i\)​行，第\(j\)​列（为了方便从\(0\)​开始）为\(d_{ij}\)​​

则\(d_{ij}=\omega_n^{ij}\)，我们不妨试试看构造矩阵\(B\)，它的第\(i\)行，第\(j\)列（为了方便从\(0\)开始）为\(e_{ij}=\omega_n^{-ij}\)​

则\(AB\)​的第\(i\)​行，第\(j\)​列是

\[\begin{aligned} x_{ij}&=\sum\limits_{k=0}^{n-1}d_{ik}e_{kj}\\ &=\sum\limits_{k=0}^{n-1}\omega_n^{(i-j)k}\\ &= \begin{cases} n,&i=j\\ \frac{1-\omega_n^{(i-j)n}}{1-\omega_n^{i-j}},&i\neq j \end{cases}\\ &=\begin{cases} n,&i=j\\ \frac{1-1}{1-\omega_n^{i-j}},&i\neq j \end{cases}\\ &=\begin{cases} n,&i=j\\ 0,&i\neq j \end{cases} \end{aligned} \]

干得漂亮，前面再除一个\(n\)就是单位矩阵了

所以\(A\)的逆矩阵，就是将每一个元素取一个共轭，前面再除以\(n\)，然后利用\(A^{-1}\)对\(Y\)​进行离散傅里叶变换即可

模板

const double pi = acos(-1);
struct Complex {
    double re, im;
    Complex operator+(const Complex& rhs) const {
        return {re + rhs.re, im + rhs.im};
    }
    Complex operator-(const Complex& rhs) const {
        return {re - rhs.re, im - rhs.im};
    }
    Complex operator*(const Complex& rhs) const {
        return {re * rhs.re - im * rhs.im, re * rhs.im + im * rhs.re};
    }
    Complex operator/(double rhs) const { return {re / rhs, im / rhs}; }
};

int k, l, n;
void init() {
    k = 1, l = 0;
    while (k < n) k <<= 1, l++; // k=不小于n的最小的2^l
    rev[0] = 0;
    for (int i = 1; i < k; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
}


// f=1，正向，f=-1，逆向
void FFT(Complex* a, int n, int f) {
    for (int i = 0; i < n; i++)
        if (i < rev[i]) swap(a[i], a[rev[i]]); // 防止多次交换
    for (int i = 1; i < n; i <<= 1) {
        Complex d_omega = {cos(pi / i), f * sin(pi / i)};
        for (int j = 0; j < n; j += i << 1) {
            Complex omega = {1, 0};
            for (int k = 0; k < i; k++) {
                Complex x = a[k + j], y = omega * a[k + j + i];
                a[k + j] = x + y;
                a[k + j + i] = x - y;
                omega = omega * d_omega;
            }
        }
    }
    if (f == -1)
        for (int i = 0; i < n; i++) a[i] = a[i] / n;
}

这篇关于FFT——快速傅里叶变换的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！