FFT算法讲解——麻麻我终于会FFT了!
FFT——快速傅里叶变换
这块不写东西空荡荡的,我决定还是把FFT的定义给贴上吧
FFT(Fast Fourier Transformation)是离散傅氏变换(DFT)的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。
这三段话其实一点用也没有
FFT是干什么的
FFT在算法竞赛中就有一个用途:加速多项式乘法(暴言)
简单来说,形如 a0X0+a1X1+a2X2+⋯+anXna_0X^0+a_1X^1+a_2X^2+⋯+a_nX^na0X0+a1X1+a2X2+⋯+anXn 的代数表达式叫做多项式,可以记作f(X)=a0X0+a1X1+a2X2+⋯+anXnf(X)=a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^nf(X)=a0X0+a1X1+a2X2+⋯+anXn,其中a0,a1,⋯,ana_0,a_1,⋯,a_na0,a1,⋯,an叫做多项式的系数,XXX是一个不定元(就是不可以合并),不表示任何值,不定元在多项式中最大项的次数称作多项式的次数
如果我们当前有两个多项式f(X),g(X)f(X),g(X)f(X),g(X),现在要把他们乘起来(求卷积),最朴素的做法就是
∑i=02n−1(∑j=0iaj∗bi−j)∗xi\sum \limits _{i=0}^{2n-1} (\sum \limits _{j=0} ^{i} a_j*b_{i-j})*x^ii=0∑2n−1(j=0∑iaj∗bi−j)∗xi
这样的复杂度是Θ(n2)\varTheta(n^2)Θ(n2)的,十分不美观,FFT就是要将这个过程优化为Θ(nlogn)\varTheta(n \log n)Θ(nlogn)
前置技能
多项式
见上文
复数
复数形如a+bia+bia+bi,其中i=−1i=\sqrt {-1}i=−1
aaa叫作复数的实部,bibibi叫做复数的虚部
复数(a1+b1i)∗(a2+b2i)(a_1+b_1i)*(a_2+b_2i)(a1+b1i)∗(a2+b2i)相乘的值,即a1a2−b1b2+(a1b2+a2b1)ia_1a_2-b_1b_2+(a_1b_2+a_2b_1)ia1a2−b1b2+(a1b2+a2b1)i也是一个复数,同时我们也得到了复数的乘法法则
复数c+dic+dic+di可以用这种方式表示出来
复数乘法的在复平面中表现为辐角相加,模长相乘
单位根
复数www满足wn=1w^n=1wn=1称作www是nnn次单位根,下图包含了所有的888次单位根(图中圆的半径是1)
同样的,下图是所有的4次单位根
聪明的你也许已经发现了单位根的些许性质,即
w2n2m=wnmw_{2n}^{2m}=w_{n}^{m}w2n2m=wnm
wnm=−wnm+n2w_n^m=-w_n^{m+\frac{n}{2}}wnm=−wnm+2n
这两个要记住,一会很有用
多项式的系数表达法
我们有多项式f(X)=a0X0+a1X1+a2X2+⋯+anXnf(X)=a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^nf(X)=a0X0+a1X1+a2X2+⋯+anXn,令a⃗=(a0,a1,a2,...,an)\vec{a}=(a_0,a_1,a_2,...,a_n)a=(a0,a1,a2,...,an),则称A(X)A(X)A(X)为多项式f(X)f(X)f(X)的系数表示法
在系数表示法下,计算多项式乘法是Θ(n2)\varTheta (n^2)Θ(n2)的
多项式的点值表达法
任取n+1n+1n+1个互不相同的S={p1,p2,...,pn+1}S=\{p_1,p_2,...,p_{n+1}\}S={p1,p2,...,pn+1},对f(X)f(X)f(X)分别求值得到f(p1),f(p2),...,f(pn+1)f(p_1),f(p_2),...,f(p_{n+1})f(p1),f(p2),...,f(pn+1),那么称A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}A(X)=\{(p_1,f(p_1)),(p_2,f(p_2)),...,(p_{n+1},f(p_{n+1}))\}A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}为多项式f(X)f(X)f(X)在SSS下的点值表示法
可以把多项式想象成一个nnn次函数,点值表示法就是取SSS下每一个横坐标时对应的点,因为nnn次函数可以由n+1n+1n+1个点确定下来(可以将每一个点列一个nnn次方程),所以nnn维点值与nnn维系数一一对应
更重要的一点,点值表示法下的乘法运算获得了简化
两个多项式PPP,QQQ分别取点(x,y1)(x,y_1)(x,y1)和(x,y2)(x,y_2)(x,y2),P∗QP*QP∗Q就会取到点(x,y1∗y2)(x,y_1*y_2)(x,y1∗y2);
令C=P∗QC=P*QC=P∗Q,因为C(X)=P(X)∗Q(X)C(X)=P(X)*Q(X)C(X)=P(X)∗Q(X),所以C(x)=P(x)∗Q(x)C(x)=P(x)*Q(x)C(x)=P(x)∗Q(x),即C(x)=y1∗y2C(x)=y_1*y_2C(x)=y1∗y2
所以在点值表示法下,计算多项式乘法是Θ(n)\varTheta(n)Θ(n)的
FFT的具体过程
FFT就是将系数表示法转化成点值表示法相乘,再由点值表示法转化为系数表示法的过程,第一个过程叫做求值(DFT),第二个过程叫做插值(IDFT)
求值
还记得我们之前提到的单位根吗?再回顾一下:
w2n2m=wnmw_{2n}^{2m}=w_{n}^{m}w2n2m=wnm
wnm=−wnm+n2w_n^m=-w_n^{m+\frac{n}{2}}wnm=−wnm+2n
想要求出一个多项式的点值表示法,需要选出n+1n+1n+1个数分别带入到多项式里面,带入一个数的复杂度是Θ(n)\varTheta(n)Θ(n)的,那么总复杂度就是Θ(n2)\varTheta(n^2)Θ(n2)的,因为单位根有上面两个优美的性质,所以我们尝试可以取nnn次单位根组成SSS,看看能不能加速我们的运算
设A0(X)A_0(X)A0(X)为A(X)A(X)A(X)偶次项的和,设A1(X)A_1(X)A1(X)为A(X)A(X)A(X)奇次项的和,即
A0(X)=a0x0+a2x1+...+an−1xn/2A_0(X)=a_0x^0+a_2x^1+...+a_{n-1}x^{n/2}A0(X)=a0x0+a2x1+...+an−1xn/2
A1(X)=a1x0+a3x1+...+an−2xn/2A_1(X)=a_1x^0+a_3x^1+...+a_{n-2}x^{n/2}A1(X)=a1x0+a3x1+...+an−2xn/2
因为A(wnm)=a0wn0+a1wnm+a2wn2m+a3wn3m+...+an−1wn(n−1)∗m+anwnnmA(w_n^m)=a_0w_n^0+a_1w_n^m+a_2w_n^{2m}+a_3w_n^{3m}+...+a_{n-1}w_n^{(n-1)*m}+a_nw_n^{nm}A(wnm)=a0wn0+a1wnm+a2wn2m+a3wn3m+...+an−1wn(n−1)∗m+anwnnm
所以有
也就是说,只要有了A0(X)A_0(X)A0(X)和A1(X)A_1(X)A1(X)的点值表示,就能在Θ(n)\varTheta(n)Θ(n)时间算出A(X)A(X)A(X)的点值表示,对于当前层确定的位置iii,就可以用下一层的两个值更新当前的值,我们称这个操作为“蝴蝶变换”
因为这个过程一定要求每层都可以分成两大小相等的部分,所以多项式最高次项一定是2p2^p2p(p∈Np\in Np∈N)次方,如果不是的话,直接在最高次项补零就可以啦
于是我们有了递归的写法:
void FFT(Complex* a,int len){if(len==1) return;Complex* a0=new Complex[len/2];Complex* a1=new Complex[len/2];for(int i=0;i<len;i+=2){a0[i/2]=a[i];a1[i/2]=a[i+1];}FFT(a0,len/2);FFT(a1,len/2);Complex wn(cos(2*Pi/len),sin(2*Pi/len));Complex w(1,0);for(int i=0;i<(len/2);i++){a[i]=a0[i]+w*a1[i];a[i+len/2]=a0[i]-w*a1[i];w=w*wn;}return;
}
但递归版的FFT常数巨大,实现起来比较复杂,于是又有了迭代的写法
重新考虑下递归FFT的过程,在第iii次求解中,我们将所有元素二进制iii位为000的放在了左面,iii位为111的放在了右面,事实上,每个元素最终到的是他二进制颠倒过来的位置
再拿一张别人的图
迭代写法:
inline void DFT(Complex a[]){for(int i=0;i<len;i++)///pos[i]代表反转后的位置if(i<pos[i])swap(a[i],a[pos[i]]);for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){///len代表多项式最高次项///第一层i是枚举合并到了哪一层。Complex wm(cos(2.0*pi/i),sin(2.0*pi/i));for(int j=0;j<len;j+=i){///第二层j是枚举合并区间。Complex w(1,0);for(int k=j;k<j+mid;k++,w=w*wm){///第三层k是枚举区间内的下标。Complex l=a[k],r=w*a[k+mid];a[k]=l+r;a[k+mid]=l-r;}}}return;
}
插值
有人证出来插值只要将所有wnmw_n^mwnm换成wnm+n2w_n^{m+\frac{n}{2}}wnm+2n,也就是所有的虚部取相反数,再将最终结果除以lenlenlen,就是插值的过程了
究竟为什么?我觉得人生一定要有点遗憾可以参考这里,我就不多说了
支持插值的迭代写法:
const double DFT=2.0,IDFT=-2.0;///进行求值,第二个参数传DFT,插值传IDFT
inline void FFT(Complex a[],double mode){for(int i=0;i<len;i++)if(i<pos[i])swap(a[i],a[pos[i]]);for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){Complex wm(cos(2.0*pi/i),sin(mode*pi/i));for(int j=0;j<len;j+=i){Complex w(1,0);for(int k=j;k<j+mid;k++,w=w*wm){Complex l=a[k],r=w*a[k+mid];a[k]=l+r;a[k+mid]=l-r;}}}if(mode==IDFT)for(int i=0;i<len;i++)a[i].x/=len;return;
}
现在,你已经会写FFT了!
生成函数
小AAA有aia_iai个价值为AiA_iAi的物品,小BBB有bib_ibi个价值为AiA_iAi的物品,求用两个组成价值为cic_ici的方案数
生成函数可以解决上面的这个问题,构造两个多项式,第iii项表示价值为iii的物品有多少个,对两个人分别构造,乘在一起的多项式就代表所有的方案数了
NTT——快速数论变换
当FFT需要取模时,复数的存在就特别尴尬了,有没有什么东西可以代替单位根呢?
原根性质:假设ggg是素数ppp的一个原根,则g1,g2,...,g(p−1)g^1,g^2,...,g^{(p-1)}g1,g2,...,g(p−1)在模ppp意义下两两不同,且结果恰好为111~p−1p-1p−1
wnn≡g(p−1)≡1(modp)w_n^n≡g^{(p-1)}≡1(mod\ p)wnn≡g(p−1)≡1(mod p)(中间那个是费马小定理)
所以可以把g(p−1)/ng^{(p-1)/n}g(p−1)/n看成wn1w_n^1wn1的等价。
但这种质数必须是NTT质数(费马质数),即(p−1)(p-1)(p−1)有超过序列长度的2的正整数幂因子的质数
具体证明还是看这里吧
参考资料:
Pick’s blog
Fenghr的博客
Miskcoo’s Space
FFT算法讲解——麻麻我终于会FFT了!相关推荐
- FFT算法再学以及终于理解
前言 人生如逆旅,我亦是行人. 一.FFT FFT(Fast Fourier Transformation),中文名快速傅里叶变换,用来 加速多项式乘法 ,就是用来降低算法的时间复杂度的,将时间复杂度 ...
- 数字信号处理(DTFT与DFT、DFS的详细讲解以及FFT算法)
DTFT与DFT.DFS的详细讲解以及FFT算法 DTFT与DFT.DFS的区别在哪里呢? 离散傅里叶级数DFS 离散傅里叶变换DFT 有限长序列的线性卷积和循环卷积 利用DFT做连续信号的频谱分析 ...
- 【总结】FFT算法在信息竞赛中的应用
FFT算法本身就是一种优化,优化(类似)卷积运算的时间复杂度 (卷积:∑i,jai∗bj−i\sum_{i,j}a_i*b_{j-i}∑i,jai∗bj−i). FFT的本质,其实是利用复数的一 ...
- FFT算法8点12位硬件实现 (verilog)
FFT算法8点12位硬件实现 (verilog) 1 一.功能描述: 1 二.设计结构: 2 三.设计模块介绍 3 1.蝶形运算(第一级) 3 2.矢量角度旋转(W) 4 3.CORDIC 结果处理 ...
- FFT算法的完整DSP实现
傅里叶变换或者FFT的理论参考: [1] http://www.dspguide.com/ch12/2.htm The Scientist and Engineer's Guide to Digita ...
- 大整数乘法---FFT算法
//迭代FFT的乘法方法 // POJ 1405 Heritage /** * input data mode: * the number array 1,2,3,4 use base U = ...
- matlab fft函数说明_【V2.0更新】基于FFT算法的MTALAB傅里叶级数3D可视化
最近这份代码受到很多朋友关注,在此一并感谢!由于之前写1.0版本代码时的时间有限,当时的知识储备也不甚完善,所以只是做出了基本功能. 近期对这份代码进行了更新,优化了代码结构,并加入了FFT等算法让其 ...
- JavaScript实现快速傅立叶变换FFT算法(附完整源码)
JavaScript实现快速傅立叶变换FFT算法(附完整源码) radianToDegree.js完整源代码 ComplexNumber.js完整源代码 bitLength.js完整源代码 fastF ...
- c代码实现 ifft运算_fft算法c语言_matlab fft算法_ifft c语言
FFT快速算法C程序_工学_高等教育_教育专区.电子信息工程综合课程设计报告书 DSP 课程设计 报告 题学 目: 院: FFT 快速算法 C 程序 计算机与信息工程学院 09 ... fft算法代码 ...
- 模意义下的FFT算法
//写在前面 单就FFT算法来说的话,下面只给出个人认为比较重要的推导,详细的介绍可参考 FFT算法学习笔记 令v[n]是长度为2N的实序列,V[k]表示该实序列的2N点DFT.定义两个长度为N的实序 ...
最新文章
- 90后教授:回国是用数学计算出的“最优解”
- ABAP入门程序,你会了嘛?
- qq分享提示设备未授权_QQ帐号已经可以注销了,过去几天,第一批尝试的人已经放弃了!...
- 从源码深处体验Spring核心技术--面试中IOC那些鲜为人知的细节
- 德 梅齐里亚克的砝码问题matlab,德梅齐里亚克砝码问题之解
- leetcode129. 求根到叶子节点数字之和(dfs)
- 番外篇:根据学习程度划分程序员的级别
- 【CUDA学习】计时方法
- 23.6. Functions
- Python的文件读取操作
- 前端获取后台布尔类型_教育平台项目前端:视频讲解
- win10彻底关闭自动更新
- 阿里云 EMAS Serverless 重磅发布
- html移动端注册流程,登录和注册移动端.html
- mysql java驱动 ibm_Java 通过JDBC连接Mysql数据库
- PHP 环境搭建(win7+php5.6+apache或nginx)
- 手机裂脑纪:中国式审美还有救吗?
- 信息学竞赛中可能用到的初等数论
- js一维数组,api,二维数组
- dell服务器启动顺序如何设置_Dell PowerEdge 服务器启动指南
热门文章
- 用matlab进行函数插值的几种方法
- macOS录制系统声音及麦克风的三种方法
- PHP字符串解析函数
- linux打开dwg格式文件怎么打开软件,DWG 文件扩展名: 它是什么以及如何打开它?...
- 软件显示服务器端没有启动,打开软件显示无连接服务器,双击服务器显示如图...
- ios12越狱自签需要联网_从越狱的iOS切换到Android? 这是你需要知道的
- 新手经常忽略的嵌入式基础知识点,你都掌握了吗?
- 远程小组软件开发过程(3):人
- 和差测角天线方向图仿真matlab 含代码
- 1050ti显卡安装cuda