快速傅里叶变换FFT学习小记

FFT学得还是有点模糊,原理那些基本还是算有所理解了吧,不过自己推这个推不动。

看的资料主要有这两个:

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

https://www.zybuluo.com/397915842/note/37965

这儿简单做做笔记。

多项式点值表示

首先$FFT$用来快速计算两个多项式的乘积。

一个$n$次多项式(最高次为$n$),可以用系数表示法表示和点值表示法表示。

  • 系数表示法:$A(x)=\sum_{i=0}^{n}c_ix^i$,可以知道用系数表示法进行多项式乘法时间复杂度是$\Theta (n^2)$。
  • 点值表示法:用$n+1$个点$(x_k,y_k),k=0\dots n$,其中$x_k$互不相同,$y_k=A(x_k)$。

这$n+1$个点是可以唯一表示一个$n$次多项式的,因为:$$\begin{bmatrix} 1 & x_0 & x_0^2 & ... & x_0^n \\ 1 & x_1 & x_1^2 & ... & x_1^n \\\vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & ... & x_n^n \\\end{bmatrix}\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\c_n \\\end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\y_n \\\end{bmatrix}$$

前面那个方阵是范德蒙德矩阵,其行列式可以知道是非0,所以该方阵存在逆矩阵,方程有解且有唯一解。

现在考虑两个n次多项式$A(x)$和$B(x)$,其乘积的多项式结果是$C(x)$,那么$C(x)$的次数为$2n$。

不妨给$A$和$B$取$2n+1$个点,这样求得两个多项式的点值表示为$(x_k,A(x_k))、(x_k,B(y_k)),k=0\dots 2n$,而$C(x)$的点值表示在$\Theta (n)$就能直接求出,即为$(x_k,A(x_k)B(y_k))$。

这样进行多项式乘法运算,原本用系数表示多项式需要$\Theta (n^2)$的时间复杂度,而点值则只需要$\Theta (n)$。

不过将一个多项式转化成点值的形式,需要枚举$2n$个点,然后对于每一个点$\Theta (n)$求得多项式的值,这样时间复杂度是$\Theta (n^2)$的;另外,求得结果后也是一个点值的表示,需要还原回系数形式,直接解上面矩阵表示的方程,这要立方的时间复杂度。

这时就是$FFT$做的事了。前面转化为点值形式称为$DFT$,后面逆转化回系数形式称为$IDFT$。

DFT

对于多项式$A(x)$,比如我们要取$n$个点$(x_k,y_k)$,$x_k$选取的是$n$次单位根(满足$z^n=1$的复数解$z$),即:$$e^{\frac{2\pi ki}{n}}, k = 0\cdots n - 1, i=\sqrt{-1}$$

另外,计算这个可以用欧拉公式:$$e^{\theta i}=\cos\theta + i\sin\theta$$

记$\omega_n=e^{\frac{2\pi i}{n}}$,那么点值表示即为$$(\omega_n^0, A(\omega_n^0)), (\omega_n^1, A(\omega_n^1)),\cdots ,(\omega_n^{n-1}, A(\omega_n^{n-1}))$$

现在如何求$A(\omega_n^0), A(\omega_n^1), \cdots, A(\omega_n^{n-1})$?

看下面:$$A(x) = c_0 + c_1x + c_2x^2 + c_3x^3 + \cdots + c_nx^{n-1}\\A^{[0]}(x) = c_0 + c_2x + c_4x^2 + \cdots + c_{n - 2}x^{n/2 - 1} \\A^{[1]}(x) = c_1 + c_3x + c_5x^2 + \cdots + c_{n - 1}x^{n/2 - 1}$$

然后,可以发现$$A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$$

而这边是将$\omega_n^k,k=0\cdots n-1$代入$x$,其中$x^2$也就是$(\omega_n^k)^2$可以化简:$$(\omega_n^k)^2=\left(e^{\frac{2\pi ki}{n}}\right)^2=e^{\frac{2\pi ki}{n/2}}=\omega_{\frac{n}{2}}^k$$而$\omega_{\frac{n}{2}}^k$复数根有$n/2$个。

那么可以看到,把$A(x)$分成两个子问题,而这两个子问题选取的点也缩减了一半。这样时间复杂度就是$\Theta (nlogn)$了。于是可以用递归去实现这个分治的算法。另外一开始$n$拓展成$2$的次幂,多出来的补$0$即可。

不过,事实上可以改用迭代实现,而不用递归。画出递归的每一次划分。。$此处没图$。。

然后对比最初和最后一层各个下标的二进制,可以发现相当于是每一个下标的二进制反转了。

那么可以一开始直接反转,从最后一层开始,一层层一个个迭代上去去求得结果。

另外,有个现成的$Rader$雷德算法可以直接求得反转结果。

IDFT

上面就在$\Theta (nlogn)$下将系数表示转化成点值表示;而现在需要逆过来,把点值表示还原成所需要的结果,系数表示。$$ \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots &(\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix}$$

其实就是求$c_k$。

而那个范德蒙德矩阵的逆也有结论的,就是:$$\frac{1}{n}\begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix}$$

证明过程。。两个矩阵相乘。。我就不搬运了= =。。

这样的话等式两边同时乘上方阵的逆:$$\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix} = \frac{1}{n} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix}$$

那么可以发现,这可以用$DFT$一样的过程去求了,不同的是选的点是$\omega_n^{-k}$,最后结果除以$n$。

总之

总的过程大概就是:

  1. 进行一些预处理,把次数拓展成2的次幂,高位多出来的设置成0;
  2. 用$DFT$在$\Theta (nlogn)$分别将两个多项式$A(x)$和$B(x)$转化成点值的形式;
  3. 利用点值表示在$\Theta (n)$直接就能求得$A(x)\times B(x)$的点值表示;
  4. 用$IDFT$在$\Theta (nlogn)$将乘积结果的点值表示还原成系数表示。

代码的实现照着别人的,稍微改了下点,写了一遍:

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define MAXN 133333
const double PI=acos(-1.0);

struct Complex{
	double real,imag;
	Complex(){}
	Complex(double _real,double _imag):real(_real),imag(_imag){}
	Complex operator-(const Complex &cp) const{
		return Complex(this->real-cp.real,this->imag-cp.imag);
	}
	Complex operator+(const Complex &cp) const{
		return Complex(this->real+cp.real,this->imag+cp.imag);
	}
	Complex operator*(const Complex &cp) const{
		return Complex(this->real*cp.real-this->imag*cp.imag,this->real*cp.imag+this->imag*cp.real);
	}
	void setValue(double _real=0,double _imag=0){
		real=_real; imag=_imag;
	}
};

int len;
Complex wn[MAXN],wn_anti[MAXN];

void FFT(Complex y[],int op){
	// Rader, 位逆序置换
	for(int i=1,j=len>>1,k; i<len-1; ++i){
		if(i<j) swap(y[i],y[j]);
		k=len>>1;
		while(j>=k){
			j-=k;
			k>>=1;
		}
		if(j<k) j+=k;
	}
	// h=1,Wn^0=1
	for(int h=2; h<=len; h<<=1){
		// Wn = e^(2*PI*i/n),如果是插值Wn = e^(-2*PI*i/n),i虚数单位
		Complex Wn = (op==1 ? wn[h] : wn_anti[h]);
		for(int i=0; i<len; i+=h){
			Complex W(1,0);
			for(int j=i; j<i+(h>>1); ++j){
				Complex u=y[j],t=W*y[j+(h>>1)];
				y[j]=u+t;
				y[j+(h>>1)]=u-t;
				W=W*Wn;
			}
		}
	}
	if(op==-1){
		for(int i=0; i<len; ++i) y[i].real/=len;
	}
}
void Convolution(Complex A[],Complex B[],int n){
	// 初始化
	for(len=1; len<(n<<1); len<<=1);
	for(int i=n; i<len; ++i){
		A[i].setValue();
		B[i].setValue();
	}
	// e^(θi) = cosθ+isinθ -> Wn = cos(2*PI/n)+isin(2*PI/n)
	for(int i=0; i<=len; ++i){ // 预处理
		wn[i].setValue(cos(2.0*PI/i),sin(2.0*PI/i));
		wn_anti[i].setValue(wn[i].real,-wn[i].imag);
	}
	FFT(A,1); FFT(B,1);
	for(int i=0; i<len; ++i){
		A[i]=A[i]*B[i];
	}
	FFT(A,-1);
}
时间: 2024-10-04 05:16:31

快速傅里叶变换FFT学习小记的相关文章

【笔记篇】(理论向)快速傅里叶变换(FFT)学习笔记w

现在真是一碰电脑就很颓废啊... 于是早晨把电脑锁上然后在旁边啃了一节课多的算导, 把FFT的基本原理整明白了.. 但是我并不觉得自己能讲明白... Fast Fourier Transformation, 快速傅里叶变换, 是DFT(Discrete Fourier Transform, 离散傅里叶变换)的快速实现版本. 据说在信号处理领域广泛的应用, 而且在OI中也有广泛的应用(比如SDOI2017 R2至少考了两道), 所以有必要学习一波.. 划重点: 其实学习FFT最好的教材是<算法导论

[学习笔记] 多项式与快速傅里叶变换(FFT)基础

引入 可能有不少OIer都知道FFT这个神奇的算法, 通过一系列玄学的变化就可以在 $O(nlog(n))$ 的总时间复杂度内计算出两个向量的卷积(或者多项式乘法/高精度乘法), 而代码量却非常小. 博主一年半前曾经因COGS的一道叫做"神秘的常数 $\pi$"的题目而去学习过FFT, 但是基本就是照着板子打打完并不知道自己在写些什么鬼畜的东西OwO 不过...博主这几天突然照着算法导论自己看了一遍发现自己似乎突然意识到了什么OwO然后就打了一道板子题还1A了OwO再加上午考试差点AK

快速傅里叶变换FFT

快速傅里叶变换FFT DFT是信号分析与处理中的一种重要变换.但直接计算DFT的计算量与变换区间长度N的平方成正比,当N较大时,计算量太大,直接用DFT算法进行谱分析和信号的实时处理是不切实际的. 1.直接计算DFT 长度为N的有限长序列x(n)的DFT为: 2.减少运算量的思路和方法 思路:N点DFT的复乘次数等于N2.把N点DFT分解为几个较短的DFT,可使乘法次数大大减少.另外,旋转因子WmN具有周期性和对称性. (考虑x(n)为复数序列的一般情况,对某一个k值,直接按上式计算X(k)值需

浅谈范德蒙德(Vandermonde)方阵的逆矩阵的求法以及快速傅里叶变换(FFT)中IDFT的原理

浅谈范德蒙德(Vandermonde)方阵的逆矩阵与拉格朗日(Lagrange)插值的关系以及快速傅里叶变换(FFT)中IDFT的原理 只要稍微看过一点线性代数的应该都知道范德蒙德行列式. \[V(x_0,x_1,\cdots ,x_{n-1})=\begin{bmatrix} {1}&{1}&{\cdots}&{1}\{x_{0}}&{x_{1}}&{\cdots}&{x_{n-1}}\{x_{0}^2}&{x_{1}^2}&{\cdots

【转】快速傅里叶变换(FFT)详解

目录 前言 多项式 系数表示法 点值表示法 复数 向量 圆的弧度制 平行四边形定则 复数 运算法则 单位根 单位根的性质 快速傅里叶变换 快速傅里叶逆变换 理论总结 递归实现 迭代实现 本文只讨论FFT在信息学奥赛中的应用 文中内容均为个人理解,如有错误请指出,不胜感激 回到顶部 前言 先解释几个比较容易混淆的缩写吧 DFT:离散傅里叶变换->O(n2)O(n2)计算多项式乘法 FFT:快速傅里叶变换->O(n?log(n)O(n?log?(n)计算多项式乘法 FNTT/NTT:快速傅里叶变换

图像傅里叶变换(快速傅里叶变换FFT)

学习DIP第7天,图像傅里叶变换 习惯性,开篇废话 今天公司的网不知怎么了,死活打不开CSDN,公司有100多架客机,也有极限速度60kb/s的网速,还有3K的工资. 图像FFT 上篇已经介绍了关于2D FFT的相关知识,这篇只介绍在图像中的应用,对于一幅图像,做二维FFT后,即可得到其傅里叶变换,傅里叶变换后是二维复数矩阵,因为二维数组,如果是实数,是可以通过变换到0~255通过灰度图像显示出来,而变换结果是复数,所以我们通过显示其幅度,即复数的模,来显示傅里叶谱(幅度谱),不废话,上图: 原

几种快速傅里叶变换(FFT)的C++实现

DFT的的正变换和反变换分别为(1)和(2)式.假设有N个数据,则计算一个频率点需要N次复数乘法和N-1次复数加法,整个DFT需要N*N次复数乘 法和N(N-1)次复数加法:由于一次的复数乘法需要进行4次的实数乘法和2次的复数加法,一次的复数加法需要两次的实数加法,因此整个DFT需要 4*N*N次的实数乘法和2*N(N-1)+2*N*N≍4*N*N次的复数加法.当N比较大时,所需的计算工作量相当大,例如N=1024时大约需要 400万次乘法运算,对于实时信号处理来说,将对计算设备提出非常苛刻的要

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

//有多少项N取多至少8~10倍 #include <algorithm> #include <cmath> #include <complex> #include <cstdlib> using namespace std; typedef complex<double> cn; const int N=1e7+1; const double Pi=3.1415926535;//Pi通常取10位即可(随意) cn a[N],b[N]; int

快速傅里叶变换(FFT):COGS 2216. 你猜是不是KMP

2216. 你猜是不是KMP ★★★☆   输入文件:guess.in   输出文件:guess.out   简单对比时间限制:1 s   内存限制:256 MB [题目描述] XX在玩两个串的游戏.首先,他拿出了两个字符串 S 和 T,XX想知道 T在 S 中出现了几次,分别在哪些位置出现.注意 T 中可能有“?”字符,这个字符可以匹配任何字符. [输入格式] 两行两个字符串,分别代表 S 和 T [输出格式] 第一行一个正整数 k,表示 T 在 S 中出现了几次. 接下来 k 行正整数, 分