浅谈FFT-快速傅里叶变换

学前须知：

作为一名巨弱的数学竞赛生&高数爱好者，数论知识无疑是我在oi最擅长的领域（没有之一）了。那么我来结合网上的现有资料，以及我的个人见解，书写一篇关于快速傅里叶变换的博客吧。

关于FFT我大约半年前掌握了，现有些许生疏，而且最近学了数学中有关拓扑学的DFT，有了些新的见解，所以写了这篇。

此博客的作用是为了让不会的同学快速入门学习，以及在我本人写的过程中提升自我，无任何商业目的。PS：本人原创意识薄弱可能会从网上找一些现有资料。

此博客重在为初学者提供技巧以及结论的运用，至于理由我会粗略介绍，我认为还是要先知其然，学到手后再考虑证明什么的吧

如果你实在看了还不会的话背板子就行了，不用太过纠结。还有关于关于第二板块的数学知识大家没学过的话跳过就行就是科普一下。

另外此篇是给有数学基础高中生已经自学了高中数学的同学看的。如果对虚数三角函数等知识不熟悉可以自行学习，此论文中不再赘述了。

1.FFT的作用以及功能

相信大家都曾听闻过FFT，它是一种可以优化高精度乘法的高逼格算法。

首先我们要知道FFT是由DFT以更高效，快速计算的方式得到的。在介绍FFT之前我会再讲一下关于DTF（离散傅里叶变换）的数学意义，这是其他的博客都不介绍的内容。为什么介绍呢，这是因为可以让你了解数学中离散傅里叶变换是如何把信号从时间域转到频率域的，正如我们的现在的计算机的DFT也是用的这种转换。即把系数转成点值（过会会介绍具体）。如果数学没有到大学数学本科or研究生水平的这段可以跳过了，感兴趣的可以了解下。

2.数学上的DFT（可以忽略）

关于DFT：

如果你看到了这里，首先我默认你会狄拉克函数的基础内容

1.狄拉克函数

此函数数学表达式为：

其时域波是周期为Ts的单位脉冲序列，也称之为理想抽样函数。

由于满足了周期性，我们可以用傅里叶级数求其频谱，傅里叶级数系数为：

其傅里叶展开为：

可见其频率分布值在的离散频率点，频率分量幅值均为1/Ts。

引进冲击函数的概念后，一般的周期函数也可求傅里叶变换，如下（cn为周期函数的傅里叶级数）

用狄拉克梳状函数表示离散信号和离散频谱，将给信号分析和处理带来极大的方便。

2.离散时间傅里叶（DTFT）

既然你都看到这里了，我又默认你会傅里叶变换（调和分析范畴，非战斗人员请撤离）

既然我们会了处理连续系统的，那离散系统怎么处理呢？

进行采样，频率为fs，冲击采样序列为：

取样后：

连续信号的傅里叶变换定义：

采样后：

同理于：

由于狄拉克函数筛选性质可知：

连续信号的傅里叶变换如下：

3.离散傅里叶变换（DFT）

我们容易知道，在当前的世界，信号的频谱是不连续的，所以在DTFT中我们的频率也需要离散化处理，所需要的就是DFT了。即具有周期特性离散信号的傅里叶级数（就是将无限长的离散限号进行截短至N个采样点，然后将这个N个采样点进行周期延拓 ）

DFT：

且

IDFT（逆变换）：

且

这样DFT就介绍完了，至于证明还是算了（毕竟这里不是重点）。其实FFT本质也是DFT，只不过是利用DFT内部蕴藏的规律的一种简化算法罢了。

3.关于多项式

首先想想我们要解决的问题是什么----两个多项式想乘

容易知道：这两个多项式我们可以有不同的表示方法

1.系数表示法

十分好理解的最朴素的表达方法

如果我们用这种表示法暴力算的话，复杂度肯定是O（n^2）的，显然TLE。

2.点值表示法

我们不妨设两多项式分别为f(x）与g(x)

具体表示比如，我选一个x0分别代入f（x）与g（x），那么f（x0）的值也就是原多项式a1的第0项的值，g（x0）的值也就是原多项式a2的第零项的值。如果不理解的话我列一下图表吧。

f[x]={(x0,f[x0]),(x1,f[x1]).....(xn,f(xn)]

g[x]={x0,g[x0]),(x1,g[x1])......(xn,g[xn])

也就是说，当我们去一个x时，经过两种不同的函数关系，可以得到不同的y值。那么我们只需要y值对于相乘就好了。复杂度O（n）

4.关于oi中的DFT

对于一般的多项式来说呢，我们当然可以随便去n个x值代入计算。

也就是暴力计算出每个x的值，然后再代入y中。显然这样转换复杂度是o（n^2）

怎么办呢？

考虑DFT了。

我们看到我写的第二部分的数学中的DFT可以发现，其本质原理就是就是将无限长的离散限号进行截短至N个采样点，然后将这个N个采样点进行周期延拓。如果不知道的话没关系，其实就是去一些“神奇”的x代入。

他规定了点值表示中的n个x为n个模长为一的复数，如图：

然后要把这个单位圆N等分，如图当N=8时

这些是当n=8时n等分后要取的点。记第k个点代表的复数值是，可以知道。

那么此时称为n次单位根。而且每个单位根都可以解出来：

那么这些Wn0，Wn1.......Wn n-1就是我们要代入的x0，x1 .......xn-1。（这里的W就是欧米伽蛤）

5.DFT优化后的FFT-快速傅里叶变换

你会发现其实你这样暴力计算这些单位根后复杂度还是O（n^2）.....

通过单位根的性质，考虑分治优化。

考虑多项式A（x）=a0+a1x+a2x^2+a3x^3+.....+an-1^(n-1)

我们可以根据下标的奇数性分成两组，注意是“下标”的。

可以发现两边非常相似。

令A0（x）=a0 x^0+a2x+a4x^w+......+a(n-2)x^(n/2-1)

令A1（x）=a1 x^0+a3x+a5x^2+.....+a(n-1)x^(n/2-1)

显然可以得到A（x）=A0（x^2）+x*A1（x^2）

蝴蝶变换：当我们关于A1（x）与A2（x）在的点值求出后，可以O（n）求出A（x）在的点值了。

我们在处理过程中不断递归分治，知道n=1时return。可以知道复杂度是o（nlog2n）的。

于是我在网上找了一个朴素的FFT板子

#include<complex>
#define cp complex<double>

void fft(cp *a,int n,int inv)
{
    if (n==1)return;
    int mid=n/2;
    static cp b[MAXN];
    fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1];
    fo(i,0,n-1)a[i]=b[i];
    fft(a,mid,inv),fft(a+mid,mid,inv);
    fo(i,0,mid-1)
    {
        cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));
        b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid];
    }
    fo(i,0,n-1)a[i]=b[i];
}

嗯，TLE了。又双叒叕优化。

6.迭代版的FFT

可以发现，我们在进行FFT时会把各个系数分别放在两侧。一个系数的原来位置后最终位置的规律：

考虑二进制。

容易发现每个位置分治后其实就是二进制翻转过来的。所以我们把每个数放在最后的位置，然后不断向上还原同时求出点值就可以了。

复杂度根号O（nlog2n）

FFT板子：

void fft(cp *a,int n,int inv)
{
    int bit=0;
    while ((1<<bit)<n)bit++;
    fo(i,0,n-1)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
    }
    for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
    {
        cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根，pi的系数2已经约掉了
        for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度，i是合并到了哪一位
        {
            cp omega(1,0);
            for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分，得到右半部分的答案
            {
                cp x=a[i+j],y=omega*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;//蝴蝶变换
            }
        }
    }
}