快速傅里叶变换FFT

多项式
转换
快速傅里叶变换
- 铺垫
- 定理
- 算法构建
- IFFT
- 递归版FFT&IFFT
迭代版FFT&IFFT
- 蝴蝶效应
Code
后记

多项式

假设有nnn次多项式:
∑i=0naixi\sum_{i=0}^na_ix^ii=0∑naixi
如何表示这个多项式.
首先想到的是用所有系数表示，即:
{an,an−1,an−2,...,a2,a1,a0}\{a_n,a_{n-1},a_{n-2},...,a_2,a_1,a_0\}{an,an−1,an−2,...,a2,a1,a0}
我们将这样的表示方法为系数表示法.
但是，系数表示法在表示多项式相加和相乘很不方便.
譬如又有nnn次表达式:
∑i=0nbixi\sum_{i=0}^nb_ix^ii=0∑nbixi
那么将两式子相乘的时间复杂度是O(n2)O(n^2)O(n2)的.
如果∣a∣|a|∣a∣=10510^5105，这种算法就不是那么优越.
所幸，我们还有另一种表示方法——点值表示法.
设想，如果我们用足够多的值对(x,y) 表示这个函数:
y=∑i=0naixiy=\sum_{i=0}^na_ix^iy=i=0∑naixi
表达式的系数是可以确定的.
更美妙的是，在点值表达式下，多项式加法和乘法是可以在O(n)O(n)O(n)实现的!
比如，有点对
(a1,X1)(a2,X2)...(an,Xn)(a_1,X_1)(a_2,X_2)...(a_n,X_n)(a1,X1)(a2,X2)...(an,Xn)
(a1,Y1)(a2,Y2)...(an,Yn)(a_1,Y_1)(a_2,Y_2)...(a_n,Y_n)(a1,Y1)(a2,Y2)...(an,Yn)
那么相乘得
(a1,X1Y1)(a2,X2Y2)...(an,XnYn)(a_1,X_1Y_1)(a_2,X_2Y_2)...(a_n,X_nY_n)(a1,X1Y1)(a2,X2Y2)...(an,XnYn)
妙啊

转换

但是，我们在实际生活中并不经常使用点值表示法，而是经常使用系数表示法.
于是，便有了两者之间的转换.
系数表示法和点值表示法可以用如下矩阵乘法转换:

[1x1x12...x1n1x2x22...x2n.....1xnxn2...xnn][a0a1...an]=[X1X2...Xn]\left[\begin{array}{cc} 1\ \ x_1\ \ x_1^2\ ...\ x_1^n\\ \\ 1\ \ x_2\ \ x_2^2\ ...\ x_2^n\\ \\ .\ \ \ .\ \ \ .\ \ \ .\ \ \ .\\ \\ 1\ \ x_n\ \ x_n^2\ ...\ x_n^n\\ \end{array}\right] \left[\begin{array}{cc} a_0\\ \\ a_1\\ \\ ...\\ \\ a_n \end{array}\right]= \left[\begin{array}{cc} X_1\\ \\ X_2\\ \\ ...\\ \\ X_n \\ \end{array}\right] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡1 x1 x12 ... x1n1 x2 x22 ... x2n. . . . .1 xn xn2 ... xnn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡a0a1...an⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡X1X2...Xn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
这种变换，叫作离散傅里叶变换(DFT).
反之，由点值表示法退回系数表示法，叫作离散傅里叶逆变换(IDFT).
花了好久点个赞呗
但是发现，这样转换时间复杂度还是O(n2)O(n^2)O(n2)的!
所以，就有了快速傅里叶变换(FFT).

快速傅里叶变换

又称(Fast Fourier Transform)(FFT)，详情见Mr.Du.

快速傅里叶变换是离散傅里叶变换的快速版本.这还用说吗
它可以以优秀的O(nlogn)O(nlogn)O(nlogn)解决这个问题.
举个例子:
我们代入a=1a=1a=1这个特殊值，那么可以得到X=∑i=0naiX=\sum_{i=0}^na_iX=∑i=0nai.
我们代入a=−1a=-1a=−1、a=0a=0a=0也可以得到特殊值.
但是，这样的特殊值总是有限的.
所以，我们将眼光转向别处——复数域.
比如说有一个复数坐标系:

这个圆叫作单位圆没错我们见过，上面的点都可以乘方到1.
如果想要乘方nnn次，那么这些值就是这样分布的:

这是n=8n=8n=8的情况.最好画了不是吗
//下面所有nnn都是222的幂，注意！
我们称这样的数列为{ω1,ω2,...,ωn}\{\omega_1,\omega_2,...,\omega_n\}{ω1,ω2,...,ωn}.(当然是乘方n次后为1的 )
并称这些为 111的nnn次单位根 .

铺垫

由神奇而美妙的欧拉公式:

eix=cos⁡x+isin⁡xe^{ix}=\cos x+i\sin xeix=cosx+isinx

还有复数的极角表达式:

(a,θ1)(b,θ2)=(ab,θ1+θ2)(a, \theta_1)(b,\theta_2)=(ab,\theta_1+\theta_2)(a,θ1)(b,θ2)=(ab,θ1+θ2)

以及复数运算满足交换律和结合律.
我们就可以以这些为铺垫，进行接下来的研究了 ?

定理

下面的定理不再证明，可以结合奆老博客和Mr.Du证明:
//提示:欧拉定理和复平面(上图)是个好东西.

ωn1=cos⁡(2π/n)+isin⁡(2π/n)\omega_n^1=\cos(2\pi/n)+i\sin(2\pi/n)ωn1=cos(2π/n)+isin(2π/n)
对于ωn\omega_nωn序列，ωk=ω1k\omega_k=\omega_1^kωk=ω1k(很实用)
ω2n2k=ωnk\omega_{2n}^{2k}=\omega_n^kω2n2k=ωnk
ωnk+n2=−ωnk\omega_n^{k+\frac{n}{2}}=-\omega_n^kωnk+2n=−ωnk
ωnk+n=ωnk\omega_n^{k+n}=\omega_n^kωnk+n=ωnk

算法构建

经过这么一波操作之后，好像什么都没做…
因为直到现在，我们还是没有在这个数据上构造优化算法.
构造那么神奇的数据，不实现又有何用
要不我们先将一个ωn\omega_nωn代入试试.
比如我们设定
A(n)=∑i=0nanωniA(n)=\sum_{i=0}^na_n\omega_n^iA(n)=i=0∑nanωni
看上去没什么变化…
但是我把它的奇数指数和偶数指数分开，分别写成:
A1(n)=a0+a2ωn+...+anωnn2A_1(n)=a_0+a_2\omega_n+...+a_n\omega_n^\frac{n}{2}A1(n)=a0+a2ωn+...+anωn2n
和
A2(n)=a1ωn+a3ωn2+...+an−1ωnn2A_2(n)=a_1\omega_n+a_3\omega_n^2+...+a_{n-1}\omega_n^\frac{n}{2}A2(n)=a1ωn+a3ωn2+...+an−1ωn2n
//nnn是222的幂.
所以我们震惊的发现，竟然有这样的表达式！
A(x)=A1(x2)+ωnA2(x2)A(x)=A_1(x^2)+\omega_nA_2(x^2)A(x)=A1(x2)+ωnA2(x2)
哇我真棒
但是，这还是没有用…
在仔细研究为什么需要使用复数.
因为，ωnk+n2=−ωnk\omega_n^{k+\frac{n}{2}}=-\omega_n^kωnk+2n=−ωnk!
所以，我们取x=ωnk+n2x=\omega_n^{k+\frac{n}{2}}x=ωnk+2n得:
A(x)=A1(x2)−ωnA2(x2)A(x)=A_1(x^2)-\omega_nA_2(x^2)A(x)=A1(x2)−ωnA2(x2)
整理一下得(k<n2k<\frac{n}{2}k<2n):
A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)A(\omega_n^k)=A_1(\omega_\frac{n}{2}^k)+\omega_n^kA_2(\omega_\frac{n}{2}^k)A(ωnk)=A1(ω2nk)+ωnkA2(ω2nk)
A(ωnk+n2)=A1(ωn2k)−ωnkA2(ωn2k)A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_\frac{n}{2}^k)-\omega_n^kA_2(\omega_\frac{n}{2}^k)A(ωnk+2n)=A1(ω2nk)−ωnkA2(ω2nk)
出现了!FFT!

IFFT

看似好像不那么容易.
但是我们已经会FFT了！
在分支向下的时候，保存每一项的系数,
具体怎么保存呢?
在转换的时候，我们乘的ωn\omega_nωn.
转换回来的时候，就乘上它的共轭复数 What?你不会?
于是，我们就很快活

递归版FFT&IFFT

题目传送门
//为什么叫BF很快揭晓
先随便打个真·暴力上去…

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5;
long long a[N],b[N],c[N];
int n,m;
int main()
{//True BF//O(mn)...scanf("%d%d",&m,&n);for (int i=0;i<=m;i++) scanf("%d",&a[i]);for (int i=0;i<=n;i++) scanf("%d",&b[i]);for (int i=0;i<=m;i++)for (int j=0;j<=n;j++)c[i+j]+=a[i]*b[j];for (int i=0;i<=n+m;i++)printf("%d ",c[i]);return 0;
}

Result:

//哇竟然能有55分！
而上面的诸多公式告诉我们，FFT可以完美解决这好像就是一道FFT模板题…
根据我们的FFT算法，很容易想到递归实现FFT:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);
const int N=5e6+5;
struct complex
{double x,y;complex(double _x=0,double _y=0) {x=_x,y=_y;}
}a[N],b[N];
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);}
void FFT(int n,complex *a,int inv)
{if (n==1) return;int k=n>>1;complex a1[k],a2[k];for (int i=0;i<n;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];//A1和A2FFT(k,a1,inv);//FFT(k,a2,inv);//分治complex Wn=complex(cos(2.*pi/n),inv*sin(2.*pi/n)),w=complex(1,0);//omegaN和单位元for (int i=0;i<k;i++,w=w*Wn){a[i]=a1[i]+w*a2[i];a[i+k]=a1[i]-w*a2[i];//一次推两个}
}
int main()
{int n,m;scanf("%d%d",&n,&m);for (int i=0;i<=n;i++) scanf("%lf",&a[i].x);for (int i=0;i<=m;i++) scanf("%lf",&b[i].x);int p=1;while (p<=n+m) p<<=1;FFT(p,a,1),FFT(p,b,1);//系数转点值for (int i=0;i<=p;i++) a[i]=a[i]*b[i];FFT(p,a,-1);//点值转系数for (int i=0;i<=n+m;i++) printf("%d ",(int)(abs)(a[i].x/p+0.5));return 0;
}

然而我竟然AC了!这是怎么回事?
Result:

这还怎么提出迭代版FFT!别提了

迭代版FFT&IFFT

接下来，我们将它优化一下:
这里是经过FFT分治后的序列

初始序列(id)	0	1	2	3	4	5	6	7
1次操作	0	2	4	6	1	3	5	7
2次操作	0	4	2	6	1	5	3	7
2进制编码	000	100	010	110	001	101	011	111
2进制逆序编码	000	001	010	011	100	101	110	111

我真聪明
可以发现，最后的编码是按逆序升序排列的.
换句话说:每个下标位置处理后的位置是它二进制逆序后的位置!

用这个短小精悍的代码就可以解决一切问题:

//名叫Rader算法
int p=1,L=0;
while (p<=n+m) p<<=1,L++;
for (int i=0;i<p;i++)order[i]=(order[i>>1]>>1)|((i&1)<<(L-1));//(i&1)<<(L-1) 是为了考虑后面的1<<(L-1)个数

真有趣 ?
接下来，我们就更快活开心了 ?

蝴蝶效应

属于小优化的范畴.
主要是因为复数乘法过于缓慢，所以可以预先处理好.
而已qwq.不要想到龙卷风那里去了

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=5e6+5;
const double pi=acos(-1.);
struct complex
{double x,y;complex(double _x=0,double _y=0) {x=_x,y=_y;}
}a[N],b[N];
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+b.x*a.y);}
int order[N];
int p=1,L;
void FFT(complex *a,int inv)
{for (int i=0;i<p;i++) if (i<order[i]) swap(a[i],a[order[i]]);for (int l=1;l<p;l<<=1){complex Wn(cos(pi/l),inv*sin(pi/l));for (int R=l<<1,j=0;j<p;j+=R){complex w(1,0);for (int k=0;k<l;k++,w=w*Wn){complex x=a[j+k],y=w*a[j+l+k];//Ba[j+k]=x+y;a[j+l+k]=x-y;}}}
}
int main()
{int n,m;scanf("%d%d",&n,&m);for (int i=0;i<=n;i++) scanf("%lf",&a[i].x);for (int i=0;i<=m;i++) scanf("%lf",&b[i].x);while (p<=n+m) p<<=1,L++;for (int i=0;i<p;i++)order[i]=(order[i>>1]>>1)|((i&1)<<(L-1));FFT(a,1),FFT(b,1);for (int i=0;i<=p;i++) a[i]=a[i]*b[i];FFT(a,-1);for (int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/p+0.5));return 0;
}

Result:

感谢奆老关注 qwq ?

后记

FFT，又称Fast-Fast-TLE…
这里为什么不直接从<complex><complex><complex>这个库里调用complexcomplexcomplex?
因为…太慢了…

[模板] 快速傅里叶变换(FFT)相关推荐

[模板]快速傅里叶变换(FFT)
Miskcoo大佬的多项式全家桶传送门 rvalue大佬的FFT讲解传送门用途将多项式快速(nlogn)变成点值表达,或将点值表达快速变回系数表达(逆变换),(多数时候)来达到求卷积的目的做法 ...
基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算 ...
基于python的快速傅里叶变换FFT（一）
基于python的快速傅里叶变换FFT(一) FFT可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了.这就是很多信号分析采用FFT变换的原因. ...
MIT 线性代数 Linear Algebra 26:复矩阵，傅里叶矩阵, 快速傅里叶变换 FFT
这一讲我们来讲一下复矩阵.线性代数中,复矩阵是避免不了的话题,因为一个简单实矩阵都有可能有复数特征值. 复矩阵我们着重看一下复矩阵和实矩阵在运算上的区别. 距离首先,一个复数向量的的距离求法发生了 ...
Java中实现快速傅里叶变换FFT
Java中实现快速傅里叶变换FFT 一.概述 1.傅里叶变换(FT) 2.离散傅里叶变换(DFT) 3.快速傅里叶变换(FFT) 1)单位根 2)快速傅里叶变换的思想 3)蝶形图 4)快速傅里叶变换的 ...
OpenCV快速傅里叶变换(FFT)用于图像和视讯流的模糊检测
OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测翻译自[OpenCV Fast Fourier Transform (FFT) for blur detection in images ...
Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析
文章目录 1. 定义 2. 变换和处理 3. 函数 4. 实例演示例1:单频正弦信号(整数周期采样) 例2:单频正弦信号(非整数周期采样) 例3:含有直流分量的单频正弦信号例4:正弦复合信号例5 ...
快速傅里叶变换FFT进行频谱分析（matlab）
快速傅里叶变换FFT进行频谱分析(matlab) 本章摘要:FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了 ...
Java编程实现快速傅里叶变换FFT
快速傅里叶变换的时间复杂度分析 1 快速傅里叶变换FFT 1.1 理论分析 1.1.1 离散傅里叶变换 1.1.2 快速傅里叶变换 1.2 编程实现 1.2.1 算法思想 1.2.2 实验结果 1 快 ...

[模板] 快速傅里叶变换(FFT)