快速傅里叶变换FFT简明教程

87 阅读 0 评论 58 点赞

我是靠谱客的博主美满烧鹅，最近开发中收集的这篇文章主要介绍快速傅里叶变换FFT简明教程，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

问题

给出长度为 $n$ 的多项式 $f(x)=sum_{i=0}^{n-1}a_{i}x^{i}$ 和长度为 $m$ 多项式 $g(x)=sum_{i=0}^{m-1}b_{i}x^{i}$ 。求解多项式 $h (x) = f (x) \times g (x)$ 。
$10^{5})$

系数表示法和点值表示法

多项式的两种表示方法，系数表示法和点值表示法。
系数表示法：我们存储 $f (x)$ 的每一项 $x^{i}$ 的系数 $a_{i}$ 。
点值表示法：由大家都知道知， $N$ 个点可以确定一个最高次为 $N - 1$ 多项式，所以我们可以存储多项式上任意 $N$ 个点，来确定这个多项式。

点值表示法有个性质是，对于多项式 $f$ 的一个点值 $x,y_{1})$ ，和多项式 $g$ 的一个点值 $x,y_{2})$ ，可以得到两个多项式相乘的多项式 $h$ 在横坐标为 $x$ 处的点为 $(x,y_{1}times y_{2})$ 。

所以在对于长度为 $n$ 的多项式 $f$ 和长度为 $m$ 的多项式 $g$ 做乘法时，我们可以算出 $f$ 的 $n + m - 1$ 个点值和 $g$ 的 $n + m - 1$ 个点值，且要求两个多项式点值所取的 $x$ 相同。
我们就可以通过遍历算出这 $n + m - 1$ 个横坐标的点值纵坐标相乘，来得到相乘后多项式的 $n + m - 1$ 个点值，也就是新多项式的点值表示法。
(注意：虽然我们存储 $f$ 只需要 $n$ 个点值，但是为了计算出长度为 $n + m - 1$ 的 $h$ ，我们仍要先得到 $f$ 的(n+m-1)个点值)

因此，假如我们可以通过 $o (n l o g n)$ 的办法，将系数表示法转化为点值表示法，然后将两个多项式点值 $o (n)$ 相乘，最后通过 $o (n l o g n)$ 的办法将点值表示法转化回系数表示法，就可以解决这个问题了。

后面来介绍如何将系数表示发和点值表示法进行转换。

离散傅里叶变换（DFT）

对于没有周期的函数，我们可以认为其周期为 $\infty$ ，用无穷个周期函数去拟合。
同理，对于一个长度为 $N$ 的离散多项式 $x (n) (n = 0, 1, 2 . . N - 1)$ ，可以通过 $N$ 个周期函数来拟合。（ $N$ 个周期函数只能拟合 $N$ 个点，不能用来拟合连续多项式函数）
由此我们得到离散傅里叶级数展开式如下。
$X(k)=sum_{n=0}^{N-1}x(n)e^{-2pi ikfrac{n}{N}}$
其中 $e^{-2pi ikfrac{n}{N}}=cos(-2pi kfrac{n}{N})+isin(-2pi kfrac{n}{N})$ 是关于 $n$ 的周期为 $N$ 的函数。（由欧拉公式 $e^{ix}=cos x+isin x$ 得）

离散傅里叶逆变换（IDFT）

对于变换后得到的 $X (n)$ ，我们可以通过离散傅里叶逆变换公式得到原函数 $x (n)$ 。与 $D F T$ 公式差不多，在原来的基础上 $i$ 变成了 $- i$ ，前面乘上一个 $1 N frac{1}{N}$ 。公式如下。
$x(n)=frac{1}{N}sum_{k=0}^{N-1}X(k)e^{2pi ikfrac{n}{N}}$

快速傅里叶变换（FFT）

由前面讲过的，我们要解决的核心问题就是，算出多项式的 $N$ 个点值。
大致做法是，我们求出函数在离散傅里叶变换后的N个点值，由点值乘出新多项式后，再将其逆变换回系数表示法。

我们先来观察式子 $D F T$ 的式子。
我们为了简化式子，我们把 $e^{frac{2pi ik}{N}}$ 用 $W_{N}$ 代替。( $W_{N}$ 是一个关于 $k$ 的函数)
得到式子如下。
$X(k)=sum_{n=0}^{N-1}x(n)W_{N}^{n}(k)$
然后观察 $W_{N}$ 如下。
$W_{N}(k)=e^{frac{2pi i}{N}}=cos(frac{2pi k}{N})+isin(frac{2pi k}{N})$
发现 $W_{N}$ 是一个周期为 $N$ 的函数。

由此，我们发现变换后的函数 $X (k)$ 是由 $N$ 个周期为 $N$ 的函数的幂次，乘上一个系数，相加而成。且这 $N$ 个函数的实部和虚部都是三角函数（ $c o s$ 和 $s i n$ ）。

我们很轻易地发现了 $W_{N}(k)=W_{N}(k+N)$ ，但是这对我们来说不重要。更重要的性质源于三角函数的另一个性质——对称性。
我们把 $[0, 2 π]$ 分成两部分 $[0, π)$ ，和 $(π, 2 π]$ ，可以发现我们算出左边 $t$ 处的值，就可以得到右边 $π + t$ 处的值。
说白了就是 $s i n (t) = - s i n (π + t)$ 和 $c o s (t) = - c o s (π + t)$ 。

对应到上面就得到了一个结论，我们只要算出了 $W_{N}(t)$ ，就可以得到 $W_{N}(frac{N}{2}+t)$ 。即 $W_{N}(t)=-W_{N}(frac{N}{2}+t)$ 。

上面只是简单的性质推导，正式的傅里叶变换在下面。

总而言之就是要分治，这个式子可以和分治有机结合。（我也不知道他怎么想出来的，公式推导在《数字信号处理》第四版的P122）
我们把枚举的奇数项和偶数项分离开。
$X(k)=sum_{n为偶数}x(n)W_{N}^{n}(k)+sum_{n为奇数}x(n)W_{N}^{n}(k) \=sum_{r=0}^{frac{N}{2}-1}x(2r)W_{N}^{2r}(k)+sum_{r=0}^{frac{N}{2}-1}x(2r+1)W_{N}^{2r+1}(k) \=sum_{r=0}^{frac{N}{2}-1}x_{1}(r)W_{N}^{2r}(k)+W_{N}(k)sum_{r=0}^{frac{N}{2}-1}x_{2}(r)W_{N}^{2r}(k) \=sum_{r=0}^{frac{N}{2}-1}x_{1}(r)W_{frac{N}{2}}^{r}(k)+W_{N}(k)sum_{r=0}^{frac{N}{2}-1}x_{2}(r)W_{frac{N}{2}}^{r}(k) \=X_{1}(k)+W_{N}(k)X_{2}(k)$
其中 $X_{1}$ 是 $x_{1}$ 的傅里叶变换， $X_{2}$ 是 $x_{2}$ 的傅里叶变换。
其中有一步 $W_{N}^{2r}(k)=W_{frac{N}{2}}^{r}$ 前面忘记证了，想必读者应该一下子就能证出来了。

到此位置，所有准备工作都做完了，容易发现只要分治做好两个小多项式的点值，大的自然就出来了。所以可以分治做。

最后式子如下。
$X(k)=X_{1}(k)+W^{k}_{N}X_{2}(k),k=0,1,…,frac{N}{2}-1 \X(k+frac{N}{2})=X_{1}(k)-W^{k}_{N}X_{2}(k),k=0,1,…,frac{N}{2}-1$

但是有一个细节问题，计算 $W_{N}(k)$ 的复杂度是 $o (l o g n)$ 。而在分治的每一层我们需要求 $N$ 次 $W$ 函数，总共分治 $l o g n$ 层，最终导致我们的复杂度变成了 $o (n l o g n l o g n)$ 。
那么解决办法是前面提的 $W_{N}(t)=-W_{N}(frac{N}{2}+t)$ ，我们每层只需要计算一半的 $W$ ，也就是 $N 2 frac{N}{2}$ 个 $W$ ，于是在计算 $W$ 上的总开销就变成了 $o ( n 2 l o g n l o g n ) o(frac{n}{2}lognlogn)$ 。有锅巴用。

然后又发现一个可以优化的地方。既然 $W_{N}(k)=-W_{N}(frac{N}{2}+k)$ 来减少 $W$ 的计算次数（主要是减少 $s i n$ 和 $c o s$ 的计算次数），主要是来源于三角函数知道 $[0, π]$ 可以推演出 $[π, 2 π]$ 。这样只需要在每层算 $N 2 frac{N}{2}$ 次 $c o s$ 和 $s i n$ ，所以我们再往下细化一层，使得每层只需要计算 $N 4 frac{N}{4}$ 次 $c o s$ 和 $s i n$ 。
那么同理，三角函数知道 $[ 0 , π 2 ] [0,frac{pi}{2}]$ 就可以推算出 $[ π 2 , 2 π ] [frac{pi}{2},2pi]$ 。
所以，我们开始推 $W_{N}(k+frac{N}{4})$ 。
$W_{N}(k+frac{N}{4})=cos(frac{2pi k}{N}+frac{pi}{2})+isin(frac{2pi k}{N}+frac{pi}{2})=-sin(frac{2pi k}{N})+icos(frac{2pi k}{N})$
然后，我们发现，我们对于某个 $k$ 计算出的 $s i n$ 和 $c o s$ 函数，可以沿用至 $k + N 4 k+frac{N}{4}$ ， $k + N 2 k+frac{N}{2}$ ， $k + 3 N 4 k+frac{3N}{4}$ 。
所以，在每一层，我们只需要计算 $N 4 frac{N}{4}$ 次 $W$ 。
然后复杂度被优化到了 $o ( n 4 l o g n l o g n ) o(frac{n}{4}lognlogn)$ 。仍然没有任何卵用。

其实我们只需要算 $W_{N}(1)$ ，有 $W_{N}(k+1)=W_{N}(k)times W_{N}(1)$ ，这才是真的有用。证明可以考虑这个复数的几何意义，就是把位于圆心的单位圆划等分成了 $N$ 个部分，取出第一个部分就相当于是一个旋转矩阵，幂次就是旋转的次数。
于是，我们只要在开始用三角函数算出 $W_{N}(1)$ ，后面一直乘上它就行了。

然后就做完了。

关于转换回系数，也就是逆变换，发现式子和这个差不多，所以做法也基本一样。记得最后要除以 $N$ 。这里就不重复介绍了。

为了方便起见，一般先把多项式长度拉到2的幂次，然后再进行变换，方便递归到底层为1。反过来说，要使用FFT板子时，注意长度应为2的幂次。
虽然是分治，但为了常数考虑，一般都用非递归版的写法。

递归版
递归版常数很大，强烈建议找个常数小的非递归版代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const double pi = acos(-1);
const int MAXN = 3000005;
struct cp
{
    double s,v;
    cp(double ss=0,double vv=0){
        s=ss;v=vv;
    }
};
cp operator * (cp a,cp b)
{
    return cp(a.s*b.s-a.v*b.v,a.s*b.v+a.v*b.s);
}
cp operator + (cp a,cp b)
{
    return cp(a.s+b.s,a.v+b.v);
}
cp operator - (cp a,cp b)
{
    return cp(a.s-b.s,a.v-b.v);
}
cp operator / (cp a,int b)
{
    return cp(a.s/b,a.v/b);
}

//预处理出2的幂次的cos和sin
double coss[MAXN],sinn[MAXN];

void fft(cp *a,int len,int opt)
{
    if(len==1)return ;
    static cp b[MAXN];
    int mid=len>>1,mid_2=mid>>1;
    for(int i=0;i<len;i++)b[(i>>1)|(i&1?mid:0)]=a[i];
    for(int i=0;i<len;i++)a[i]=b[i];
    fft(a,mid,opt);fft(a+mid,mid,opt);
    cp uni(coss[mid],opt*sinn[mid]),w(1,0);
    for(int k=0;k<len;k++)
    {
        if(k<mid)b[k]=a[k]+w*a[k+mid];
        else b[k]=a[k-mid]+w*a[k];
        w=w*uni;
    }
    for(int i=0;i<len;i++)a[i]=b[i];
}
void mul(cp *a,cp *b,int len)
{
    fft(a,len,1);
    fft(b,len,1);
    for(int i=0;i<len;i++)a[i]=a[i]*b[i];
    fft(a,len,-1);
    for(int i=0;i<len;i++)a[i]=a[i]/len;
}
string sa,sb,sans;
cp a[MAXN],b[MAXN];
int ans[MAXN];
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin>>sa>>sb;
    int n=sa.size(),m=sb.size();
    int len=1;
    while(len<n+m-1)coss[len]=cos(pi/len),sinn[len]=sin(pi/len),len<<=1;
    for(int i=0;i<len;i++)
    {
        a[i]=i<n?cp(sa[n-i-1]-'0',0):cp(0,0);
        b[i]=i<m?cp(sb[m-i-1]-'0',0):cp(0,0);
    }
    
    mul(a,b,len);
    
    for(int i=0;i<len;i++)
    {
        ans[i]+=round(a[i].s);
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    int flag=0;
    for(int i=len-1;i>=0;i--)
    {
        if(ans[i]==0&&flag==0)continue;
        else
        {
            flag=1;
            sans+=(char)(ans[i]+'0');
        }
    }
    cout<<sans<<endl;
    return 0;
}