FFT（快速傅里叶变换）——从入门到放弃

· · 个人记录

FFT（Fast Fourier Transformation）就是“快速傅里叶变换”的意思，它是一种用来计算DFT（离散傅里叶变换）和IDFT（离散傅里叶逆变换）的一种快速算法。这种算法运用了一种高深的数学方式、把原来复杂度为O(n^2)的朴素多项式乘法转化为了O(nlogn)的算法。

最常用的FFT算法是Cooley和Turkey提出的Radix-2 FFT（基-2FFT），把一个DFT拆成两个一半大小的DFT。这篇文章只介绍Radix-2 FFT。此方法因为实现简单倍受欢迎。Radix-2 FFT适用于N为2的整数次方的FFT。

此外还有Radix-4 FFT（基-4FFT），Radix-8 FFT（咳咳FFT），它们效率更高，但适用范围太窄。

综合了Radix-2的适用范围和Radix-4的效率的算法是SplitRadix FFT（分裂基FFT），名字听上去就很帅，实际玩起来也挺帅的……效率比Radix2快一点吧大概，但实现起来复杂很多。

一、朴素的多项式乘法

设f(x)和g(x)各表示一个n-1次多项式即：

f(x)=\sum^{n-1}_{i=0}a_i*x^i g(x)=\sum^{n-1}_{i=0}b_i*x^i

将所有的两项系数相乘后再相加，利用这种方法计算多项式乘法复杂度为O(n^2)。对于1e6的数据范围来说，它的时间复杂度是无法接受的。

二、多项式的点值表示

将n个互不相同的x带入多项式，会得到n个不同的取值y。则该多项式被这n个点(x1,y1),(x2,y2),…,(xn,yn)唯一确定。其中

y_i= \sum ^{n-1}_{j=0}a_j*x^j_i

已知f(x)与g(x)的点值表示，可以在O(n)的时内求出f(x)*g(x)的点值表示。

但是对于系数表示的多项式，在转换为点值表示时，由于要选n个点，每个点O(n)，这种表示法计算多项式乘法的时间复杂度仍然为O(n^2)，而如果利用高斯消元转换回系数表示法，所需的时间还要更多。

可以优化吗？

三、复数

可以借助复数优化。

i的定义

### 复数的表示 形如$a+bi$的数叫复数。在复平面中，$x$轴代表实数，$y$轴（除原点外的点）代表虚数，从原点$(0,0)$到$(a,b)$的向量表示复数$a+bi$。这个向量的模长就是这个复数的模长。以逆时针为正方向，从$x$轴正半轴到已知向量的转角的有向角叫做幅角。 ``` struct complex { double x,y; complex (double a=0,double b=0){x=a,y=b;} }; ``` ### 复数的运算 #### 加法、减法 与向量相同。 ``` complex operator + (complex a,complex b) { return complex(a.x+b.x , a.y+b.y); } complex operator - (complex a,complex b) { return complex(a.x-b.x , a.y-b.y); } ``` #### 乘法 由定义知，两个复数$a+bi$与$c+di$的乘积： $$(a+bi)*(c+di)$$ $$=ac+adi+bci-bd$$ $$=(ac-bd)+(bc+ad)i$$ 所以： ``` complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); } ``` ### 单位根 #### 单位根的定义 在复平面上，以原点为圆心，$1$为半径作圆，所得的圆叫单位圆。可以发现，单位圆上的两个复数相乘，结果仍在单位圆上，且幅角等于原来两个复数的幅角之和。 证明：设复数$cos\theta+sin\theta i$与$cos\lambda+sin\lambda i (cos\theta+sin\theta i)*(cos\lambda+sin\lambda i) =cos\theta cos\lambda-sin\theta sin\lambda+cos\theta sin\lambda i+sin\theta cos\lambda i =cos(\theta+\lambda)+sin(\theta+\lambda)i

单位根的求法

所以，以原点为起点，圆的n等分点为终点，做n个向量，设幅角为正且最小的向量对应的复数为\omega _n，称为n次单位根。显然，其余n-1个复数为\omega^2_n,\omega^3_n,…,\omega^n_n(\omega^0_n=\omega^n_n=1)。

对三角函数稍有常识的人就能看出，\omega^k_n=cos\frac{2\pi*k}{n}+sin\frac{2\pi*k}{n}i。

单位根的性质

1.\omega^{2k}_{2n}=\omega^k_n

证明：

\omega^{2k}_{2n}=cos\frac{2\pi*2k}{2n}+sin\frac{2\pi*2k}{2n}i =\omega^k_n=cos\frac{2\pi*k}{n}+sin\frac{2\pi*k}{n}i=\omega^k_n

2.\omega^{k+\frac{n}{2}}_n=-\omega^k_n

证明：

\omega^{k+\frac{n}{2}}_n=cos\frac{2\pi*(2k+n)}{2n} =cos(\frac{2\pi*2k}{2n}+\pi)=-\omega^k_n

四、DFT（离散傅里叶变换）

我们可以把n次单位根代入两个多项式，获得两个多项式的点值表示。

且慢，一共n个单位根，每个单位根代入n次，这岂不是O(n^2)的时间复杂度吗？

其实是可以优化的：

f(x)=a_0+a_1*x+a_2*x^2+......+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}

按奇偶分组：

f(x)=(a_0+a_2*x^2+......+a_{n-2}*x^{n-2})+(a_1*x+a_3*x^3+......+a_{n-1}*x^{n-1})

设：

f_0(x)=a_0+a_2*x+......+a_{n-2}*x^{\frac{n-2}{2}} f_1(x)=a_1+a_3*x+......+a_{n-1}*x^{\frac{n-2}{2}}

显然，

f(x)=f_0(x^2)+xf_1(x^2)

将\omega^k_n代入：

f(\omega^k_n)=f_0(\omega^{2k}_n)+\omega^k_nf_1(\omega^{2k}_n)

再将\omega^{k+\frac{n}{2}}_n代入：

f(\omega^{k+\frac{n}{2}}_n)=f_0(\omega^{2k+n}_n)+\omega^{k+\frac{n}{2}}_nf_1(\omega^{2k+n}_n) =f_0(\omega^{2k}_n)-\omega^k_nf_1(\omega^{2k}_n)

两者实部相同，只有虚部相反。因此在算第一个时可以顺便O(1)算出第二个。问题就这样缩小了一半。而由于缩小后的问题仍能这样解决，我们就可以递归下去。这样我们就获得了一个O(nlogn)的做法。需要注意的是，我们要把多项式的次数补成2^n-1次，使得分治可以正常进行。

可以改成非递归版吗？

观察分组过程：

原数组 0 1 2 3 4 5 6 7
分组 0 2 4 6 1 3 5 7
再分组 0 4 2 6 1 5 3 7

转换成二进制：

000 001 010 011 100 101 110 111
000 010 100 110 001 011 101 111
000 100 010 110 001 101 011 111

在这个过程结束后，下标的二进制恰好与原来相反！

那么我们就可以用下面的方法O(n)求出所有反转后的二进制下标了：

for(int i=0;i<limit;i++)
        r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ;

这也就是所谓的蝴蝶运算。

原理：

我们可以把一个二进制数看成两部分，它的前bit-1位是一部分，它的最后一位是一部分。

一个数的二进制翻转就相当于是把它的最后一位当成首位，然后在后面接上它前bit-1为的二进制翻转。而且在这个循环中我们能保证，在计算“i”的二进制翻转之前1~i-1中的所有数的二进制翻转都已经完成。“i”的前bit-1位的数值其实就是i>>1的值，直接调用i>>1的二进制翻转的结果就相当于调用了“i”的前bit-1位二进制翻转的结果。

至此，我们已经求出了两个多项式的点值表示。然后即可在O(n)时间内求出两者乘积的点值表示。

五、IDFT（离散傅里叶逆变换）

可是，我们要的是系数啊！

所以我们要考虑如何把点值表示法转换为系数表示法，这个过程就是IDFT。利用单位根的一些性质，可以加速这个转换的过程。

先给出结论：

a_k=\frac{\sum^{n-1}_{i=0}f(\omega^i_n)*(\omega^{-k}_n)^i}{n}

其中a_k是k次项系数，f(x)是多项式在x处的点值表示。

我们试着来证明一下：考虑式子右上边的

\sum^{n-1}_{i=0}f(\omega^i_n)*(\omega^{-k}_n)^i

将f(x)展开：原式=

\sum^{n-1}_{i=0}(\sum^{n-1}_{j=0}a_j(\omega^i_n)^j)(\omega^{-k}_{n})^i =\sum^{n-1}_{i=0}(\sum^{n-1}_{j=0}a_j(\omega^j_n)^i)(\omega^{-k}_{n})^i

显然(\sum^{n-1}_{j=0}a_j(\omega^j_n)^i)(\omega^{-k}_{n})^i=\sum^{n-1}_{j=0}(a_j(\omega^j_n)^i(\omega^{-k}_{n})^i)。

所以原式=

\sum^{n-1}_{i=0}\sum^{n-1}_{j=0}(a_j(\omega^j_n)^i(\omega^{-k}_{n})^i) =\sum^{n-1}_{i=0}\sum^{n-1}_{j=0}(a_j(\omega^{j-k}_n)^i) =\sum^{n-1}_{j=0}(\sum^{n-1}_{i=0}(\omega^{j-k}_n)^i)a_j

设S(x)=\sum^{n-1}_{i=0}x_i

首项为1，公比\omega^k_n，由等比数列求和公式得：

S(\omega^k_n)=\frac{(1-(\omega^k_n)^n)}{1-\omega^k_n}

显然(\omega^k_n)^n=(\omega_n^1)^{nk}=(\omega^n_n)^k=1，于是分子=1-1=0

我们惊讶地发现：S(\omega^k_n)=0！

但是别忘了：我们没有考虑k=0时的情况。此时的每一项都是1，因此和为n。

将这个结果回带到\sum^{n-1}_{j=0}(\sum^{n-1}_{i=0}(\omega^{j-k}_n)^i)a_j里面：

只有j=k时\sum^{n-1}_{i=0}(\omega^{j-k}_n)^i值为n，否则值为0。也就是说原式=na_k。

因此：

a_k=\frac{\sum^{n-1}_{i=0}f(\omega^i_n)*(\omega^{-k}_n)^i}{n}

只要把多项式的点值表示当做系数，代入相应单位根的倒数，用DFT的方法求出结果，所得的结果再除以n，就可以得到多项式的系数表示！

如何求出单位根的倒数呢？只要找到另一个复数，让两个复数模长相同，幅角相加等于2\pi即可。画图可知，若单位根\omega^k_n=cos\theta+sin\theta i，那么它的倒数就等于cos\theta-sin\theta i。

六、总结

七、扩展：傅里叶变换应用初探

MP3作为一种有损压缩格式，通过剔除掉音乐中特定人耳听觉不敏感的频段，减少数据量，获得更大的压缩比。这个“剔除”的过程就用了DFT。

基于Vocaloid歌声合成软件的初音未来，其合成过程用到了大量的DFT。Vocaloid的合成原理基于STFT(短时傅立叶变换)，单音轨合成，每秒至少要进行344次2048点FFT，从而修改信号的频谱特征，Vocaloid合成时大部分时间花在FFT上。所以，你的老婆就是DFT！

汤姆猫之类的变声器使用DFT改变声音的频谱特征。

均衡器，它通过DFT增强或削弱声音中某个特定频率范围的强度。提高高频和低频部分可以让音乐更明亮有节奏感。

降噪器通过DFT减弱声音中含有噪音的频段，从而实现降噪。我们的手机、电脑的声卡也一般具有降噪器，让你的麦克风输出质量更高。

DFT不仅可以用于音频处理，还可以用于图像处理：给图片加一个水印，直接贴在图片上，很容易被抹掉。一种合理的方法是DFT，在频谱上添加高频水印，再IDFT。

八、代码

#include<cstdio>
#include<cmath>
#include<cctype>

const int MAXN=1e7+10;

inline int Read()
{
    char c=getchar();
    int x=0;
    while(!isdigit(c))
    {
        c=getchar();
    }
    while(isdigit(c))
    {
        x=(x<<3)+(x<<1)+(c^48);
        c=getchar();
    }
    return x;
}
const double Pi=acos(-1.0);
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
}a[MAXN],b[MAXN];

void swap(complex &a,complex &b)
{
    complex c=a;
    a=b;
    b=c;
}

complex operator + (complex a,complex b)
{
    return complex(a.x+b.x , a.y+b.y);
}
complex operator - (complex a,complex b)
{
    return complex(a.x-b.x , a.y-b.y);
}
complex operator * (complex a,complex b)
{
    return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);
}
int N,M;
int l,r[MAXN];
int limit=1;

void FFT(complex *A,int opt)
{
    for(int i=0;i<limit;i++) 
    {
        if(i<r[i])
        {
            swap(A[i],A[r[i]]);
        }
    }
    for(int mid=1;mid<limit;mid<<=1)
    {
        complex Wn( cos(Pi/mid) , opt*sin(Pi/mid) );
        for(int R=mid<<1,j=0;j<limit;j+=R)
        {
            complex w(1,0);
            for(int k=0;k<mid;k++,w=w*Wn) 
            {
                complex x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
        }
    }
}

int main()
{
    int N=Read(),M=Read();
    for(int i=0;i<=N;i++) a[i].x=Read();
    for(int i=0;i<=M;i++) b[i].x=Read();
    while(limit<=N+M) 
    {
        limit<<=1;
        ++l;
    }
    for(int i=0;i<limit;i++)
    {
        r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ;
    }
    FFT(a,1);
    FFT(b,1);
    for(int i=0;i<=limit;i++)
    {
        a[i]=a[i]*b[i];
    }
    FFT(a,-1);
    for(int i=0;i<=N+M;i++)
    {
        printf("%d ",(int)(a[i].x/limit+0.5));
    }
    return 0;
}