快速傅里叶变换FFT

最新推荐文章于 2025-06-30 16:28:41 发布

weixin_30679823

最新推荐文章于 2025-06-30 16:28:41 发布

阅读量76

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/maoruimas/p/9726524.html

本文深入探讨了快速傅立叶变换(FFT)算法，详细解释了其基于分治思想对离散傅立叶变换(DFT)和逆离散傅立叶变换(IDFT)的优化原理。FFT将复杂数运算的复杂度从O(n^2)降至O(nlogn)，适用于多项式快速乘法。文章还提供了FFT的代码实现，包括递归和模拟递归两种方法，并讨论了在多项式乘法中的应用。

FFT是基于分治思想对DFT和IDFT的优化。用于多项式的快速乘法

DFT: \(\mathbb{C} ^{N}\rightarrow \mathbb{C} ^{N}, \left( x_{1},\ldots ,x_{N}\right) \mapsto \left( y_{1},\ldots ,y_{N}\right)\)
IDFT是DFT的逆映射
记\(\omega _{N}=e^{\dfrac {2\pi i}{N}}\)，则
\[DFT: y_{k}=\sum ^{N-1}_{n=0}\omega ^{nk}_{N}x_{n}\\ IDFT: x_{k}=\sum ^{N-1}_{n=0}\dfrac {\omega ^{-nk}_{N}}{N}y_{n}\]
于是，DFT与IDFT都是N维复空间上的线性变换，其矩阵为Vandermonde矩阵。由于是矩阵乘法，复杂度为O(n^2)。注意到复单位根的性质使得DFT与IDFT高度相似，因此，只需要考虑降低DFT的复杂度即可。

定义：
\[E_{k}=\sum ^{N/2-1}_{m=0}\omega ^{2mk}_{N}x_{2m}=\sum ^{N/2-1}_{m=0}\omega ^{mk}_{N/2}x_{2m}\\ O_{k}=\sum ^{N/2-1}_{m=0}\omega ^{2mk}_{N}x_{2m+1}=\sum ^{N/2-1}_{m=0}\omega ^{mk}_{N/2}x_{2m+1}\]
则\[y_{k}=E_{k}+\omega ^{k}_{N}O_{k}\\ y_{k+N/2}=E_{k}-\omega ^{k}_{N}0_{k}\]
显然，如此递归可以将复杂度降至O(nlogn)。另外，N必须是2的幂次。

应用
由多项式的知识知\[\left\{ f\left( x\right) |\deg f=n-1\right\} \cong \left\{ \left( a_{0},\ldots ,a_{n-1}\right) |a_{n-1}\neq 0\right\} \\ \cong \left\{ \left\{ \left( x_{i},f\left( x_{i}\right) \right) \right\} |x_{i}\neq x_{j},\forall i,j\right\}\]
于是，通过FFT可以快速将多项式函数的“系数表示”转化为“点值表示”。
当计算多项式的乘积时，我们知道用“点值表示”计算是方便的。因此可以先将待乘多项式转化为“点值表示”，再将结果转回“系数表示”。

代码（递归）

void separate (complex<double>* a, int n) {
    complex<double>* b = new complex<double>[n/2];
    for(int i=0; i<n/2; i++) 
        b[i] = a[i*2+1];
    for(int i=0; i<n/2; i++) 
        a[i] = a[i*2];
    for(int i=0; i<n/2; i++) 
        a[i+n/2] = b[i];
    delete[] b;
}
// N must be a power-of-2
//DFT: sig=1
//IDFT: sig=-1
// N input samples in X[] are FFT'd and results left in X[]
void fft2 (complex<double>* X, int N,int sig) {
    if(N < 2) {
        // bottom of recursion.
    } else {
        separate(X,N);
        fft2(X,     N/2,sig);
        fft2(X+N/2, N/2,sig);
        // combine results of two half recursions
        for(int k=0; k<N/2; k++) {
            complex<double> e = X[k    ];   // even
            complex<double> o = X[k+N/2];   // odd
                         // w is the "twiddle-factor"
            complex<double> w = exp( complex<double>(0,sig*2.*M_PI*k/N) );
            X[k    ] = e + w * o;
            X[k+N/2] = e - w * o;
        }
    }
}

模拟递归。
参考网上的博客，这种模拟首先需要做一种置换：对指标二进制表示下对称的元素做对换。（至于如何实现这一置换，可以反向做二进制加法。例如1100的后一位是0010，正如0011的后一位是0100一样。过程见代码。）接着，如同递归函数从栈顶到栈底的过程一样，做模拟的递归过程。

void rearrange(complex<double> x[],int n)
{
    for(int i=1,j=n/2;i<n;++i)
    {
        if(i<j)swap(x[i],x[j]);
        int tmp=n/2;
        while(tmp&&j>=tmp){j-=tmp;tmp/=2;}
        if(j<tmp)j+=tmp;
    }
}
void fft(complex<double> x[],int n,int sig)
{
    rearrange(x,n);
    for(int i=2;i<=n;i*=2)
    {
        for(int j=0;j<n;j+=i)
        {
            for(int k=j;k<j+i/2;++k)
            {
                complex<double> e=x[k],o=x[k+i/2],w=exp(complex<double>(0,sig*2.*pi*(k-j)/i));
                x[k]=e+w*o;
                x[k+i/2]=e-w*o;
            }
        }
    }
    if(sig==-1)
    {
        for(int i=0;i<n;++i)x[i]/=n;
    }
}

转载于:https://www.cnblogs.com/maoruimas/p/9726524.html