FFT

时间:2020-01-06
本文章向大家介绍FFT,主要包括FFT使用实例、应用技巧、基本知识点总结和需要注意事项,具有一定的参考价值,需要的朋友可以参考一下。

\[ FFT(快速傅里叶变换)\]

多项式

定义:多项式是指由变量、系数以及它们之间的加、减、乘、幂运算(非负整数次方)得到的表达式。

  • :多项式中的每个单项式叫做多项式的项。

  • 多项式的次数:这些单项式中的最高项次数,就是这个多项式的次数。

一个\(n\)次多项式可以表示为\(\sum\limits_{i=0}^na_ix^i\),其中\(a_i\)为系数。

复数

定义:我们把形如\(z=a+bi\)\(a\),\(b\)均为实数)的数称为复数,其中\(a\)称为实部,\(b\)称为虚部,\(i\)称为虚数单位,(\(i^2=-1\)

\(z_1=a+bi\),\(z_2=c+di\)

\[z_1\pm z_2=(a\pm c)+(b\pm d)i\]

\[z_1\times z_2=(ac-bd)+(ad+bc)i\]

\[z_1 \div z_2=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i\]

单位根

定义:记方程\(\omega^n=1(n\in N^{+})\)\(n\)个复数解为\(n\)次单位根。

可以解出\(\omega=e^{\frac{2k\pi i}{n}}\)

定义主次单位根\(w_n=e^{\frac{2\pi i}{n}}\)

对于所有单位根,记作:\(w_n^k=e^{\frac{2k\pi i}{n}}\)

多项式乘法卷积

对于\(n\)次多项式\(A,B\)\(A=\sum\limits_{i=0}^na_ix^i\),\(B=\sum\limits_{i=0}^nb_ix^i\)

\(A\)\(B\)的乘法卷积为\(C\)\(C=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{n}a_ib_jx^{i+j}\)

多项式的点值表示

众所周知,一个\(n\)次函数可以由\(n+1\)个点对\((x_i,y_i)\)唯一确定(其中\(x_i\)互不相同)

因此可以用\(n+1\)个点对表示出一个\(n\)次的多项式。

\(A=((x_1,y_1),(x_2,y_2),\dots,(x_n,y_n))\)\(B=((x_1,y'_1),(x_2,y_2'),\dots,(x_n,y_n'))\)

\(A*B=((x_1,y_1y_1'),(x_2,y_2y_2'),\dots,(x_n,y_ny_n'))\)

所以如果多项式的点值已知,那么多项式的乘法是\(O(n)\)的。

若把点对的数量扩大至\(2n+1\),则可以唯一确定\(C\)

DFT & IDFT

DFT(离散傅里叶变换) 是系数数列 \(a_0,a_1,\cdots,a_n\) 转为点值数列的过程,IDFT(逆离散傅里叶变换)DFT 的逆运算。

注意到单位根存在以下性质

\[\omega_{2n}^{k}+\omega_{2n}^{n+k}=0\]

\[(\omega_{2n}^k)^2=\omega_{n}^k\]

考虑用单位根的性质优化 DFT。

对于多项式 \(A(x)=\sum_{i=0}^{2n-1}a_ix^i\),其中 \(2n=2^k(k\in N^{+})\),我们将 \(\omega_{2n}^0,\omega_{2n}^1,\cdots,\omega_{2n}^{{2n}-1}\) 代入公式计算点值。

现在重定义两个多项式

\[A_0(x)=a_0+a_2x+a_4x^2+\cdots+a_{2n}x^n\]

\[A_1(x)=a_1+a_3x+a_5x^2+\cdots+a_{2n-1}x^n\]

显然

\[A(x)=A_0(x^2)+xA_1(x^2)\]

将单位复根代入上式

\[\begin{aligned}A(\omega_{2n}^k) & =A_0((\omega_{2n}^k)^2)+\omega_{2n}^kA_1((\omega_{2n}^k)^2)\\&=A_0(\omega_n^k)+\omega_{2n}^kA_1(\omega_n^k)\\A(\omega_{2n}^{n+k})&=A_0((\omega_{2n}^{n+k})^2)+\omega_{2n}^{n+k}A_1((\omega_{2n}^{n+k})^2)\\&=A_0(\omega_n^k)-\omega_{2n}^kA_1(\omega_n^k)\end{aligned}\]

发现对于 \(k\in[0,1,\cdots,n-1]\) \(A(\omega_{2n}^k)\)\(A(\omega_{2n}^{n+k})\) 是可以在一起计算的。

于是有以下伪代码

function FFT(complex A[], int Length)
    if Length == 1 return
    complex A0[Length / 2], A1[Length / 2]
    for int i = 0 to Length / 2 - 1 with i += 1
        A0[i] = A[i * 2]
        A1[i] = A[i * 2 + 1]
    FFT(A0, Length / 2)
    FFT(A1, Length / 2)
    complex wn = (cos(2 * Pi / Length), sin(2 * Pi / Length))
    complex w = (1, 0)
    for int i = 0 to Length / 2 - 1 with i += 1, w *= wn
        A[i] = A0[i] + A1[i] * w
        A[i + Length / 2] = A0[i] - A1[i] * w

考虑 IDFT,IDFT 是 DFT 的逆运算,令 \(DFT(\{a_i\})=\{b_i\}\),已知

\[b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\]

存在结论

\[a_k=\frac 1 n\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\]

证明,将前式代入后式

\[\begin{aligned}a_k&=\frac 1 n\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\\&=\frac 1 n\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}\\&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_n^{-ki}\omega_{n}^{ij}\\&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{i(j-k)}\end{aligned}\]

考虑 \(\sum_{i=0}^{n-1}\omega_n^{i(j-k)}\)

\(j=k\)
\[\sum_{i=0}^{n-1}\omega_n^{i(j-k)}=\sum1=n\]

\(j\neq k\)
\[\begin{aligned}\sum_{i=0}^{n-1}\omega_n^{i(j-k)}&=\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\&=\frac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}\\&=\frac{1-(\omega_n^n)^{j-k}}{1-\omega_n^{j-k}}=\frac{1-1}{1-\omega_n^{j-k}}=0\end{aligned}\]

所以,

\[\begin{aligned}a_k&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{i(j-k)}\\&=\frac 1 n\cdot na_k=a_k\end{aligned}\]

发现\(DFT\)\(IDFT\)的公式形式一样,所以可以用一个函数实现。

蝴蝶定理

考虑 FFT 的分治过程,以 \(n=16\) 为例

\[\{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15\}\]

\[\{0,2,4,6,8,10,12,14\},\{1,3,5,7,9,11,13,15\}\]

\[\{0,4,8,12\},\{2,6,10,14\},\{1,5,9,13\},\{3,7,11,15\}\]

\[\{0,8\},\{4,12\},\{2,10\},\{6,14\},\{1,9\},\{5,13\},\{3,11\},\{7,15\}\]

其二进制下表示,

\[\{0000,1000\},\{0100,1100\},\{0010,1010\},\{0110,1110\},\]

\[\{0001,1001\},\{0101,1101\},\{0011,1011\},\{0111,1111\}\]

观察发现,若干次蝴蝶操作(奇偶分治)后,其数位二进制翻转后是连续的。

对于二进制翻转,可用递推计算,rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0)

考虑用蝴蝶定理使 \(FFT\) 的过程避免递归,可以先将 \(\{a_i\}\)\(rev\) 序重新排列,在分治树上从下往上,

function FFT(complex A[], int Length, int on)
    for int i = 0 to Length - 1 with i += 1
        if i < Rev[i]
            swap(A[i], A[Rev[i]])
    for int k = 1 to n - 1 with k *= 2
        complex wn = (cos(Pi / k), sin(on * Pi / k))
        for int i = 0 to n with i += k * 2
            complex w = (1, 0)
            for int j = i to i + k - 1 with j += 1
                complex x = a[j]
                complex y = a[j + k]
                a[j] = x + y
                a[j + k] = x - y
    if on == -1
        for int i = 0 to n - 1 with i += 1
            a[i] /= n

然后枚举当前合并的长度 \(2k\),枚举合并区间开始位置 \(i\),枚举区间中的元素 \(a_j\)


struct node{
    double x,y;
}A[MAXN],B[MAXN];
node operator+ (node p,node q){return (node){p.x+q.x,p.y+q.y};}
node operator- (node p,node q){return (node){p.x-q.x,p.y-q.y};}
node operator* (node p,node q){return (node){p.x*q.x-p.y*q.y,p.x*q.y+p.y*q.x};}
void FFT(node *X,int flag){
    for (register int i=0;i<len;i++) if (i<rev[i]) swap(X[i],X[rev[i]]);
    for (register int l=1;l<len;l<<=1){
        node T=(node){cos(Pi/l),flag*sin(Pi/l)};
        for (int s=0;s<len;s+=(l<<1)){
            node t=(node){1,0};
            for (int k=0;k<l;k++,t=t*T){
                node Nx=X[s+k],Ny=t*X[s+k+l];
                X[s+k]=Nx+Ny,X[s+k+l]=Nx-Ny;
            }
        }
    }
    for (register int i=0;i<len;i++) A[i].x/=len;
    return;
}
int main(){
    n=qr(),m=qr();
    for (int i=0;i<=n;i++) A[i].x=qr();
    for (int i=0;i<=m;i++) B[i].x=qr();
    while (len<=n+m) len<<=1,++L;
    for (int i=0;i<len;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    FFT(A,1);FFT(B,1);
    for (int i=0;i<=len;i++) A[i]=A[i]*B[i];
    FFT(A,-1);
    for (int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].x+0.5));
    return 0;
}

原文地址:https://www.cnblogs.com/bestlxm/p/12156371.html