FFT算法讲解——麻麻我终于会FFT了!

article/2025/10/19 1:04:25

FFT——快速傅里叶变换

这块不写东西空荡荡的,我决定还是把FFT的定义给贴上吧

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

这三段话其实一点用也没有

FFT是干什么的

FFT在算法竞赛中就有一个用途:加速多项式乘法(暴言)

简单来说,形如 a 0 X 0 + a 1 X 1 + a 2 X 2 + ⋯ + a n X n a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^n a0X0+a1X1+a2X2++anXn 的代数表达式叫做多项式,可以记作 f ( X ) = a 0 X 0 + a 1 X 1 + a 2 X 2 + ⋯ + a n X n f(X)=a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^n f(X)=a0X0+a1X1+a2X2++anXn,其中 a 0 , a 1 , ⋯ , a n a_0,a_1,⋯,a_n a0,a1,,an叫做多项式的系数 X X X是一个不定元(就是不可以合并),不表示任何值,不定元在多项式中最大项的次数称作多项式的次数

如果我们当前有两个多项式 f ( X ) , g ( X ) f(X),g(X) f(X),g(X),现在要把他们乘起来(求卷积),最朴素的做法就是

∑ i = 0 2 n − 1 ( ∑ j = 0 i a j ∗ b i − j ) ∗ x i \sum \limits _{i=0}^{2n-1} (\sum \limits _{j=0} ^{i} a_j*b_{i-j})*x^i i=02n1(j=0iajbij)xi

这样的复杂度是 Θ ( n 2 ) \varTheta(n^2) Θ(n2)的,十分不美观,FFT就是要将这个过程优化为 Θ ( n log ⁡ n ) \varTheta(n \log n) Θ(nlogn)

前置技能

多项式

见上文

复数

复数形如 a + b i a+bi a+bi,其中 i = − 1 i=\sqrt {-1} i=1

a a a叫作复数的实部, b i bi bi叫做复数的虚部

复数 ( a 1 + b 1 i ) ∗ ( a 2 + b 2 i ) (a_1+b_1i)*(a_2+b_2i) (a1+b1i)(a2+b2i)相乘的值,即 a 1 a 2 − b 1 b 2 + ( a 1 b 2 + a 2 b 1 ) i a_1a_2-b_1b_2+(a_1b_2+a_2b_1)i a1a2b1b2+(a1b2+a2b1)i也是一个复数,同时我们也得到了复数的乘法法则

复数 c + d i c+di c+di可以用这种方式表示出来

(c+di)

复数乘法的在复平面中表现为辐角相加,模长相乘

单位根

复数 w w w满足 w n = 1 w^n=1 wn=1称作 w w w n n n次单位根,下图包含了所有的 8 8 8次单位根(图中圆的半径是1)

8次单位根

同样的,下图是所有的4次单位根

4次单位根

聪明的你也许已经发现了单位根的些许性质,即

w 2 n 2 m = w n m w_{2n}^{2m}=w_{n}^{m} w2n2m=wnm

w n m = − w n m + n 2 w_n^m=-w_n^{m+\frac{n}{2}} wnm=wnm+2n

这两个要记住,一会很有用

多项式的系数表达法

我们有多项式 f ( X ) = a 0 X 0 + a 1 X 1 + a 2 X 2 + ⋯ + a n X n f(X)=a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^n f(X)=a0X0+a1X1+a2X2++anXn,令 a ⃗ = ( a 0 , a 1 , a 2 , . . . , a n ) \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)的系数表示法

在系数表示法下,计算多项式乘法是 Θ ( n 2 ) \varTheta (n^2) Θ(n2)

多项式的点值表达法

任取 n + 1 n+1 n+1互不相同 S = { p 1 , p 2 , . . . , p n + 1 } S=\{p_1,p_2,...,p_{n+1}\} S={p1,p2,...,pn+1},对 f ( X ) f(X) f(X)分别求值得到 f ( p 1 ) , f ( p 2 ) , . . . , f ( p n + 1 ) f(p_1),f(p_2),...,f(p_{n+1}) f(p1),f(p2),...,f(pn+1),那么称 A ( X ) = { ( p 1 , f ( p 1 ) ) , ( p 2 , f ( p 2 ) ) , . . . , ( p n + 1 , f ( p n + 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) S S S下的点值表示法

可以把多项式想象成一个 n n n次函数,点值表示法就是取 S S S下每一个横坐标时对应的点,因为 n n n次函数可以由 n + 1 n+1 n+1个点确定下来(可以将每一个点列一个 n n n次方程),所以 n n n维点值与 n n n维系数一一对应

更重要的一点,点值表示法下的乘法运算获得了简化

两个多项式 P P P, Q Q Q分别取点 ( x , y 1 ) (x,y_1) (x,y1) ( x , y 2 ) (x,y_2) (x,y2) P ∗ Q P*Q PQ就会取到点 ( x , y 1 ∗ y 2 ) (x,y_1*y_2) (x,y1y2)

C = P ∗ Q C=P*Q C=PQ,因为 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 ) = y 1 ∗ y 2 C(x)=y_1*y_2 C(x)=y1y2

所以在点值表示法下,计算多项式乘法是 Θ ( n ) \varTheta(n) Θ(n)

FFT的具体过程

FFT就是将系数表示法转化成点值表示法相乘,再由点值表示法转化为系数表示法的过程,第一个过程叫做求值(DFT),第二个过程叫做插值(IDFT)

求值

还记得我们之前提到的单位根吗?再回顾一下:

w 2 n 2 m = w n m w_{2n}^{2m}=w_{n}^{m} w2n2m=wnm

w n m = − w n m + n 2 w_n^m=-w_n^{m+\frac{n}{2}} wnm=wnm+2n

想要求出一个多项式的点值表示法,需要选出 n + 1 n+1 n+1个数分别带入到多项式里面,带入一个数的复杂度是 Θ ( n ) \varTheta(n) Θ(n)的,那么总复杂度就是 Θ ( n 2 ) \varTheta(n^2) Θ(n2)的,因为单位根有上面两个优美的性质,所以我们尝试可以取 n n n次单位根组成 S S S,看看能不能加速我们的运算

A 0 ( X ) A_0(X) A0(X) A ( X ) A(X) A(X)偶次项的和,设 A 1 ( X ) A_1(X) A1(X) A ( X ) A(X) A(X)奇次项的和,即

A 0 ( X ) = a 0 x 0 + a 2 x 1 + . . . + a n − 1 x n / 2 A_0(X)=a_0x^0+a_2x^1+...+a_{n-1}x^{n/2} A0(X)=a0x0+a2x1+...+an1xn/2

A 1 ( X ) = a 1 x 0 + a 3 x 1 + . . . + a n − 2 x n / 2 A_1(X)=a_1x^0+a_3x^1+...+a_{n-2}x^{n/2} A1(X)=a1x0+a3x1+...+an2xn/2

因为 A ( w n m ) = a 0 w n 0 + a 1 w n m + a 2 w n 2 m + a 3 w n 3 m + . . . + a n − 1 w n ( n − 1 ) ∗ m + a n w n n m A(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+...+an1wn(n1)m+anwnnm

所以有

​在这里插入图片描述

在这里插入图片描述

也就是说,只要有了 A 0 ( X ) A_0(X) A0(X) A 1 ( X ) A_1(X) A1(X)的点值表示,就能在 Θ ( n ) \varTheta(n) Θ(n)时间算出 A ( X ) A(X) A(X)的点值表示,对于当前层确定的位置 i i i,就可以用下一层的两个值更新当前的值,我们称这个操作为“蝴蝶变换”

因为这个过程一定要求每层都可以分成两大小相等的部分,所以多项式最高次项一定是 2 p 2^p 2p( p ∈ N p\in N pN)次方,如果不是的话,直接在最高次项补零就可以啦

于是我们有了递归的写法:

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的过程,在第 i i i次求解中,我们将所有元素二进制 i i i位为 0 0 0的放在了左面, i i i位为 1 1 1的放在了右面,事实上,每个元素最终到的是他二进制颠倒过来的位置

这里写图片描述

再拿一张别人的图

这里写图片描述

迭代写法:

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;
}

插值

有人证出来插值只要将所有 w n m w_n^m wnm换成 w n m + n 2 w_n^{m+\frac{n}{2}} wnm+2n,也就是所有的虚部取相反数,再将最终结果除以 l e n len len,就是插值的过程了

究竟为什么?我觉得人生一定要有点遗憾可以参考这里,我就不多说了

支持插值的迭代写法:

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了!

生成函数

A A A a i a_i ai个价值为 A i A_i Ai的物品,小 B B B b i b_i bi个价值为 A i A_i Ai的物品,求用两个组成价值为 c i c_i ci的方案数

生成函数可以解决上面的这个问题,构造两个多项式,第 i i i项表示价值为 i i i的物品有多少个,对两个人分别构造,乘在一起的多项式就代表所有的方案数了

NTT——快速数论变换

当FFT需要取模时,复数的存在就特别尴尬了,有没有什么东西可以代替单位根呢?

原根性质:假设 g ​ g​ g是素数 p ​ p​ p的一个原根,则 g 1 , g 2 , . . . , g ( p − 1 ) ​ g^1,g^2,...,g^{(p-1)}​ g1,g2,...,g(p1)在模 p ​ p​ p意义下两两不同,且结果恰好为 1 ​ 1​ 1~ p − 1 ​ p-1​ p1

w n n ≡ g ( p − 1 ) ≡ 1 ( m o d p ) w_n^n≡g^{(p-1)}≡1(mod\ p) wnng(p1)1(mod p)(中间那个是费马小定理)

所以可以把 g ( p − 1 ) / n g^{(p-1)/n} g(p1)/n看成 w n 1 w_n^1 wn1的等价。

但这种质数必须是NTT质数(费马质数),即 ( p − 1 ) (p-1) (p1)有超过序列长度的2的正整数幂因子的质数

具体证明还是看这里吧

参考资料:
Pick’s blog
Fenghr的博客
Miskcoo’s Space


http://chatgpt.dhexx.cn/article/pKSMMM80.shtml

相关文章

怎么计算网站流量?

如何计算网站的流量呢&#xff1f;在这华仔给大家分享一个如何计算流量的算法&#xff1a; 举个栗子&#xff1a;1G1024M&#xff0c;10G就是101024M10240M. 一个1M的文件被下载1000次的流量约为1G;被下载10000次的流量约为10G. 假如你每月的网站流量为10G&#xff0c;那10G的流…

CDN流量是什么,怎么计算?

CDN流量通俗来讲就是使用CDN加速时&#xff0c;网络加速会产生一个数据使用量&#xff0c;到达某一个时段&#xff0c;统计出这个时段使用的量&#xff0c;也就是CDN流量&#xff0c;和网站流量的使用很像。 随着互联网的发展&#xff0c;用户在使用网络时对网站的浏览速度和…

流量计算器

为所做的工作处于初始开发的阶段&#xff0c;所以数据一直在变化&#xff0c;导致1s的流量大小一直也在变化。每次都需要手动根据新的参数进行计算&#xff0c;真的好烦。 所以呢&#xff0c;作为一个程序员&#xff0c;能让程序做的事情&#xff0c;自己就不要动手了呀&#x…

Flink1.14相关-3.数据流量计算

1. 前言 使用环境 Flink1.14.6Centos7.9Java8 Flink1.14.6安装部署测试参考&#xff1a;参考链接 2 .数据采集 2.1 信号监听 需要监听文件夹下的新文件产生&#xff0c;并且数据库中的值未更新时才发送消息通知后续模块开始采集数据。【后面需要重新搭建监听部分和判断部分…

适应多场景的客流量统计-人流量统计算法

在商场、展厅、景区等受人流量影响较大的场所&#xff0c;流量统计算法可以快速获取流量数据和动态趋势&#xff0c;辅助评估店铺和部分活动的效果&#xff0c;帮助商业决策。另外&#xff0c;在地铁站、火车站、机场等公共场所。实时检测人数可以及时预警高密度人群&#xff0…

网站的服务器带宽计算,服务器带宽和流量计算方式

服务器带宽和流量计算方式。网站流量和带宽该怎么计算&#xff0c;给一些参照数据信息&#xff0c;由于不清楚流量和带宽是怎么计算的&#xff0c;因此不清楚用多宽的带宽更有效、更节省资产! 网站服务器上最好是安装流量统计手机软件(强烈推荐应用DUMeter)&#xff0c;假如流量…

博途1200/1500PLC累计流量计算FB(SCL算法详解+优化)

在编写这篇博客之前其实已经写过一篇SMART S7-200PLC的流量累计的应用文章,由于很多朋友咨询博途PLC下的流量累计实现算法。今天我们以博途PLC的开发环境为例详细讲解算法的实现原理和注意事项同时给出算法的优化方法。受水平和能力所限,文中难免出现错误和不足之处,诚恳的欢…

阿里云轻量应用服务器流量计算方法

阿里云轻量应用服务器套餐有峰值带宽限制每月流量的&#xff0c;还有固定带宽的&#xff0c;阿里云轻量应用服务器流量是怎么计算的&#xff1f;阿里云轻量应用服务器来说说不同套餐下轻量服务器流量计算方法&#xff1a; 轻量应用服务器带宽套餐和流量套餐 阿里云轻量应用服…

物联网GPRS模块流量计算

物联网GPRS模块流量计算 MQTT(消息队列遥测传输) 是ISO 标准下一个基于TCP/IP的消息发布/订阅传输协议。 一、TCP消耗流量计算 以太网数据包结构&#xff1a; 以太网首部 IP首部 TCP首部 APPL首部 用户数据 以太网尾部 以太网首部为14个字节 IP首部为20个字节 TCP首部…

恒容容器放气的瞬时流量的计算

有时候&#xff0c;你会遇到一个问题&#xff0c;该问题的描述如下&#xff1a; 你有一个已知体积的容器&#xff0c;设容器体积为&#xff0c;里面装有一定压力(初始压力)的气体&#xff0c;如空气或氢气等&#xff0c;设初始压力为&#xff0c;容器出口连接着一个阀门开关&am…

阿里云服务器公网带宽流量是怎么计算的?

阿里云服务器流量如何计算&#xff1f;云服务器出流量和入流量都要计算吗&#xff1f;不是&#xff0c;只计算云服务器公网出方向流量&#xff0c;阿里云服务器内网流量和公网出方向流量都是免费的&#xff0c;护云盾来详细说下阿里云服务器流量计算及流量收费说明&#xff1a;…

计算机网络-流量强度

若R链路带宽&#xff08;链路宽度&#xff09;&#xff0c;L分组长度&#xff08;一个分组的大小&#xff09;&#xff0c;a分组到达队列的平均速率&#xff08;分组数量&#xff09;&#xff0c;流量强度公式 &#xff1a;I La/R 举例理解&#xff1a; 当汽车排队从关卡上桥…

腾讯云服务器带宽计费模式按流量内网、外网出入流量计算说明

腾讯云服务器公网带宽计费模式按使用流量计费&#xff0c;云服务器对内或对外产生的流量如何计算&#xff1f;云服务器出方向&#xff08;下行流量&#xff09;和入方向&#xff08;上行流量&#xff09;怎么计算&#xff1f;腾讯云百科来详细说下腾讯云服务器按使用流量计算说…

计算视频流量

码率也可以叫比特率&#xff0c;就是一种音乐每秒播放的数据量&#xff0c;单位用bit表示&#xff0c;也就是二进制位。 bps就是比特率。b就是比特&#xff08;bit&#xff09;&#xff0c;s就是秒&#xff08;second&#xff09;&#xff0c;p就是每&#xff08;per&#xff0…

电缆载流量计算对照表

10.6/1KV聚氯乙烯绝缘电力电缆载流量 常用型号VV22、VLV22聚氯乙烯绝缘钢带铠装聚氯乙烯护套电力电缆载流量 注&#xff1a;以上电缆载流量计算条件 1. 线芯长期工作温度&#xff1a;70℃&#xff1b; 2. 环境温度&#xff1a;25℃ &#xff1b; 3. 埋地深度&#xff1a;10…

明渠如何快速估算水流量(明渠流量计算)

明渠流量估算一般采用速度面积法估算&#xff0c;如果你有流速仪可以测量渠道的流速&#xff0c;如果没有&#xff0c;可以通过漂浮物与秒表来估算&#xff0c;漂浮物容易受到风速影响&#xff0c;风大了是不行的&#xff0c;比如通过一个漂浮物&#xff0c;在时间t内漂浮的距离…

IPCam网络摄像头

文章目录 软件安装及编译环境搭建及代码获取1、于VirtualBoxVM安装Ubuntu2、Ubuntu开机设定3、MobaXterm安装及开发板连接4、套件安装以及SDK编译5、 如何获取代码6、如何更新代码 IPCam网络摄像头7、如何编译IPcam8、配置板端资源以及环境9、参数配置及运行效果9.1参数配置文件…

如何低成本化实时网络摄像头监控直播?

大众直播时代&#xff0c;处处有直播&#xff0c;直播已经在方方面面改变着人们的生活和工作。随着网络直播应用生态的越发完善&#xff0c;你会发现&#xff0c;很多传统监控升级为互联网直播的应用越来越多。那么&#xff0c;如何将常规监控摄像头转为互联网直播&#xff1f;…

【解决方案】5G时代浪潮来袭,EasyNVR助力5G厂区视频监控安防采集可视化展示

智慧工厂被认为是5G技术的重要应用场景之一&#xff0c;利用5G网络将生产设备无缝连接&#xff0c;并进一步打通设计、采购、仓储、物流等环节&#xff0c;满足工业环境下设备互联和远程交互应用需求。TSINGSEE青犀视频面向工厂智能化升级需求&#xff0c;推出5G智慧工厂方案&a…

多摄像机网络智能视频监控系统设计与实现

1、智能视频监控系统项目背景介绍 2、系统的需求分析 在多摄像机视频监控系统中,通过前端结构化信息提取设备获取单个摄像机的结构化信息以后,需要将数据传输到云端服务器当中进行统一的存储和管理,从而为再识别等服务提供可靠的数据来源,满足更加丰富的用户需求。 本章将…