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

给定两个多项式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^s$$w_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电路中会更加清楚：

// 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)$