📜  多项式乘法的迭代快速傅立叶变换

📅  最后修改于: 2021-04-24 21:53:32             🧑  作者: Mango

给定两个多项式A(x)和B(x),找到乘积C(x)= A(x)* B(x)。在上一篇文章中,我们讨论了解决O(nlogn)复杂度问题的递归方法。
例子:

输入:A [] = {9,-10,7,6} B [] = {-5,4,0,-2}输出:C(x)=  -12x^6 -14x^5 +44x^4 -20x^3 - 75x^2 +86x -45

在现实生活中的应用中,例如信号处理,速度非常重要,本文探讨了一种高效的FFT实现。本文重点介绍在O(nlogn)时间运行的FFT算法的迭代版本,但与递归版本相比,其隐藏常量较低,此外,它还节省了递归堆栈空间。
先决条件:FFT的递归算法。
在以下内容的for循环评估中,回想起前一篇文章的递归FFT伪代码:  y_k  w_n^k 计算两次。我们可以将循环更改为仅计算一次,将其存储在临时变量t中。因此,它变成了
\leftarrow 0至n / 2 – 1
不要 \leftarrow \omega y^{\left [ 1 \right ]}
 y_{k}\leftarrow \omega y_{k}^{\left [ 0 \right ]} + t
 y_{k + \left ( n/2 \right )}\leftarrow _{k}^{\left [ 0 \right ]} - t
 \omega \leftarrow \omega \omega _{n}

该循环中的运算乘以旋转因子w =  w_n^k 经过 y_k^[^1^] ,将乘积存储到t中,然后从中减去t  y_k^[^0^] ,称为蝶式运算。
如图所示,这是蝶形运算的样子:
蝴蝶手术
让我们考虑n = 8并继续进行迭代fft算法的形成。查看上面的递归树,我们发现,如果将初始系数矢量a的元素排列为它们在叶子中出现的顺序,则可以跟踪Recusive-FFT过程的执行,但是从下至上而不是从上到下下。首先,我们将元素成对,使用一次蝶形运算计算每对元素的DFT,然后用其DFT替换该元素对。然后,向量保存n / 2个2元素DFT。接下来,我们将这些n / 2个DFT成对使用,并通过执行两次蝶形运算(用两个4个元素的DFT替换两个2个元素的DFT)来计算它们来自四个矢量元素的DFT。然后,该向量包含n / 4个4元素DFT。我们以这种方式继续进行,直到向量拥有两个(n / 2)个元素DFT,我们使用n / 2个蝶形运算将其合并为最终的n个元素DFT。
递归树
要将这种自下而上的方法转换为代码,我们使用数组A [0…n],该数组最初按输入向量a的元素在树的叶子中出现的顺序保存这些元素。因为我们必须在树的每个级别上结合DFT,所以我们引入变量s来计算级别,范围从1(在底部,当我们将对组合成2元素DFT时)到lgn(在顶部,当将两个n / 2个元素DFT组合在一起以产生最终结果时)。因此,算法为:

1.对于s = 1到lgn 2.对k = 0到n-1进行 2^s 3.结合两者 2^s-1  A [k … k +]中的元素DFT  2^s-1 -1]和A [k + 2^s-1 … k + 2^s -1]放入A [k … k + 2^s -1]

现在,为了生成代码,我们按叶的顺序排列系数向量元素。示例-叶子0、4、2、6、1、5、3、7中的顺序是索引的一点反转。从000,001,010,011,100,101,110,111开始,反转每个二进制数的位以获得000,100,010,110,001,101,011,111,迭代FFT的伪代码:

BIT-REVERSE-COPY(a,A)n =长度[a],对于k = 0到n-1,做A [rev(k)] = a [k]迭代FFT BIT-REVERSE-COPY(a,A) n =长度(a),其中s = 1以记录n做m = 2^sw_m = e^(2*PI*i/m)对于j = 0至m / 2-1,对于k = j至n-1,由m做t = w A [k + m / 2] u = A [k] A [k] = u + t A [k + m / 2] = ut w = w*w_n返回A

从下面的并行FFT电路中会更加清楚:
并行FFT

// CPP program to implement iterative
// fast Fourier transform.
#include 
using namespace std;
  
typedef complex cd;
const double PI = 3.1415926536;
  
// Utility function for reversing the bits
// of given index x
unsigned int bitReverse(unsigned int x, int log2n)
{
    int n = 0;
    for (int i = 0; i < log2n; i++)
    {
        n <<= 1;
        n |= (x & 1);
        x >>= 1;
    }
    return n;
}
  
// Iterative FFT function to compute the DFT
// of given coefficient vector
void fft(vector& a, vector& A, int log2n)
{
    int n = 4;
  
    // bit reversal of the given array
    for (unsigned int i = 0; i < n; ++i) {
        int rev = bitReverse(i, log2n);
        A[i] = a[rev];
    }
  
    // j is iota
    const complex J(0, 1);
    for (int s = 1; s <= log2n; ++s) {
        int m = 1 << s; // 2 power s
        int m2 = m >> 1; // m2 = m/2 -1
        cd w(1, 0);
  
        // principle root of nth complex 
        // root of unity.
        cd wm = exp(J * (PI / m2)); 
        for (int j = 0; j < m2; ++j) {
            for (int k = j; k < n; k += m) {
  
                // t = twiddle factor
                cd t = w * A[k + m2]; 
                cd u = A[k];
  
                // similar calculating y[k]
                A[k] = u + t; 
  
                // similar calculating y[k+n/2]
                A[k + m2] = u - t; 
            }
            w *= wm;
        }
    }
}
  
int main()
{
    vector a{ 1, 2, 3, 4 };
    vector A(4);
    fft(a, A, 2);
    for (int i = 0; i < 4; ++i)
        cout << A[i] << "\n";
}
Input:  1 2 3 4
Output:
(10, 0)
(-2, -2)
(-2, 0)
(-2, 2)

时间复杂度分析:
复杂度为O(nlgn)。为了说明这一点,我们将在O(nlgn)时间中运行的最里面的循环显示为:
L(n) = \sum_{s=1}^{lg \ n}\sum_{j=0}^{2^{s-1}-1}\frac{n}{2^{s}}
=\sum_{s=1}^{lg \ n}\frac{n}{2^{s}}.2^{s-1}
 =\sum_{s=1}^{lg \ n}\frac{n}{2}
 =\Theta (n \ lg \ n)