美文网首页
多项式与快速傅里叶变换

多项式与快速傅里叶变换

作者: 漫游之光 | 来源:发表于2019-10-26 21:07 被阅读0次

    本来这里我想写仔细一点,但我想,傅里叶变换确实比较复杂,况且,我也无法在没有任何参照的情况下写出下面的代码,所以我也就贴上代码,然后写写我在实现过程中的一些想法。下面的代码基本上都是和《算法导论》上面的实现一样的。

    多项式的两种表示方法

    说实话,第一次看到多项式式的点值形式的时候,我是比较惊讶的,因为在我的第一感觉是,点值无法确定一个多项式,但书中用了一个很简单的说明就让我明白了我的直觉是错的。对于多项式点值表达式的乘法也让我印象深刻。其实比较简单,如果C(x) = A(x) * B(x),那么对于点x_0,必然存在C(x_0) = A(x_0) * B(x_0),也就是对应两个点的直接相乘。

    以前学《信号与系统》的时候,我就知道时域卷积对应于频域相乘,但是我还是不太清楚为什么,这里给了我一点直观的感觉。

    算法的基本思路

    多项式相乘,如果使用直接的方法,那么需要O(n^2)的时间,而点值形式的多项式乘法只需要O(n),这时候我们很自然的想到,能否找到一种快速的方法,实现两种形式的转化。

    先考虑从系数形式到点值形式的转化,根据霍纳法则,我们可以在O(n)的时间内求得一个点对应的值,但一共有n个点,也需要O(n^2)。这时候,我们想到了一些算法设计的技巧,比如分治法。我们可以把多项式按照奇偶划分成两部分,然后分别求解,然后合并。则看上去很不错,但写出递归式:
    T(n,n)=2*T(n/2,n)+O(n,n) = O(n^2)
    这里不好的地方在于,递归过程中,我们要求的n个点没有减少,如果我们能在递归过程中,减少对点的求解,也就是说,让根之间有关系,减少计算,如果能减少到一半,那么递归式变成了:
    T(n,n)=2*T(n/2,n/2)+O(n,n/2) = O(nlogn)
    这样,这个算法就比直接求解要优秀。

    这里,选择的特别的根是单位复数根,它有一个很重要的性质。
    \begin{array}{l}{f\left(\omega_{n}^{i+n / 2}\right)=f\left(-\omega_{n}^{i}\right)=} \\ {a_{0}+a_{2} \omega_{n / 2}^{i}+a_{4}\left(\omega_{n / 2}^{i}\right)^{2} \ldots+a_{n-2}\left(\omega_{n / 2}^{i}\right)^{(n-2) / 2}} \\ {-\omega_{n}^{i}\left(a_{1}+a_{3} \omega_{n / 2}^{i}+a_{5}\left(\omega_{n / 2}^{i}\right)^{2} \ldots+a_{n-1}\left(\omega_{n / 2}^{i}\right)^{(n-2) / 2}\right)}\end{array}

    利用这个性质,可以在计算的过程中,实现减少点的计算。

    计算完了,我们还要转回来,根据复数根的性质,我们发现,这个矩阵的逆矩阵很好求,就是虚部取反,除以n,然后问题有回到了前面。

    #include <iostream>
    #include <complex>
    #include <vector>
    #include <cmath>
    #include <algorithm>
    using namespace std;
    
    const double PI = acos(-1);
    
    void fft(vector<complex<double> > &a ,int n, int op){
        if(n == 1){
            return;
        }
        //divide
        vector<complex<double> > a0(n/2,0);
        vector<complex<double> > a1(n/2,0);
        for(int i=0;i<n/2;i++){
            a0[i] = a[i*2];
            a1[i] = a[i*2+1];
        }
        //conquer
        fft(a0,n/2,op);
        fft(a1,n/2,op);
        //merge
        complex<double> w(1,0); //通过巧妙的选取根来减少计算
        complex<double> wn(cos(2 * PI / n),sin(2 * PI * op /n));
        for(int i=0;i<n/2;i++){  //如果选取的根不特殊的话,需要对n个值分别进行计算
            a[i] = a0[i] + w*a1[i];
            a[i+n/2] = a0[i] - w * a1[i]; //根据单位负数根的性质可以很快求出另一个根对应的值 
            w = w * wn;
        }
    }
    
    void fft2(vector<complex<double> > &a ,int n, int op){
        //开始输入的时候位置为对应计数比特翻转之后的值
        int len = (int)log2(n) - 1;
        vector<int> rev(n,0);
        for(int i=0;i<n;i++){ //计算翻转位置
            //利用和位置为i/2的翻转数的关系,提高速度
            rev[i] = (rev[i>>1]>>1) | ((i&0x01)<<len);
        }
        //改变输入数组的顺序
        for(int i=0;i<n;i++){
            if(i < rev[i]){
                swap(a[i],a[rev[i]]);
            }
        }
    
        //模拟递归的计算顺序,从底层开始,一个log2(n)层
        for(int i = 1;i< n; i *= 2){
            complex<double> wn(cos( PI / i),sin(PI * op / i));
            for(int k=0;k<n;k+= 2 * i){
                complex<double> w(1,0);
                for(int j=0;j<i;j++){
                    complex<double> t = w * a[k + j + i];
                    complex<double> u = a[k + j];
                    a[k+j] = u + t;
                    a[k+j+i] = u - t;
                    w *= wn;
                }
            }
        }
        if(op == -1){
            for(int i=0;i<n;i++){
                a[i] /= n;
            }
        }
    }
    
    /*
    2 2
    1 2 1
    2 0 1
    2 4 3 2 1
    */
    
    int main(int argc, char const *argv[])
    {
        int n,m;
        cin>>n>>m;
        int n1= 1;
        while(n1<=m+n){
            n1 <<= 1;
        }
        vector<complex<double> > v1(n1,0);
        vector<complex<double>> v2(n1,0);
        vector<complex<double>> res(n1,0);
        int tmp;
        for(int i=0;i<=n;i++){
            cin>>v1[i];
        }
        for(int i=0;i<=m;i++){
            cin>>v2[i];
        }
        //将系数形式转化为点值形式
        fft(v1,n1,1);
        fft(v2,n1,1);
        
        //利用点值形式计算多项式乘法
        for(int i=0;i<n1;i++){
            res[i] = v1[i] * v2[i];
            // cout<<v1[i]<<" "<<v2[i]<<" "<<res[i]<<endl;
        }
        
        //将点值表达转化为系数形式
        fft(res,n1,-1);
        for (int i = 0; i <= m+n; i++){
            printf("%.0lf ", floor(res[i].real()/n1));
        } 
        system("pause");
        return 0;
    }
    

    代码的注释写得挺多的,我个人只能说对这个算法有一个基本的理解,主要是为了记录学习的过程。如果有错误,请多多包涵。

    相关文章

      网友评论

          本文标题:多项式与快速傅里叶变换

          本文链接:https://www.haomeiwen.com/subject/rvxnpctx.html