十分简明易懂的FFT（快速傅里叶变换）

最新推荐文章于 2024-09-06 17:21:28 发布

Thomas会写字

最新推荐文章于 2024-09-06 17:21:28 发布

阅读量3.3w

点赞数 88

分类专栏：数学计算文章标签： c# 算法

原文链接：https://blog.csdn.net/enjoy_pascal/article/details/81478582/

版权

数学计算专栏收录该内容

16 篇文章 1 订阅

订阅专栏

FFT是一种用于计算离散傅里叶变换的高效算法，由库利和图基在1965年提出。文章介绍了FFT的基本概念，包括系数表示法、点值表示法，以及复数和单位根在FFT中的作用。通过DFT和IFFT的概念，解释了如何使用FFT进行快速多项式乘法，并提供了FFT的分治实现和优化。

摘要由CSDN通过智能技术生成

FFT前言

快速傅里叶变换 (fast Fourier transform)，即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。

FFT（Fast Fourier Transformation） 是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

——百度百科

FFT（Fast Fourier Transformation），中文名快速傅里叶变换，用来加速多项式乘法朴素高精度乘法时间 $O(n^2)$ )，但FFT能 $O(n log_2 n)$ 的时间解决。
FFT名字逼格高，也难懂，其他教程写得让人看不太懂，于是自己随便写一下

建议对复数、三角函数相关知识有所耳闻（不会也无所谓）

下面难懂的点我会从网上盗

多项式的系数表示法和点值表示法FFT其实是一个用 $O(n log_2 n)$ 的时间将一个用系数表示的多项式转换成它的点值表示的算法
多项式的系数表示和点值表示可以互相转换

系数表示法

一个n-1次n项多项式 $f(x)$ 可以表示为 $f(x)=\sum^{n-1}_{i=0}a_ix^i$

也可以用每一项的系数来表示 $f(x)$ ，即 $f(x)=\{a_0,a_1,a_2,...,a_{n-1} \}$

这就是系数表示法，也就是平时~~数学课上~~用的方法

点值表示法

把多项式放到平面直角坐标系里面，看成一个函数
把 $n$ 个不同的 $x$ 代入，会得出 $n$ 个不同的 $y$ ，在坐标系内就是 $n$ 个不同的点
那么这 $n$ 个点唯一确定该多项式，也就是有且仅有一个多项式满足 $\forall k,f(x_k)=y_k$
理由很简单，把n nn条式子联立起来成为一个有n条方程的n元方程组，每一项的系数都可以解出来

那么 $f(x)$ 还可以用 $f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}$ 来表示，这就是点值表示法。

高精度乘法下两种多项式表示法的区别

对于两个用系数表示的多项式，我们把它们相乘，设两个多项式分别为 $A(x),B(x)$ ，我们要枚 $A$ 每一位的系数与 $B$ 每一位的系数相乘，那么系数表示法做多项式乘法时间复杂度 $O(n^2)$ 。

但两个用点值表示的多项式相乘，只需要 $O(n)$ 的时间。

什么意思呢？

设两个点值多项式分别为：
$f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}$

$g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),...,(x_{n-1},g(x_{n-1}))\}$

设它们的乘积是 $h(x)$ ，那么
$h(x)=\{(x_0,f(x_0)·g(x_0)),(x_1,f(x_1)·g(x_1)),...,(x_{n-1},f(x_{n-1})·g(x_{n-1}))\}$

所以这里的时间复杂度只有一个枚举的 $O(n)$

突然感觉高精度乘法能 $O(n)$ 解决一堆题？
但是朴素的系数表示法转点值表示法的算法还是 $O(n^2)$ 的复杂操作
朴素系数转点值的算法叫DFT（离散傅里叶变换），点值转系数叫IDFT（离散傅里叶逆变换）
难道高精度乘法只能 $O(n^2)$ 了吗？

DFT前置知识&技能

复数

毕竟高中有所以不多说

我们把形如 $z=a+bi$ （ $a,b$ 均为实数）的数称为复数，其中 $a$ 称为实部， $b$ 称为虚部， $i$ 称为虚数单位。当虚部等于零时，这个复数可以视为实数；当 $z$ 的虚部不等于零时，实部等于零时，常称 $z$ 为纯虚数。复数域是实数域的代数闭包，也即任何复系数多项式在复数域中总有根。复数是由意大利米兰学者卡当在十六世纪首次引入，经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作，此概念逐渐为数学家所接受。
——百度百科

初中数学老师会告诉你没有 $\sqrt{-1}$ ，但仅限 $R$

扩展至复数集 $C$ ，定义 $i^2=-1$ ，一个复数 $z$ 可以表示为 $z=a+bi(a,b\in R)$

其中 $a$ 称为实部， $b$ 称为虚部， $i$ 称为虚数单位

在复数集中就可以用i ii表示负数的平方根，如 $\sqrt{-7}=\sqrt{7}i$

还可以把复数看成复平面直角坐标系上的一个点，比如下面

$x$ 轴就是实数集中的坐标轴， $y$ 轴就是虚数单位 $i$ 轴

这个点 $(2,3)$ 表示的复数就是 $2+3i$ ，或者想象它代表的向量为 $(2,3)$

其实我们还可以自己想象 ~~（其实没有这种表达方式）~~ 它可以表示为 $(\sqrt{13},\theta)$

一个复数 $z$ 的模定义为它到原点的距离，记为 $|z|=\sqrt{a^2+b^2}$

一个复数 $z=a+bi$ 的共轭复数为 $a-bi$ （虚部取反），记为 $\overline{z}=a-bi$

复数的运算

复数不像点或向量，它和实数一样可以进行四则运算

设两个复数分别为 $z_1=a+bi,z_2=c+di$ ，那么

$z_1+z_2=(a+c)+(b+d)i$

$z_1z_2=(ac-bd)+(ad+bc)i$

复数相加也满足平行四边形法则

这张是从网上盗的

即 $AB+AD=AC$

复数相乘还有一个值得注意的小性质

$(a_1,\theta_1)*(a_2,\theta_2)=(a_1a_2,\theta_1+\theta_2)$

即模长相乘，极角相加

DFT（离散傅里叶变换）

一定注意从这里开始所有的 $n$ 都默认为 2 的整数次幂

对于任意系数多项式转点值，当然可以随便取任意 $n$ 个 $x$ 值代入计算，但是暴力计算 $x_k^0,x_k^1,...,x_k^{n-1}(k\in[0,n))$ 当然是 $O(n^2)$ 的时间。其实可以代入一组神奇的 $x$ ，代入以后不用做那么多的次方运算，这些 $x$ 当然不是乱取的，而且取这些 $x$ 值应该就是傅里叶的主意了。

考虑一下，如果我们代入一些 $x$ ，使每个 $x$ 的若干次方等于 1，我们就不用做全部的次方运算了 $\pm 1$ 是可以的，考虑虚数的话 $\pm i$ 也可以，但只有这四个数远远不够

傅里叶说：这个圆圈上面的点都可以做到

以原点为圆心，画一个半径为 1 的单位圆，那么单位圆上所有的点都可以经过若干次次方得到 1

傅里叶说还要把它给 $n$ 等分了，比如此时 $n=8$

橙色点即为 $n=8$ 时要取的点，从 $(1,0)$ 点开始，逆时针从 0 号开始标号，标到 7 号，记编号为 $k$ 的点代表的复数值为 $\omega_n^k$ ，那么由模长相乘，极角相加可知 $(\omega_n^1)^k=\omega_n^k$ ，其中 $\omega_n^1$ 称为 $n$ 次单位根，而且每一个 $\omega$ 都可以求出 ~~（我三角函数不好）~~

$\omega_n^k=\cos{k\over n}2\pi +i\sin{k\over n} 2\pi$

或者说也可以这样解释这条式子

注意 $sin^2\theta+cos^2\theta=1$ 什么的，就容易理解了

那么 $\omega^0_n,\omega^1_n,...,\omega^{n-1}_n$ 即为我们要代入的 $x_0,x_1,...,x_{n-1}$

单位根的一些性质

推FFT的过程中需要用到 $\omega$ 的一些性质

$\omega^k_n=\omega^{2k}_{2n}$

它们表示的点（或向量）表示的复数是相同的
证明
$\omega^k_n=\cos{k\over n}2\pi+i\sin{k\over n} 2\pi=\cos{2k\over 2n}2\pi+i\sin{2k\over 2n} 2\pi=\omega^{2k}_{2n}$

$\omega^{k+{n \over 2}}_n=-\omega_n^k$

它们表示的点关于原点对称，所表示的复数实部相反，所表示的向量等大反向
证明
$\omega^{n\over 2}_n=\cos{{n\over 2}\over n}2\pi+i\sin{{n\over 2}\over n}2\pi=\cos\pi+i\sin\pi=-1$
（这个东西和 $e^{ix}=\cos x+i\sin x$ 与 $e^{i\pi}+1=0$ 有点关系，~~我不会就不讲了~~）

$\omega^0_n=\omega^n_n$

都等于 $1$ ，或 $1+0i$

FFT（快速傅里叶变换）

虽然 DFT 搞出来一堆很牛的 $\omega$ 作为代入多项式的 $x$ 值

但只是代入这类特殊 $x$ 值法的变换叫做 DFT 而已，还是要代入单位根暴力计算

DFT 还是暴力 $O(n^2)$ …

但 DFT 可以分治来做，于是 FFT（快速傅里叶变换） 就出来了

首先设一个多项式 $A(x)$

$A(x)=\sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$

按 $A(x)$ 下标的奇偶性把 $A(x)$ 分成两半，右边再提一个 $x$

$A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})$

$=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})$

两边好像非常相似，于是再设两个多项式 $A_1(x),A_2(x)$ ，令

$A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{n\over 2}-1}$

$A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{n \over 2}-1}$

比较明显得出

$A(x)=A_1(x^2)+xA_2(x^2)$

再设 $k< {n\over 2}$ ，把 $\omega^k_n$ 作为 $x$ 代入 $A(x)$ ~~（接下来几步变换要多想想单位根的性质）~~

$A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2)$

$=A_1(\omega^{2k}_n)+\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})+\omega^k_nA_2(\omega^k_{n\over 2})$

那么对于 $A(\omega^{k+{n\over2}}_n)$ 的话，代入 $\omega^{k+{n \over 2}}_n$

$A(\omega^{k+{n\over 2}}_n)=A_1(\omega^{2k+n}_n)+\omega^{k+{n\over 2}}_nA_2(\omega^{2k+n}_n)$

$=A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n)$

$=A_1(\omega^{2k}_n)-\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})-\omega^k_nA_2(\omega^k_{n\over2})$

发现了什么？

$A(\omega^k_n)$ 和 $A(\omega^{k+{n\over2}}_n)$ 两个多项式后面一坨东西只有符号不同
就是说，如果已知 $A_1(\omega^k_{n\over 2})$ 和 $A_2(\omega^k_{n\over 2})$ 的值，我们就可以同时知道 $A(\omega^k_n)$ 和 $A(\omega^{k+{n\over2}}_n)$ 的值

现在我们就可以递归分治来搞 $FFT$ 了

每一次回溯时只扫当前前面一半的序列，即可得出后面一半序列的答案

$n==1$ 时只有一个常数项，直接 $return$

时间复杂度 $O(n\log_2n)$

IFFT（快速傅里叶逆变换）

想一下，我们不仅要会 $FFT$ ，还要会IFFT（快速傅里叶逆变换）

我们把两个多项式相乘 （也叫求卷积），做完两遍 $FFT$ 也知道了积的多项式的点值表示

可我们平时用系数表示的多项式，点值表示没有意义

怎么把点值表示的多项式快速转回系数表示法？
$IDFT$ 暴力 $O(n^2)$ 做？其实也可以用 $FFT$ 用 $O(n\log_2n)$ 的时间搞

你有没有想过为什么傅里叶是把 $\omega^k_n$ 作为 $x$ 代入而不是别的什么数？

原因是有的~~但是有我也看不懂~~

所以只用记住一个结论

一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以 $n$ 即为原多项式的每一项系数

意思就是说 $FFT$ 和 $IFFT$ 可以一起搞

朴素版FFT板子

$c++$ 有自带的复数模板 $complex$ 库
$a.real()$ 即表示复数 $a$ 的实部

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

void fft(cp *a,int n,int inv)//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));//inv取决是否取共轭复数
        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];
}

两个多项式 $a,b$ 相乘再转系数多项式 $c$ ，通常只用打这么一小段

	cp a[MAXN],b[MAXN];
	int c[MAXN];
	fft(a,n,1),fft(b,n,1);//1系数转点值
	fo(i,0,n-1)a[i]*=b[i];
	fft(a,n,-1);//-1点值转系数
	fo(i,0,n-1)c[i]=(int)(a[i].real()/n+0.5);//注意精度

很明显，FFT只能处理n nn为2 22的整数次幂的多项式

所以在 $FFT$ 前一定要把 $n$ 调到 2 的次幂

这个板子看着好像很优美，但是

递归常数太大，要考虑优化…

FFTの优化——迭代版FFT

~~这个图也是盗的~~

这个很容易发现点什么吧？

每个位置分治后的最终位置为其二进制翻转后得到的位置

这样的话我们可以先把原序列变换好，把每个数放在最终的位置上，再一步一步向上合并

一句话就可以 $O(n)$ 预处理第i ii位最终的位置 $rev[i]$

fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));

~~至于蝴蝶变换它死了其实是我不会~~

真·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;//这个就是蝴蝶变换什么的
            }
        }
    }
}

~~这个板子好像不是那么好背~~

至少这个板子已经很优美了

FFT后记

~~本人版权意识薄弱……~~

本博客部分知识学习于

https://www.cnblogs.com/RabbitHu/p/FFT.html
https://www.cnblogs.com/zwfymqz/p/8244902.html?mType=Group#_label3
https://blog.csdn.net/ggn_2015/article/details/68922404