FFT

我们从一个问题进行讨论 F F T FFT FFT，我们有两个多项式 A ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 , B ( x ) = b 0 + b 1 x + b 2 x 2 + ⋯ + b m − 1 x m − 1 A(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1},\\B(x)=b_0+b_1x+b_2x^2+\dots+b_{m-1}x^{m-1} A(x)=a0+a1x+a2x2+⋯+an−1xn−1,B(x)=b0+b1x+b2x2+⋯+bm−1xm−1我们现在要求 A ( x ) ⋅ B ( x ) A(x)\cdot B(x) A(x)⋅B(x)。朴素算法是显而易见的，时间复杂度是 Θ ( n m ) \Theta(nm) Θ(nm)暴力相乘即可，但这个时间复杂度我们是接受不了的。这时我们可以引入一个点值表示法。

前置知识

点值表示法

多项式通常有两种表示方法，一种是系数表示法即形如 P ( x ) = ∑ i = 0 n − 1 a i x i P(x)=\sum_{i=0}^{n-1}a_ix^i P(x)=∑i=0n−1aixi的表示法，而点值表示法则是把 n n n个互不相同的数 x x x带入多项式 P ( x ) P(x) P(x)中的的值作为 y y y值确定 n n n个点，我们可以证明互不重合的 n n n个点可以唯一确定一个多项式。

单位根

对于 ω n = 1 \omega^n=1 ωn=1的解傅里叶得出的结论是，对于 ω n = 1 \omega^n=1 ωn=1这个方程的复数解为，在复平面上我们对单位圆从正实数轴开始 n n n等分如图，这 n n n个 n n n等分点表示的复数即为原方程的复数解。

我们可以得到以下性质

ω n k = cos ⁡ 2 k π n + i sin ⁡ 2 k π n \omega_n^k=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n} ωnk=cosn2kπ+isinn2kπ
∀ i ≠ j \forall i\neq j ∀i=j， i , j < n i,j<n i,j<n，有 ω n i ≠ ω n j \omega_n^i\neq\omega_n^j ωni=ωnj
ω 2 n 2 k = ω n k \omega_{2n}^{2k}=\omega_n^k ω2n2k=ωnk（折半引理）
ω n k = − ω n k + n 2 \omega_n^k=-\omega_n^{k+\frac n2} ωnk=−ωnk+2n
ω n 0 = ω n n \omega_n^0=\omega_n^n ωn0=ωnn
ω n i ⋅ ω n j = ω n i + j \omega_n^i\cdot\omega_n^j=\omega_n^{i+j} ωni⋅ωnj=ωni+j

FFT推导

正变换

我们要把一个多项式转为点值表示法，朴素转是 Θ ( n 2 ) \Theta(n^2) Θ(n2)，达不到我们的目的。这里我们要把 n n n转为一个 2 2 2的幂的数不够就补位系数补 0 0 0，这样我们可以把多项式分为两个部分
A ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 = ( a 0 + a 2 x 2 + ⋯ + a n − 2 x n − 2 ) + ( a 1 + a 3 x 3 + ⋯ + a n − 1 x n − 1 ) \begin{align}A(x)&=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}\\&=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1+a_3x^3+\dots+a_{n-1}x^{n-1})\end{align} A(x)=a0+a1x+a2x2+⋯+an−1xn−1=(a0+a2x2+⋯+an−2xn−2)+(a1+a3x3+⋯+an−1xn−1)
令 A 1 ( x ) = a 0 + a 2 x + ⋯ + a n − 2 x n 2 − 1 A_1(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac n2-1} A1(x)=a0+a2x+⋯+an−2x2n−1， A 2 ( x ) = a 1 + a 3 x 2 + ⋯ + a n − 1 x n 2 − 1 A_2(x)=a_1+a_3x^2+\dots+a_{n-1}x^{\frac n2-1} A2(x)=a1+a3x2+⋯+an−1x2n−1

可得 A ( x ) = A 1 ( x 2 ) + x A 2 ( x 2 ) A(x)=A_1(x^2)+xA_2(x^2) A(x)=A1(x2)+xA2(x2)

我们可以带入 n n n个 n n n次单位根
k ∈ [ 0 , n 2 − 1 ] , A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) k ∈ [ n 2 − 1 , n − 1 ] , A ( ω n k + n 2 ) = A 1 ( ω n 2 k + n ) + ω n k + n 2 A 2 ( ω n 2 k + n ) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) k\in\bigg[0,\frac n2-1\bigg],A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})=A_1(\omega_{\frac n2}^{k})+\omega_n^kA_2(\omega_{\frac n2}^{k})\\ k\in\bigg[\frac n2-1,n-1\bigg],A(\omega_n^{k+\frac n2})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac n2}A_2(\omega_n^{2k+n})=A_1(\omega_{\frac n2}^{k})-\omega_n^kA_2(\omega_{\frac n2}^{k}) k∈[0,2n−1],A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)=A1(ω2nk)+ωnkA2(ω2nk)k∈[2n−1,n−1],A(ωnk+2n)=A1(ωn2k+n)+ωnk+2nA2(ωn2k+n)=A1(ω2nk)−ωnkA2(ω2nk)
这样我们可以通过分治 log ⁡ 2 n \log_2n log2n层，这样我们可以在 n log ⁡ 2 n n\log_2n nlog2n的时间里求出 n n n个单位根的多项式值，我们把两个多项式的每个单位根对应的多项式值算出来，把同一单位根的两个多项式的值相乘，便可以得到乘积多项式的点值表达式。但在一些数据优秀的 O J OJ OJ上我们并不能通过，所以我们需要一个蝴蝶变换来优化。

蝴蝶变换

我们会发现一个规律，一个项的最后拆分位置恰好为下标的二进制反过来，我们以 8 8 8项的多项式为例，模拟拆分的过程：

位置编号（二进制）	000 000 000	001 001 001	010 010 010	011 011 011	100 100 100	101 101 101	110 110 110	111 111 111
初始序列	x 0 x_0 x0	x 1 x_1 x1	x 2 x_2 x2	x 3 x_3 x3	x 4 x_4 x4	x 5 x_5 x5	x 6 x_6 x6	x 7 x_7 x7
一次拆分后	x 0 x_0 x0	x 2 x_2 x2	x 4 x_4 x4	x 6 x_6 x6	x 1 x_1 x1	x 3 x_3 x3	x 5 x_5 x5	x 7 x_7 x7
两次拆分后	x 0 x_0 x0	x 4 x_4 x4	x 2 x_2 x2	x 6 x_6 x6	x 1 x_1 x1	x 5 x_5 x5	x 3 x_3 x3	x 7 x_7 x7
三次拆分后	x 0 ( 000 ) x_{0(000)} x0(000)	x 4 ( 100 ) x_{4(100)} x4(100)	x 2 ( 010 ) x_{2(010)} x2(010)	x 6 ( 110 ) x_{6(110)} x6(110)	x 1 ( 001 ) x_{1(001)} x1(001)	x 5 ( 101 ) x_{5(101)} x5(101)	x 3 ( 011 ) x_{3(011)} x3(011)	x 7 ( 111 ) x_{7(111)} x7(111)

我们就可以得到一个迭代版的 F F T FFT FFT这样是可以通过所有 O J OJ OJ的洗礼的。

逆变换

正变换我们其实是求了这样一个东西
[ 1 1 1 ⋯ 1 1 ω n 1 ω n 2 ⋯ ω n − 1 1 ω n 2 ω n 4 ⋯ ω 2 n − 2 ⋮ ⋮ ⋮ ⋱ ⋮ 1 ω n − 1 n − 1 ω n − 1 2 n − 2 ⋯ ω n − 1 ( n − 1 ) 2 ] [ a 0 a 1 a 2 ⋮ a n − 1 ] = [ A 0 A 1 A 2 ⋮ a n − 1 ] \begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega_n^1&\omega_n^2&\cdots&\omega^{n-1}\\ 1&\omega_n^2&\omega_n^4&\cdots&\omega^{2n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_{n-1}^{n-1}&\omega_{n-1}^{2n-2}&\cdots&\omega_{n-1}^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_{n-1} \end{bmatrix}= \begin{bmatrix} A_0\\A_1\\A_2\\\vdots\\a_{n-1} \end{bmatrix} ⎣ ⎡111⋮11ωn1ωn2⋮ωn−1n−11ωn2ωn4⋮ωn−12n−2⋯⋯⋯⋱⋯1ωn−1ω2n−2⋮ωn−1(n−1)2⎦ ⎤⎣ ⎡a0a1a2⋮an−1⎦ ⎤=⎣ ⎡A0A1A2⋮an−1⎦ ⎤
由于这个矩阵的元素非常特殊，它的逆矩阵也有特殊的性质，就是每一项取倒数再除以变换的长度就能得到它的逆矩阵。

我们可以得到
1 ω k = ω k − 1 = e − 2 π i k = cos ⁡ ( 2 π k ) + i sin ⁡ ( − 2 π k ) = cos ⁡ ( 2 π k ) − i sin ⁡ ( 2 π k ) \frac{1}{\omega_k}=\omega_k^{-1}=e^{-\frac{2\pi i}{k}}=\cos(\frac{2\pi}{k})+i\sin(-\frac{2\pi}k)=\cos(\frac{2\pi}k)-i\sin(\frac{2\pi}k) ωk1=ωk−1=e−k2πi=cos(k2π)+isin(−k2π)=cos(k2π)−isin(k2π)
只需要我们在正变换的函数里面多传一个系数，表示正/逆变换，在计算的时候乘上系数即可。

code

const double Pi=acos(-1);
struct _cp{double re,im;inline _cp operator + (const _cp b)const{return (_cp){re+b.re,im+b.im};}inline _cp operator - (const _cp b)const{return (_cp){re-b.re,im-b.im};}inline _cp operator * (const _cp b)const{return (_cp){re*b.re-im*b.im,re*b.im+im*b.re};}inline _cp operator / (const _cp b)const{_cp x=(_cp){re,im};double q=b.re*b.re+b.im*b.im;_cp y=(_cp){b.re/q,b.im/q};return x*y;}
};//自己写的复数运算
int rev[5000005];
void FFT(_cp *a,int n,int inv)//inv为逆变换系数
{int bit=0;while((1<<bit)<n)bit++;for(int i=0;i<n;i++){rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));if(i<rev[i])swap(a[i],a[rev[i]]);}for(int mid=1;mid<n;mid*=2){_cp temp=(_cp){cos(Pi/mid),inv*sin(Pi/mid)};for(int i=0;i<n;i+=mid*2){_cp omega=(_cp){1,0};for(int j=0;j<mid;j++,omega=omega*temp){_cp x=a[i+j],y=omega*a[i+j+mid];a[i+j]=x+y,a[i+j+mid]=x-y;}}}
}

FFT快速傅里叶变换相关推荐

【经典算法实现 44】理解二维FFT快速傅里叶变换及 IFFT快速傅里叶逆变换（迭代法和递归法）
[经典算法实现 44]理解二维FFT快速傅里叶变换及 IFFT快速傅里叶逆变换(迭代法和递归法) 一.二维FFTFFTFFT快速傅里叶变换公式推导二.二维FFTFFTFFT 及 IFFTIF ...
FFT快速傅里叶变换超详细的入门学习总结
FFT快速傅里叶变换说明本文创作的目的是为自己巩固该算法,加深印象并深入理解,同时也为FFT入门学者提供一份可鉴的学习总结. 原文链接:https://blog.csdn.net/qq_39565 ...
FFT快速傅里叶变换C语言实现信号处理对振动信号进行实现时域到频域的转换
FFT快速傅里叶变换C语言实现信号处理对振动信号进行实现时域到频域的转换,可实现FFT8192个点或改成其他FFT1024.4096等等,可以直接运行,运行结果与matlab运行的一致,写好了注释, ...
c语言fft乘法步骤,C语言实现FFT(快速傅里叶变换).doc
C语言实现FFT(快速傅里叶变换) 择蚁牙幸帆揣邓淌港烬粹甩滋整维含兔忿茂慨渔下餐随扼哇房坏鹅穆礼围引介害芝共茨恿把喜恤寇杖除冕嗓停揍猫调锚遭傀个碱晓频斌硕宾撕坪莱哩腊养掘蹄轴国繁蔬虞靡砖焙倍勾呸怀怒 ...
快速傅里叶变换c语言函数,C语言实现FFT(快速傅里叶变换)
while(1); } #include #include /********************************************************************* ...
如何 FFT(快速傅里叶变换) 求幅度、频率（超详细含推导过程）
目录如何 FFT(快速傅里叶变换) 求幅度.频率(超详细含推导过程) 一. 打颗栗子二. 求幅度 1. 快速傅里叶变换 2. 求出复数的绝对值 3. 归一化小结三. 求频率 1. 频率公式 ...
FFT快速傅里叶变换的应用——画单边频谱图matlab
FFT快速傅里叶变换的应用--画单边频谱图matlab 快速傅里叶变换在数字信号处理里用的十分广泛,在matlab仿真中,处理信号的时频域变换十分有效,这里结合两个做过的仿真,来说一说fft的应用:画 ...
FFT 快速傅里叶变换初探
一直认为很高深的东西其实也并不很难. 以下内容部分来自qy大神的ppt,同时结合了自己的理解.但理解还不是很深,需要继续研究. 开头首先什么是傅里叶变换:傅立叶变换能将满足一定条件的某个函数表示成三 ...
FFT快速傅里叶变换详解
介绍简而言之,这个东西用来做多项式乘法. 比如说,有两个多项式: 3x2+2x+1,2x2+x+43x^2+2x+1~,~2x^2+x+4 3x2+2x+1 , 2x2+x+4 那么他们乘起来就是 ...
MATLAB——FFT(快速傅里叶变换)
基础知识 FFT即快速傅里叶变换,利用周期性和可约性,减少了DFT的运算量.常见的有按时间抽取的基2算法(DIT-FFT)按频率抽取的基2算法(DIF-FFT). 1.利用自带函数fft进行快速傅里叶 ...

FFT快速傅里叶变换

FFT