首页 > 其他分享 >FFT

FFT

时间:2024-09-13 21:13:13浏览次数:1  
标签:right frac FFT 点值 omega fushu left

FFT

简介

用于求卷积(\(a,b\) 已知):

\[\sum_{i=0}^n a_ib_{n-i} \]

或者多项式乘法(\(A(x),B(x)\) 已知):

\[C(x)=A(x)B(x) \]

\(A(x)=\sum_{i=0}^{n} a_i x^i\\ B(x)=\sum_{i=0}^{m} b_i x^i\)

可见 \(C(x)\) 是 \(n+m\) 次多项式。

如果我们把卷积的 \(a_i,b_i\) 看成多项式的系数,卷积就变成求:

\[[x^n] C(x)=A(x)B(x) \]

其中 \([x^n]\) 表示 \(x^n\) 的系数。

求卷积或者多项式乘法的时间复杂度是 \(O(n^2)\) 的。使用 FFT 可以做到 \(O(n\log n)\)。

大体思路

设 \(C(x)\) 的项数为 \(n\)(不是次数),若 \(A,B\) 不足 \(n\) 项就补系数 \(0\)。

显然这个过程很难优化,我们从另一个角度去想。

对于一个多项式,求其在 \(x\) 处的值的时间复杂度是 \(O(n)\) 的,我们把这个操作叫做点值(DFT)。

\(n\) 个点可以唯一确定一个 \(n\) 项多项式(即 \(n-1\) 次多项式),证明不显然但是略过我不会。而据说拉格朗日插值法给出了一种 \(O(n^2)\) 的求出这个多项式(即求出它的系数)的方法我不会,这个由点值求系数的过程叫做插值(IDFT)。

因此求多项式乘法,可以变成求任意 \(n\) 个点在两个多项式 \(A,B\) 的点值,然后由 \(C(x)=A(x)B(x) \ O(n)\) 求出多项式 \(C\) 在这 \(n\) 个点的点值,然后做一次插值求出 \(C\) 的系数。当然这个也是 \(O(n^2)\) 的,但是聪明的傅里叶给出了一种基于单位根的特殊性质的分治方法求 DFT 和 IDFT,成为快速傅里叶变换(FFT)。

单位根

写作 \(\omega_n^k\),读作 \(n\) 次单位根的 \(k\) 次方。

\(\omega_n^k\) 是在复数域上的向量,形如在复平面上分成 \(n\) 等分,其中 \(\omega_n^0=\omega_n^n=1\)。按逆时针分别为 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\)。

后面为了方便,我们会假设 \(n=2^k\)。

有几个重要的性质。

  • \(\omega_n^n=1\) 正确性显然
  • \(\omega_{an}^{ak}=\omega_n^k\) 在平面上想象一下,显然正确
  • \(\omega_n^{k-\frac{n}{2}}=-\omega_n^k\) 相当于把向量转 \(180^。\)。

点值

基于这些性质,我们求 \(A,B\) 在 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\) 处的点值。

以求 \(A(x)\) 为例,我们要求所有 \(A(\omega_n^k),0\le k<n\)。

我们把 \(A\) 按奇偶分为两部分:

\[A_1(x)=a_0+a_2 x+a_4 x^2+\dots a_{n-2} x^{\frac{n-2}{2}}\\ A_2(x)=a_1+a_3x^2+a_5x^4+\dots a_{n-1}x^{\frac{n-2}{2}}\]

因此有:

\[A(x)=A_1(x^2)+xA_2(x^2) \]

代入 \(x=\omega_n^k\),有:

\[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^k A_2(\omega_n^{2k})\\ A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^k A_2(\omega_{\frac{n}{2}}^k)\]

我们先求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值。这个范围是缩小了一半的。然后把它们相加就得到了 \(A\) 在 \(k< \frac{n}{2}\) 的点值。

然后我们求剩下一半的 \(A_1,A_2\) 的点值。仍然设 \(k< \frac{n}{2}\) 发现:

\[A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}} A_2(\omega_n^{2k+n})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^{k} A_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^{k} A_2(\omega_{\frac{n}{2}}^k)\]

因此你发现,这俩十分地相似,因此你求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值之后,可以直接 \(O(n)\) 求出 \(A\) 的 \(n\) 个点值了。然后就这样分治下去,点值时间复杂度为 \(O(n\log n)\)。

插值

求插值的过程是类似的,改几个参数就行了。我还没搞懂是怎么推出来的

然后我们就这么求两次点值,再求一次插值就做完了。

改进

code

然后你会发现过不了板子……

因为这个递归的过程常数很大,假设 FFT 的常数本身就大,然后就超时了。

考虑从下往上递推分治的过程。

\[\left(a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7}\right)\\ \left(a_{0}, a_{2}, a_{4}, a_{6}\right)\left(a_{1}, a_{3}, a_{5}, a_{7}\right)\\ \left(a_{0}, a_{4}\right)\left(a_{2}, a_{6}\right)\left(a_{1}, a_{5}\right)\left(a_{3}, a_{7}\right)\\ \left(a_{0}\right)\left(a_{4}\right)\left(a_{2}\right)\left(a_{6}\right)\left(a_{1}\right)\left(a_{5}\right)\left(a_{3}\right)\left(a_{7}\right) \]

然后你惊奇地发现最下层的顺序是 \(000,100,010,110,001,101,011,111\),刚好是 \(0\sim 7\) 的二进制的 reverse。

然后就有了如下板子:

Code

#include<bits/stdc++.h>
#define sf scanf
#define pf printf
using namespace std;
typedef long long ll;
const int N=4e6+3;
int n,m;
struct fushu{
	double x,y;
	fushu (double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
fushu operator + (fushu a,fushu b){ return (fushu){a.x+b.x,a.y+b.y}; };
fushu operator - (fushu a,fushu b){ return (fushu){a.x-b.x,a.y-b.y}; };
fushu operator * (fushu a,fushu b){ return (fushu){a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y}; }
int len;
const double pi=acos(-1.0);
int re[N];
void FFT(fushu *c,int type){
	for(int i=0;i<(1<<len);i++){
		if(i<re[i]) swap(c[i],c[re[i]]);
	}
	for(int mid=1;mid<(1<<len);mid<<=1){
		fushu Wn( cos(pi/mid) , type*sin(pi/mid) );
		for(int r=mid<<1,j=0;j<(1<<len);j+=r){
			fushu w(1,0);
			for(int i=0;i<mid;i++,w=w*Wn){
				fushu x=c[j+i],y=w*c[j+mid+i];
				c[j+i]=x+y;
				c[j+mid+i]=x-y;
			}
		}
	}
}
int main(){
	sf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) sf("%lf",&a[i].x);
	for(int i=0;i<=m;i++) sf("%lf",&b[i].x);
	while((1<<len)<=n+m) len++;
	for(int i=0;i<(1<<len);i++){
		re[i]=(re[i>>1]>>1)|((i&1)<<(len-1));
	}
	FFT(a,1);
	FFT(b,1);
	for(int i=0;i<(1<<len);i++) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=0;i<=n+m;i++){
		pf("%d ",(int)(a[i].x/(1<<len)+0.5));
	}
}

算法缺陷

由于复数域是用浮点数计算的,所以会存在掉精度问题。如果答案是对一个特别的指数取模,如著名的 \(998244353\),可以使用原根代替单位根计算,在剩余系里计算而不是在复数域计算。详见 NTT。

标签:right,frac,FFT,点值,omega,fushu,left
From: https://www.cnblogs.com/liyixin0514/p/18412826

相关文章

  • Cooley-Tukey FFT算法的非递归实现
    一维情况#include<iostream>#include<complex>#include<cmath>constdoublePI=3.14159265358979323846;//交换位置voidswap(std::complex<double>&a,std::complex<double>&b){std::complex<double>temp=a......
  • RFFT:数据与代码已开源,京东推出广告图生成新方法 | ECCV 2024
    论文将多模态可靠反馈网络(RFNet)结合到一个循环生成图片过程中,可以增加可用的广告图片数量。为了进一步提高生产效率,利用RFNet反馈进行创新的一致条件正则化,对扩散模型进行微调(RFFT),显著增加生成图片的可用率,减少了循环生成中的尝试次数,并提供了高效的生产过程,而不牺牲视觉吸引力。......
  • 如何编译FFTW3库:静态库与动态库的编译指南
    目录1.下载并解压FFTW3库2.配置编译选项3.编译并安装库4.验证编译结果5.在项目中使用FFTW3库6.总结FFTW3(FastestFourierTransformintheWest)是一个广泛使用的高性能傅里叶变换库。它支持多种优化,适用于多线程计算和SIMD指令,是处理大规模数据傅里叶变换的理......
  • 在 Ubuntu 中查找库的位置:以 FFTW3 库为例
    目录一、为什么要查找库的位置?二、查找库位置的常用方法三、实践示例四、总结在Ubuntu或其他基于Linux的操作系统中,开发者常常需要查找已安装库的位置,以便进行编译、链接或配置环境变量。本文将详细介绍如何在Ubuntu中查找库的位置,并以常用的FFTW3库为例进行......
  • 音频去噪:使用Python和FFT增强音质
    根据定义,声音去噪是从音频信号中去除不需要的噪音或干扰,以提高其质量和清晰度的过程。这涉及识别和隔离噪音成分(通常以不规则或高频元素为特征),并将其过滤掉,同时保持原始声音的完整性。声音去噪目标是改善聆听体验以及音频分析和处理的准确性。过滤掉噪音对于高保真音频来说......
  • 小学生都能看懂的FFT
    小学生都能看懂的FFT!!作者(小学生)其实学这个学了两个月,但我相信你只要努力,就能成功好的,废话不多说,正片开始FFT章节1:了解FFT是干嘛的oiwiki:FFT支持在O(n......
  • python 实现FFT快速傅立叶变换算法
    FFT快速傅里叶变换介绍FFT(快速傅里叶变换)是计算离散傅里叶变换(DFT)及其逆变换的一种高效算法。DFT是一种将信号从时域转换到频域的数学工具,而FFT通过减少计算量来加速这一过程。FFT的基本思想FFT利用了DFT中的对称性和周期性,通过分而治之的策略将DFT分解为更小的DFT,从而显......
  • Contest5385 - FFT
    ContestA签到题BFFT/NTT快速傅立叶P3803【模板】多项式乘法(FFT)&P1919【模板】高精度乘法|A*BProblem升级版参考资料:炫酷反演魔术-博客-vfleaking的博客题解P3803【【模板】多项式乘法(FFT)】-洛谷专栏FFT总体思路FFT处理循环卷积问题,而卷积问题通用的......
  • STM32H7 HAL库CubeMX 双重ADC模式同步采样详细配置+FFT计算相位差
    前言在电赛备赛期间琢磨了一下ADC同步采样的实现方式,本来是打算直接用AD7606来着,但是搞了半天也没把驱动整出来...考虑到AD7606本身采样率也拉不到太高,于是就花了几天时间把片上ADC配出来了。查资料的时候我发现关于STM32双重ADC模式的资料是真的少,用FFT算两路信号相位差的实例代......
  • MATLAB仿真:数字信号处理用FFT对信号分析
    目录1.实验目的2实验原理3.实验仪器及设备4.实验步骤及内容(1)对以下序列进行谱分析。(2)对以下周期序列进行谱分析。 (3)对模拟周期信号进行谱分析1.实验目的学习用FFT 对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析误差及其原因,以便正确应用 FFT。2......