OJ：快速傅里叶变换（FFT）学习笔记

最新推荐文章于 2020-03-19 12:42:42 发布
最新推荐文章于 2020-03-19 12:42:42 发布
阅读量452 收藏 2
点赞数
分类专栏： OJ笔记
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/WinterShiver/article/details/104251042
版权

FFT 快速傅里叶变换 学习笔记

前言

FFT是一个快速(O(nlogn))求两个多项式乘积的算法。
在阅读本文之前，请先阅读1.本文是对1的解释和补充，主要贡献有两点：

  • 把原文没说清的理论部分补上；
  • 通过递归，把原文的思路重新理了一遍，使理论和代码都更加好理解。

另外，本文的附录里面只放了源代码，且在文中有提到。

文章目录

理论：IFFT的理论基础

只要搞懂以下5个topic，就可以明白快速傅里叶变换的代码逻辑：

1.多项式的系数表示法和点值表示法
2.n次单位根
3.DFT
4.FFT
5.IFFT

理论部分的前四点构成了FFT的理论基础。 对于这四点，1都有，且很好理解。请读者阅读原文熟悉前四点理论。

下面主要讲一下第5点（IFFT的理论基础）。

FFT本质上解决了一个乘法问题

以下所有记号和1对应，并记 w i = ω n i = e x p ( i ∗ 2 π i / n ) w_i = \omega_n^i = exp(\mathbf{i}*2\pi i/n) wi​=ωni​=exp(i∗2πi/n)，对应n次单位根。（ i \mathbf{i} i是虚数单位。）

记n次单位根的幂矩阵W，系数向量C，值向量V如下：

W = [ W i , j = ( w i ) j ] n ∗ n = [ ( w 0 ) 0 ⋯ ( w 0 ) n − 1 ( w 1 ) 0 ⋯ ( w 1 ) n − 1 ⋮ ⋱ ⋮ ( w n − 1 ) 0 ⋯ ( w n − 1 ) n − 1 ] W=[W_{i,j}=(w_i)^j]_{n*n}= \left[ \begin{aligned} (w_0)^0 & \cdots & (w_0)^{n-1} \\ (w_1)^0 & \cdots & (w_1)^{n-1} \\ \vdots & \ddots & \vdots \\ (w_{n-1})^0 & \cdots & (w_{n-1})^{n-1} \\ \end{aligned} \right] W=[Wi,j​=(wi​)j]n∗n​=⎣⎢⎢⎢⎢⎢⎡​(w0​)0(w1​)0⋮(wn−1​)0​⋯⋯⋱⋯​(w0​)n−1(w1​)n−1⋮(wn−1​)n−1​⎦⎥⎥⎥⎥⎥⎤​

C = [ a 0 a 1 ⋮ a n − 1 ] C = \left[ \begin{aligned} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{aligned} \right] C=⎣⎢⎢⎢⎢⎢⎡​a0​a1​⋮an−1​​⎦⎥⎥⎥⎥⎥⎤​

V = [ A ( x 0 ) A ( x 1 ) ⋮ A ( x n − 1 ) ] V = \left[ \begin{aligned} A(x_0) \\ A(x_1) \\ \vdots \\ A(x_{n-1}) \end{aligned} \right] V=⎣⎢⎢⎢⎢⎢⎡​A(x0​)A(x1​)⋮A(xn−1​)​⎦⎥⎥⎥⎥⎥⎤​

其中 W C = V WC=V WC=V.

FFT本质上解决的问题是已知W和C，计算V的乘法问题，对这个问题给出了一个 O ( n l o g n ) O(nlogn) O(nlogn)的方法。

IFFT可以被表示成一个类似的乘法问题

IFFT可以表示成一个类似的乘法问题，并通过相似的方式解决。这一部分参考了2.

我们首先求W的逆矩阵。记 v i = c o n j ( w i ) = e x p ( − i ∗ 2 π i / n ) v_i = conj(w_i) = exp(-\mathbf{i}*2\pi i/n) vi​=conj(wi​)=exp(−i∗2πi/n)，即 w i w_i wi​的共轭。（ i \mathbf{i} i是虚数单位。）

记矩阵W’为

W ′ = [ W i , j ′ = ( v i ) j ] n ∗ n = [ ( v 0 ) 0 ⋯ ( v 0 ) n − 1 ( v 1 ) 0 ⋯ ( v 1 ) n − 1 ⋮ ⋱ ⋮ ( v n − 1 ) 0 ⋯ ( v n − 1 ) n − 1 ] W'=[W'_{i,j}=(v_i)^j]_{n*n}= \left[ \begin{aligned} (v_0)^0 & \cdots & (v_0)^{n-1} \\ (v_1)^0 & \cdots & (v_1)^{n-1} \\ \vdots & \ddots & \vdots \\ (v_{n-1})^0 & \cdots & (v_{n-1})^{n-1} \\ \end{aligned} \right] W′=[Wi,j′​=(vi​)j]n∗n​=⎣⎢⎢⎢⎢⎢⎡​(v0​)0(v1​)0⋮(vn−1​)0​⋯⋯⋱⋯​(v0​)n−1(v1​)n−1⋮(vn−1​)n−1​⎦⎥⎥⎥⎥⎥⎤​

发现 W W ′ = W ′ W = n I WW'=W'W=nI WW′=W′W=nI（易证）,于是W的逆矩阵是 W ′ n \frac{W'}{n} nW′​,有 W ′ V = n C W'V=nC W′V=nC.因此，只要在FFT的实现算法中把n次单位根改成其共轭，就可以从向量V用 O ( n l o g n ) O(nlogn) O(nlogn)时间求解到向量C，即实现IFFT算法。

综上所述，实现IFFT和实现FFT使用的是同一套逻辑，只是使用了不同的n次单位根。

FFT算法的基础实现

以下是FFT算法的基础实现，fft_logic是FFT和IFFT的共同主逻辑，fft和ifft是给fft_logic套的两个壳。结合注释阅读fft_logic的代码即可。
（完整的程序见附录1，包括测试用的主函数。complex类的使用见3.）

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支持位数
complex<double> b[N]; // 用来充当临时调整空间的数组

void fft_logic(complex<double> *a, int n, int inv){

    // 参数：
    // 当inv = 1时，a是系数多项式，n是当前数组长度（2的幂次），函数效果是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当前数组长度（2的幂次），函数效果是原地变成系数多项式，但是所得的系数是n倍，需要在包装的函数中进行调整

    if (n == 1) return; // 为什么？因为omega_1^0=1，点值多项式和系数多项式的表示完全一致。

    // 利用B暂存和写回，把a的顺序调整为 a[0] a[2] .. a[n-2] a[1] a[3] .. a[n-1]，前后两半

    for(int i = 0; i < n/2; i ++){
        b[i]       = a[i * 2];
        b[i + n/2] = a[i * 2 + 1];
    }
    for(int i = 0; i < n; i ++)
        a[i] = b[i];

    // 分治求A1和A2

    fft_logic(a, n/2, inv);
    fft_logic(a + n/2, n/2, inv);

    // 通过A1和A2，计算A

    double unit_rad = 2 * pi / n; // 单位角幅度值

    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i 
        complex<double> tmp1 = a[i];
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }

}

void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}

void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}

如同原文作者所说，这种方法的缺点就是常数太大需要优化…

FFT的迭代版本：优化常数

基于递归的版本

针对上一节中算法的不足，注意到ffp_logic函数中“每次递归都重排”，我们将其优化成“预先重排后直接按顺序处理”。

优化的思路是，针对原来算法中每一次A1和A2的穿插，其实可以通过提前重排原始数组来代替。拿来原文的一张图来说明：

n=8时，原来的顺序最终被穿插成了0 4 2 6 1 5 3 7,n=4时则是0 2 1 3.这个数列仅由n决定，所以我们可以用一个函数生成这样的穿插数列。

void fft_rearrange_decidesequence_logic(int *rev, int n){

    // 给定数组rev和数组长度n，函数的功能是将所需的顺序写入数组，比如n = 4时将顺序 0 2 1 3 写入数组rev

    if(n == 1){
        rev[0] = 0;
        return;
    }

    // 获得 n/2 时的顺序,暂时放在rev的后半

    fft_rearrange_decidesequence_logic(rev + n/2, n/2);

    // 利用 n/2 时的顺序构造 n 时的顺序

    for(int i = 0; i < n/2; i ++){
        rev[i] = 2 * rev[i + n/2];
        rev[i + n/2] = 2 * rev[i + n/2] + 1;
    }

}

fft_rearrange_decidesequence_logic函数用递归的方式，自然地把穿插的逻辑描述出来，结果得到正确的重排顺序。

void fft_rearrange_logic(complex<double> *a, int n){

    // 按照算法，把a重新排列为可以直接递归向上的顺序

    // 计算bit: 满足pow(2, bit) = n

    int bit = 0;
    while((1 << bit) < n) bit ++;
    
    // rev: 确定a最终rearrange的位置序列

    int* rev = new int[n];
    fft_rearrange_decidesequence_logic(rev, n);

    // 按照rev序列调整a的顺序

    complex<double>* tmp_a = new complex<double>[n];

    for(int i = 0; i < n; i ++) tmp_a[i] = a[rev[i]];
    for(int i = 0; i < n; i ++) a[i] = tmp_a[i];

    delete[] rev;
}

fft_rearrange_logic根据n值获得重排序列，并将输入进行重排。

void fft_merge_logic(complex<double> *a, int n, int inv){

    if(n == 1) return;

    fft_merge_logic(a,       n/2, inv);
    fft_merge_logic(a + n/2, n/2, inv);

    double unit_rad = 2 * pi / n;

    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i
        complex<double> tmp1 = a[i]; 
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }

}

void fft_logic(complex<double> *a, int n, int inv){

    // 参数：
    // 当inv = 1时，a是系数多项式，n是当前数组长度（2的幂次），函数效果是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当前数组长度（2的幂次），函数效果是原地变成系数多项式，但是所得的系数是n倍，需要在包装的函数中进行调整

    // a中元素的顺序调整为算法要求的顺序

    fft_rearrange_logic(a, n);
    
    // 调整顺序之后的a进行合并

    fft_merge_logic(a, n, inv);
}

在将重排环节统一之后，fft_logic函数如上所示。在这里，fft_merge_logic把上一节中重排后合并的逻辑抽了出来。

完整的，可运行的代码见附录2.

这一节的目的是对于原文改良速度的FFT算法写出一个易于理解，可读性强的实现。

rearrange逻辑：递归变非递归

比较一下附录2和3的代码，发现fft_rearrange_decidesequence_logic被砍掉了，换成了一个简洁的循环。
（另外，按照rev序列调整a的顺序的部分，尊重原文1修改了一下。我觉得两种方式效率差异不大，而且也都比较容易理解。）

因为有一个结论：每个位置分治后的最终位置为其二进制翻转后得到的位置（即老生常谈的蝴蝶算法），所以可以用这个公式直接生成rev，而不需要写递归。

void fft_rearrange_logic(complex<double> *a, int n){

    // 按照算法，把a重新排列为可以直接递归向上的顺序

    // 计算bit: 满足pow(2, bit) = n

    int bit = 0;
    while((1 << bit) < n) bit ++;
    
    // rev: 确定a最终rearrange的位置序列

    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }

    // 按照rev序列调整a的顺序

    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
    }

    delete[] rev;
}

完整的代码如附录3所示。

去除算法中的所有递归

我们考虑递归的结构（原谅我懒得再画一张新图）：

8依赖2个4的完成，每个4依赖其下面2个2的完成，每个2依赖下面2个1的完成。因此，可以通过一个循环，先完成8个1（不需要操作），再完成4个2，再完成2个4，最后完成8.

把这个逻辑用循环写出来，就得到了刨除递归的算法。这个算法和原文1的最终板子是对应的。

void fft_logic(complex<double> *a,int n,int inv){

    // bit: pow(2, bit) = n

    int bit = 0;
    while((1 << bit) < n) bit ++;
    
    // rev: 确定最终rearrange的位置序列

    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }

    // 按照rev序列调整a的顺序

    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
    }

    delete[] rev;
    
    // 调整顺序之后的a进行合并

    for (int mid=1;mid<n;mid*=2){
        // mid循环内，准备把两个长度为mid的序列合并成2 * mid的序列

        double unit_rad = 2 * pi / (2 * mid); // 单位角幅度值

        for (int i=0;i<n;i+=mid*2){
            // i循环内，把位置在i~i+mid*2位置的两个长度为mid的序列合并成2 * mid的序列

            for (int j = 0; j < mid; j ++){
                // j循环内的逻辑和
                complex<double> x(cos(j * unit_rad), inv*sin(j * unit_rad)); // x = omega_n^i
                complex<double> tmp1 = a[i+j]; 
                complex<double> tmp2 = x * a[i+j+mid];
                a[i+j]     = tmp1 + tmp2;
                a[i+j+mid] = tmp1 - tmp2;
            }
        }
    }
}

完整的代码如附录4所示。

附录

附录1

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支持位数
complex<double> b[N]; // 用来充当临时调整空间的数组

void fft_logic(complex<double> *a, int n, int inv){

    // 参数：
    // 当inv = 1时，a是系数多项式，n是当前数组长度（2的幂次），函数效果是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当前数组长度（2的幂次），函数效果是原地变成系数多项式，但是所得的系数是n倍，需要在包装的函数中进行调整

    if (n == 1) return; // 为什么？因为omega_1^0=1，点值多项式和系数多项式的表示完全一致。

    // 利用B暂存和写回，把a的顺序调整为 a[0] a[2] .. a[n-2] a[1] a[3] .. a[n-1]，前后两半

    for(int i = 0; i < n/2; i ++){
        b[i]       = a[i * 2];
        b[i + n/2] = a[i * 2 + 1];
    }
    for(int i = 0; i < n; i ++)
        a[i] = b[i];

    // 分治求A1和A2

    fft_logic(a, n/2, inv);
    fft_logic(a + n/2, n/2, inv);

    // 通过A1和A2，计算A

    double unit_rad = 2 * pi / n; // 单位角幅度值

    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i 
        complex<double> tmp1 = a[i];
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }

}

void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}

void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}

// 主函数测试

complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘结果

int main(){

    // 两个相乘多项式的初始化。注意：两个相乘多项式和最后的乘积多项式需要有相同的项数

    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }

    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4

    // F，乘，F逆

    fft(A1, N);
    fft(A2, N);

    /*

    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;

    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;

    */

    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];

    ifft(C, N);

    // 输出

    

    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;

    

    return 0;
}

附录2

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支持位数
complex<double> b[N]; // 用来充当临时调整空间的数组

void fft_merge_logic(complex<double> *a, int n, int inv){

    if(n == 1) return;

    fft_merge_logic(a,       n/2, inv);
    fft_merge_logic(a + n/2, n/2, inv);

    double unit_rad = 2 * pi / n;

    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i
        complex<double> tmp1 = a[i]; 
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }

}

void fft_rearrange_decidesequence_logic(int *rev, int n){

    // 给定数组rev和数组长度n，函数的功能是将所需的顺序写入数组，比如n = 4时将顺序 0 2 1 3 写入数组rev

    if(n == 1){
        rev[0] = 0;
        return;
    }

    // 获得 n/2 时的顺序,暂时放在rev的后半

    fft_rearrange_decidesequence_logic(rev + n/2, n/2);

    // 利用 n/2 时的顺序构造 n 时的顺序

    for(int i = 0; i < n/2; i ++){
        rev[i] = 2 * rev[i + n/2];
        rev[i + n/2] = 2 * rev[i + n/2] + 1;
    }

}

void fft_rearrange_logic(complex<double> *a, int n){

    // 按照算法，把a重新排列为可以直接递归向上的顺序

    // 计算bit: 满足pow(2, bit) = n

    int bit = 0;
    while((1 << bit) < n) bit ++;
    
    // rev: 确定a最终rearrange的位置序列

    int* rev = new int[n];
    fft_rearrange_decidesequence_logic(rev, n);

    // 按照rev序列调整a的顺序

    complex<double>* tmp_a = new complex<double>[n];

    for(int i = 0; i < n; i ++) tmp_a[i] = a[rev[i]];
    for(int i = 0; i < n; i ++) a[i] = tmp_a[i];

    delete[] rev;
}

void fft_logic(complex<double> *a, int n, int inv){

    // 参数：
    // 当inv = 1时，a是系数多项式，n是当前数组长度（2的幂次），函数效果是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当前数组长度（2的幂次），函数效果是原地变成系数多项式，但是所得的系数是n倍，需要在包装的函数中进行调整

    // a中元素的顺序调整为算法要求的顺序

    fft_rearrange_logic(a, n);
    
    // 调整顺序之后的a进行合并

    fft_merge_logic(a, n, inv);
}

void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}

void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}

// 主函数测试

complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘结果

int main(){

    // 两个相乘多项式的初始化。注意：两个相乘多项式和最后的乘积多项式需要有相同的项数

    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }

    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4

    // F，乘，F逆

    fft(A1, N);
    fft(A2, N);

    /*

    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;

    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;

    */

    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];

    ifft(C, N);

    // 输出

    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;

    return 0;
}

附录3

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
const int N = 1024; // 多项式的最大支持位数
complex<double> b[N]; // 用来充当临时调整空间的数组

void fft_merge_logic(complex<double> *a, int n, int inv){

    if(n == 1) return;

    fft_merge_logic(a,       n/2, inv);
    fft_merge_logic(a + n/2, n/2, inv);

    double unit_rad = 2 * pi / n;

    for(int i = 0; i < n/2; i ++){
        complex<double> x(cos(i * unit_rad), inv*sin(i * unit_rad)); // x = omega_n^i
        complex<double> tmp1 = a[i]; 
        complex<double> tmp2 = x * a[i + n/2];
        a[i]       = tmp1 + tmp2;
        a[i + n/2] = tmp1 - tmp2;
    }

}

void fft_rearrange_logic(complex<double> *a, int n){

    // 按照算法，把a重新排列为可以直接递归向上的顺序

    // 计算bit: 满足pow(2, bit) = n

    int bit = 0;
    while((1 << bit) < n) bit ++;
    
    // rev: 确定a最终rearrange的位置序列

    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }

    // 按照rev序列调整a的顺序

    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
    }

    delete[] rev;
}

void fft_logic(complex<double> *a, int n, int inv){

    // 参数：
    // 当inv = 1时，a是系数多项式，n是当前数组长度（2的幂次），函数效果是原地变成点值多项式
    // 当inv = -1时，a是点值多项式，n是当前数组长度（2的幂次），函数效果是原地变成系数多项式，但是所得的系数是n倍，需要在包装的函数中进行调整

    // a中元素的顺序调整为算法要求的顺序

    fft_rearrange_logic(a, n);
    
    // 调整顺序之后的a进行合并

    fft_merge_logic(a, n, inv);
}

void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}

void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}

// 主函数测试

complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘结果

int main(){

    // 两个相乘多项式的初始化。注意：两个相乘多项式和最后的乘积多项式需要有相同的项数

    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }

    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4

    // F，乘，F逆

    fft(A1, N);
    fft(A2, N);

    /*

    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;

    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;

    */

    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];

    ifft(C, N);

    // 输出

    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;

    return 0;
}

附录4

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
const int N = 8; // 多项式的最大支持位数
complex<double> b[N]; // 用来充当临时调整空间的数组

void fft_logic(complex<double> *a,int n,int inv){

    // bit: pow(2, bit) = n

    int bit = 0;
    while((1 << bit) < n) bit ++;
    
    // rev: 确定最终rearrange的位置序列

    int* rev = new int[n];
    for(int i = 0; i < n; i ++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }

    // 按照rev序列调整a的顺序

    for(int i = 0; i < n; i ++){
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
    }

    delete[] rev;
    
    // 调整顺序之后的a进行合并

    for (int mid=1;mid<n;mid*=2){
        // mid循环内，准备把两个长度为mid的序列合并成2 * mid的序列

        double unit_rad = 2 * pi / (2 * mid); // 单位角幅度值

        for (int i=0;i<n;i+=mid*2){
            // i循环内，把位置在i~i+mid*2位置的两个长度为mid的序列合并成2 * mid的序列

            for (int j = 0; j < mid; j ++){
                // j循环内的逻辑和
                complex<double> x(cos(j * unit_rad), inv*sin(j * unit_rad)); // x = omega_n^i
                complex<double> tmp1 = a[i+j]; 
                complex<double> tmp2 = x * a[i+j+mid];
                a[i+j]     = tmp1 + tmp2;
                a[i+j+mid] = tmp1 - tmp2;
            }
        }
    }
}

void fft(complex<double> *a, int n){
    // 输入系数多项式及其长度，原地转换为点值多项式
    fft_logic(a, n, 1);
}

void ifft(complex<double> *a, int n){
    // 输入点值多项式及其长度，原地转换为系数多项式
    fft_logic(a, n, -1);
    for(int i = 0; i < n; i ++) 
        a[i] /= n;
}

// 主函数测试

complex<double> A1[N], A2[N]; // 相乘的多项式
complex<double> C[N]; // 相乘结果

int main(){

    // 两个相乘多项式的初始化。注意：两个相乘多项式和最后的乘积多项式需要有相同的项数

    for(int i = 0; i < N; i ++){
        A1[i] = A2[i] = 0;
    }

    A1[0] = 1, A1[1] = 1, A1[2] = 3, A1[3] = 2; // A1(x) = 2x3 + 3x2 +  x + 1
    A2[0] = 4, A2[1] = 3, A2[2] = 5, A2[3] = 1; // A2(x) =  x3 + 5x2 + 3x + 4

    // F，乘，F逆

    fft(A1, N);
    fft(A2, N);

    /*

    for(int i = 0; i < N; i ++) 
        cout << A1[i] << " ";
    cout << endl;

    for(int i = 0; i < N; i ++) 
        cout << A2[i] << " ";
    cout << endl;

    */

    for(int i = 0; i < N; i ++) 
        C[i] = A1[i] * A2[i];

    ifft(C, N);

    // 输出

    

    for(int i = 0; i < N; i ++) 
        cout << C[i] << endl;

    

    return 0;
}

参考

  1. 十分简明易懂的FFT（快速傅里叶变换） ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

  2. 快速傅里叶变换（FFT）和逆快速傅里叶变换（IFFT） ↩︎

  3. C++ complex类 ↩︎

五种傅里叶变换的分析（FT、FS、DTFT、DFT、DFS）（基于MATLAB）
haohan的博客
05-30 3万+
五种傅里叶变换 FT： 傅里叶变换 Fourier Transform FS： 傅里叶级数 Fourier Series DTFT：离散时间傅里叶变换 Discrete-time Fourier Transform DFT： 离散傅里叶变换 Discrete Fourier Transform DFS： 离散傅里叶级数 Discrete Fourier Series 各种信号时域和频域的关系 时域 频域 连续 、非周期性 非周期性、连续 连续 、周期性 非周期性、离散 离散
快速傅里叶变换（FFT）学习笔记
Jun
07-27 605
Step1 模板理解: 1.背景：完成两个多项式的卷积 即：两多项式形如∑i=0naixi\sum_{i=0}^n a_{i}x^{i}∑i=0n​ai​xi ，∑i=0mbixi\sum_{i=0}^m b_{i}x^{i}∑i=0m​bi​xi 卷积后为 ∑i=0n+m(ci=∑j=0iajbi−j)xi\sum_{i=0}^{n+m}( c_{i}=\sum_{j=0}^ia_{j}b_{i...
[FFT] 快速傅里叶变换学习笔记
HeRaNO@HSAHRBNU
05-05 1万+
F(lag)F(FF)T(ransform)
hdu 4609（快速傅里叶变换）
北门的智慧
10-09 1593
题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4609； 题目大意：给出一些边长，问任取3边组成三角形的概率是多少。 分析：只要求出能组成三角形的有多少种再除以C(n,3)就好了。求得方法自然是用两最小边之和大于最长边。求和的部分可以用FFT优化。 代码如下： #include #include #include #include #
洛谷OJ P3803 【模板】多项式乘法（FFT）
baodream的博客
10-16 494
题目链接：https://www.luogu.org/problemnew/show/P3803 题目思路：FFT模板题 AC代码： // luogu-judger-enable-o2 #include &lt;stdio.h&gt; #include &lt;string.h&gt; #include &lt;iostream&gt; #include &lt;algorithm&gt;...
syzoj-hustoj:快速地搭建 SYZOJ
05-15
本仓库除了 syzoj 相关文件外均以 MIT 协议发布，syzoj 相关文件的许可证见压缩包内。 本教程仅在 Ubuntu 18.04 LTS 上测试通过，可能适用其它 Linux 发行版，但不支持 Windows 系统。如果使用 Windows，请使用...
OJ:OJ-在线法官，学习编码
05-06
**在线法官（Online Judge，OJ）是编程学习者和竞赛程序员的重要平台，它们提供了大量编程题目供用户练习和挑战。以下将详细讲解OJ、常见的OJ平台以及与C语言相关的编程知识点。** **一、OJ的概念** 在线法官（On...
基于Python的OJ刷题学习笔记设计源码
最新发布
10-02
该项目是一个基于Python的在线编程竞赛（OJ）刷题学习笔记设计源码，包含231个文件，包括202个Python源文件、10个PNG图片文件、9个Markdown文件、4个Python编译文件、2个文本文件、2个图片文件、1个Git忽略文件和1个...
oj：优化的JSON
01-31
一个快速的JSON解析器和Object marshaller作为Ruby gem。 3.0版已发布！ 3.0提供了更好的json gem和Rails兼容性。 它还提供了其他优化选项。 使用 require 'oj' h = { 'one' => 1 , 'array' => [ true , false ] }...
快速傅里叶变换fft
hudada19的博客
02-16 340
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如 果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱 提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去 做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在圈圈就根...
[UOJ 34]多项式乘法(FFT)
某沙茶的OI代码库
03-12 1500
题目链接http://uoj.ac/problem/34思路裸的FFT模板题。 具体的FFT算法请看我的另一篇闲扯http://blog.csdn.net/qpswwww/article/details/44223331代码#include <stdio.h> #include <stdlib.h> #include <string.h> #include <algorithm> #include
BNUOJ 27935 我爱背单词（FFT）
AC_Arthur的专栏
01-12 1064
题目链接：点击打开链接 思路： 该题暴力当然可以过，   如果数据量加大，  我们还有一种nlogn的算法：FFT 仔细观察这个复习单词量的累加方式可以发现， 这是一个卷积， 可以用FFT加速算法。 细节参见代码： #include #include #include #include #include #include #include #include #includ
FFT理解
ACM之路
05-30 215
https://www.zhihu.com/question/22611929
各种傅里叶变换
likejoey的专栏
03-17 761
傅里叶变换是基本，后续各种派生都是处于某个目的。Fourier 名称时域形式频域形式FT (Fourier Transform)连续非周期连续非周期FS (Fourier Series)连续周期离散非周期DTFT(Discrete Time Fourier Transform)离散非周期连续周期DFT(Disxrete Fourier Transform)离散周期离散周期   transform ...
天大OJ:卡牌大师C++代码
12-30
题目描述： 分配抽取N张扑克牌，分别从A-K,其中规定A为1点；J,Q,K分别为11,12,13点；数字牌点数与其所表示的数字一致。现在可通过移动N张牌的任意次序，并在两两之间插入预算符号+，-，*，/，（注：除法向下取整）。...
常用函数的傅里叶变换汇总
热门推荐
Varalpha的博客
03-19 15万+
快速链接: 【总目录】信号与系统 学习笔记 Signals and Systems with Python (1) 简介 Intro (2) 傅里叶 Fourier 常用函数的傅里叶变换汇总 (3) LTI 系统 与 滤波器 常用函数的傅里叶变换汇总 f(t)⟵⟶F(jω)F(jt)⟵⟶2πf(−ω)f(αt)⟵⟶1∣α∣F(jωα)a⋅f1+b⋅f2⟵⟶a⋅F1+b⋅F2f(t±...
FFT小总结
weixin_30443731的博客
09-16 204
FFT实质上做的是循环卷积，ck=sigam（ai*bj,(i+j)%n=k）,其中n是倍长后的长度，所以我们有时候需要的只是普通的卷积，我们就需要把原数组倍长，再用FFT求卷积，由于高位都是0，所以做出来就等价于普通的卷积。 转载于:https://www.cnblogs.com/OYzx/p/5877649.html...
FFT（dddd）
DS00HY的小博客
10-04 251
  首先推荐几个博客 其实应该先看一看卷积：https://blog.csdn.net/qq_39521554/article/details/79083864 关于FFT的https://www.cnblogs.com/fenghaoran/p/7107608.html https://blog.csdn.net/qq_39521554/article/details/79083864 ...
【原创】FFT（快速傅里叶变换）学习笔记 未完成
致上
11-17 324
为你撰写的故事FFT说在前面从多项式说起多项式的表示法系数表示法点值表示法多项式的运算法则及朴素算法时间复杂度分析 FFT 说在前面 写了一上午，然后浏览器卡了，触发了自动刷新功能，87%的都没有了。 然后一怒之下删掉了。 从多项式说起 多项式的表示法 系数表示法 f(x)⁡=a0x0+a1x1+a2x2+⋯+anxn \large\operatorname{f(x)}=a_0x^0+a_1x^1...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值