zoukankan      html  css  js  c++  java
  • 模板:快速傅里叶变换(FFT)

    参考:http://blog.csdn.net/f_zyj/article/details/76037583
    如果公式炸了请去我的csdn博客:http://blog.csdn.net/luyouqi233/article/details/79323568
    原文即是一篇很好的FFT入门博客,但是笔者打算为了日后的学习,则将原篇章的结构删改增添一下,如有思路上的雷同十分正常。
    “是时候打开FFT的大门了!”

    预备知识:

    1.至少知道基础数论与一定解三角形知识(大概是高中水平)。
    2.定义(i=sqrt{-1})
    3.引入复数(即形如(a+bi)(a,b均为实数)的数的集合)
    4.((cos heta+i imes sin heta)^k=cos(k heta)+i imes sin(k heta))
    5.显然我们对多项式FFT之后得到的答案不是我们想要的,那么这时候就需要反着用FFT把式子再变回去(本文记做IFFT)。

    这里证明一下第四条,用归纳法。
    显然当(k=1)时成立。
    (k)成立时,我们有:
    ((cos heta+i imes sin heta)^{k+1})
    (=(cos heta+i imes sin heta)^k imes (cos heta+i imes sin heta))
    (=(cos(k heta)+i imes sin(k heta)) imes (cos heta+i imes sin heta))
    (=cos(k heta)cos heta+i imes sin(k heta)cos heta+i imes cos(k heta)sin heta+i^2 imes sin(k heta) sin heta)
    (=cos(k heta)cos heta-sin(k heta) sin heta+i imes (sin(k heta)cos heta+cos(k heta)sin heta))
    (=cos((k+1) heta)+i imes sin((k+1) heta))
    得证。

    问题引入:

    (A(x)=sum_{i=0}^{n-1}a_ix^i,B(x)=sum_{i=0}^{n-1}b_ix^i),求(A(x) imes B(x))后的多项式系数。

    初探:

    显然我们有一个(O(n^2))的解法,但是实在是太慢了。
    考虑到一个(n-1)次多项式可以看做是定义在复数域上的函数,则我们一定可以找到n个点来唯一确定这个函数。
    当然我们也可以通过这些点来表示这个多项式。
    假设:
    (A(x))被表示为:(<(x_0,y_{a_0}),(x_1,y_{a_1}),ldots,(x_{2n-2},y_{a_{2n-2}})>)
    (B(x))被表示为:(<(x_0,y_{b_0}),(x_1,y_{b_1}),ldots,(x_{2n-2},y_{b_{2n-2}})>)
    显然(A(x) imes B(x))被表示为:(<(x_0,y_{a_0}y_{b_0}),(x_1,y_{a_1}y_{b_1}),ldots,(x_{2n-2},y_{a_{2n-2}}y_{b_{2n-2}})>)

    这里多取了点的原因在于(A(x) imes B(x))是一个(2n-2)次多项式,则至少要取(2n-1)个点才能保证正确。

    但是显然还是(O(n^2))的。

    再试:

    考虑设(A(x_i)=A_0(x_i^2)+x_iA_1(x_i^2)),其中:
    (A_0(x)=a_0+a_2x+a_4x^2+ldots+a_{n-2}x^{frac{n}{2}+1})
    (A_1(x)=a_1+a_3x+a_5x^2+ldots+a_{n-1}x^{frac{n}{2}+1})

    其实就是按照系数下标的奇偶性分类了一下。

    此时我们再令取点的(x)值为(<x_0,x_1,ldots,x_{frac{n}{2}-1},-x_0,-x_1,ldots,-x_{frac{n}{2}-1}>)
    我们发现把(x)平方后我们的取值瞬间缩小了一半,而原式唯一变化的就是(A_1(x))前的符号。
    看起来我们似乎找到了(O(nlogn))的可行方案。
    但是很可惜,这样优秀的(x)取值的性质只会保留一次,也就是说我们只是得到了一个(O(frac{n^2}{2}))
    如何才能每次将问题的规模缩小一半是我们的目标。

    插曲:

    有个人告诉你:不如试试(X_n=cosfrac{2pi}{n}+i imes sinfrac{2pi}{n})(0ldots n-1)次方作为(x)的取值。

    这块大家一直有个疑惑:这是怎么构造出来的啊?
    事实上傅里叶变换最早是应用于信号处理上的,傅里叶提出:任何连续周期信号可以由一组适当的正弦曲线组合而成。
    多项式可以看做非连续周期信号,然后通过各种奇妙的姿势让它逼近正弦曲线的组合形,详情可以看松松松WC2018的课件。
    “逼近”显然用到了微积分,不适合初学者,所以就直接跳过了。(其实我也不会……)
    (再多说一点吧,其实上面和下面的数学推理完全可以从物理层面理解,还是可以参考松松松WC2018的课件)

    继续:

    那么令取点的(x)值为(<X_n^0,X_n^1,ldots,X_n^{n-1}>)
    我们可知:
    ((X_n^{k})^2)

    (=X_n^{2k})

    (=cosfrac{2k imes 2pi}{n}+i imes sinfrac{2k imes 2pi}{n})

    (=cosfrac{2kpi}{frac{n}{2}}+i imes sinfrac{2kpi}{frac{n}{2}})

    (=X_{frac{n}{2}}^k)


    (X_n^{k})

    (=cosfrac{k imes 2pi}{n}+i imes sinfrac{k imes 2pi}{n})

    根据三角函数的周期性可知,(k)(n)取模显然不会对答案造成影响。
    于是我们有(X_n^{k}=X_n^{k\%n})

    那么显然对于(<(X_n^0)^2,(X_n^1)^2,ldots,(X_n^{n-1})^2>)

    它等效于(<X_{frac{n}{2}}^0,X_{frac{n}{2}}^1,ldots,X_{frac{n}{2}}^{frac{n}{2}-1},X_{frac{n}{2}}^0,X_{frac{n}{2}}^1,ldots,X_{frac{n}{2}}^{frac{n}{2}-1}>)

    我们好像看到了(O(nlogn))的曙光了。

    尾声:

    显然我们可以对(x)的取值折半,然后对于左右区间的(x)值递归下去即可。
    Q1:诶等等,“再试”里面的内容好像没有应用上啊……

    A1:那就转化一下,其实我们只需要求一个区间的(A_0(x))(A_1(x))值递归下去求(A(x))即可。
    也就是说其实我们是得到了:
    (<(A_0)_0,(A_0)_1,ldots,(A_0)_{frac{n}{2}-1},(A_1)_0,(A_1)_1,ldots,(A_1)_{frac{n}{2}-1}>)

    Q2:这好像是画蛇添足……

    A2:emmm……我说这个可以用于常数优化你信吗……
    显然(A(X_n^k)=(A_0)_{k\%frac{n}{2}}+X_n^k(A_1)_{k\%frac{n}{2}})

    取模是因为,不要忘了我们的取值是由两个一样的左右区间合并在一起的。

    那么我们得到了(<A_0,A_1,ldots,A_{n-1}>)

    (其中(A_k=A(X_n^k))

    我们好像把这个序列的长度减少了一半诶!那自然是快了二倍啊。

    不要忘了n要满足始终是2的倍数,所以n要取2的整数次幂,同时将没用的次幂的系数填成0。

    Q3:IFFT怎么做啊?

    A3:继续看下去……?

    补遗:

    略讲一下IFFT。
    显然我们可以把FFT的最初算法(也就是DFT)看做两个矩阵相乘。

    两个矩阵分别一个填((X_n^k)^m),一个填系数,可以上参考处原博客看矩阵。

    那么我们把第一个矩阵变成逆矩阵岂不是为IFFT?
    其实就是这样,并且事实上就是填(((X_n^{-k})^m)/n),具体证明过程看参考处原博客。
    剩下的做法就和FFT一样啦。

    谢幕:

    事实上我上述讲的内容其实没有多少用(滑稽。
    因为你理解半天也不如不理解知道怎么用然后默写下来。
    但是理解了更好背啊。

    例题:

    模板:
    HDU1402:A * B Problem Plus:http://www.cnblogs.com/luyouqi233/p/8448969.html

    应用:
    BZOJ3527:[ZJOI2014]力:http://www.cnblogs.com/luyouqi233/p/8452117.html

    +++++++++++++++++++++++++++++++++++++++++++

    +本文作者:luyouqi233。               +

    +欢迎访问我的博客:http://www.cnblogs.com/luyouqi233/+

    +++++++++++++++++++++++++++++++++++++++++++

  • 相关阅读:
    mitmproxy抓包工具
    java基础|int和Integer的区别
    Vue|退出功能
    Vue|分页处理
    apt-get本地软件源搭建
    rqt_plot报错
    创建ROS 工作空间时出现:程序“catkin_init_workspace”尚未安装,程序“catkin_make”尚未安装。
    ubuntu16.04安装ROS
    debian及Ubuntu各版本下载地址获取
    解决sudo rosdep init和rosdep update的错误
  • 原文地址:https://www.cnblogs.com/luyouqi233/p/8447569.html
Copyright © 2011-2022 走看看