FFT(快速傅里叶变换)——从入门到放弃

2018-03-29 14:46:43


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的定义

$i=\sqrt{-1}$,这就是$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$$

  1. $\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;
}