快速傅里叶变换fft_FFT快速傅里叶变换的推导

离散傅里叶的计算

从我的学习经验来看，如果不是对DFT非常熟悉的，可能会在理解FFT算法的时候有点懵逼。

首先搞清楚，DFT计算的是什么及其公式。

其中

展开后得：

…………

上面多项式可转化为下面的矩阵：

​

​​最后看到离散傅里叶就是一个复数的矩阵乘法计算。

复指数的周期说明

从上面看到矩阵里的元素全是复指数的幂，根据复指数W的周期性，可以推导出矩阵中相同值的元素。

再有：

DFT计算的拆分

把DFT计算公式按偶数和奇数拆分成两个。

​

​​此处省略了几条公式，本人懒得写了。

总结下来得到下面公式。

G(k)是偶数部分，H(k)是基数部分，并且是以周期为N/2的周期函数。

由(1)(2)(3)(4)得：

以N=4为例，

​

​​从上面可以看出，x(0)和x(2)的“输入参数”是一模一样的。

个人觉得蝶形运算单元不需要说的特别明白，只需要知道她到底是做了什么运算就可以了。

这里拿X(0)和X(2)做说明。

​

​​看出这里的优化就是多用了一个临时寄存器C来减少一次乘法运算。

下面是FFT比较关键的一部分内容。

​

​​文字写的太草，大概意思是，X(k)，G(k)，H(k)的运算时一样的，只是参数不一样，那么我以函数形式来说明一下：

XGK(x数组[长度是N], W[周期是N])
{
    //对数组x加权（权是W）求和
}

所以快速傅里叶变换是一个递归运算。

FFT算法实现

下面是重点中的重点了。书本有个叫位码倒读的优化概念，不知道是不是本人的理解力有问题，书本里说的内容我觉得从编程的角度去理解不太好弄。

我尝试用N=8的数组来解释。

从图中可以看到，输入的小写x[n]已经通过倒位码排了序，输出的大写X[k]是按顺序输出。

​

​​当N=8时，log2(8)=3，共3级，每级记做M。从图里看第一级（M=0）中有4组碟形，每组只有一个蝶形运算。

把组记做J，一个蝶形记做k，一个蝶形运算的数组间隔是L。

那么对应的规律是：

L = 1<<M；

J = N / (2 * L)；

K=L；

还有一个很重要的规律，每一组的下标间隔是l*2，即j += L*2；

更重要的来了，每一个蝶形运算的下标怎么取？

从图中可以看出：蝶形的上级是J + K，下级是上级+L。

下面上一下伪代码。

//假设X = X1 * W*X2
for (M = 0; M < log2(N); ++M)  //第M级
{
    int l = 1<<M;   
    for (j = 0; j < N; j += 2*l)   //第j组
    {
        for (k = 0; k < l; ++k)   //第k个蝶形运算
        {
            //X1起始下标：
            int index1 = j + k;
            //X2起始下标：
            int index2 = index1 + l;
            
            rightUp = x[index1] + W*x[index2];
            rightDown =  x[index1] - W*x[index2];
        }
    }
}

再看上图，输出的结果的下标和输入的下标完全对应，即书本里说的“即位运算”！

完整快速傅里叶变换代码（参考：https://blog.csdn.net/tuwenqi2013/article/details/71772841）

void fft(complex x[], int N, complex* W)
{
    complex up, down, product;
    change(x, N);  //调用变址函数  
    int size = log(N) / log(2);
    for (int M = 0; M < size; M++)           //第M级
    {
        int l = 1 << M;
        for (int j = 0; j < N; j += 2 * l)    //第j组
        {
            for (int k = 0; k < l; k++)         //第k个蝶形运算
            {
                int wIndex = N * k / 2 / l;
                int index1 = j + k;
                int index2 = index1 + l;
                product = x[index2] * W[wIndex];
                up = x[index1] + product;
                down = x[index2] - product;
                x[index1] = up;      //即位运算
                x[index2] = down;
            }
        }
    }
}