FFT算法实现

引言

网上关于FFT的文章实在不少，但是理论居多，讲算法的也多，但是可供调试的程序实在不多，笔者写了一个，给大家做一下参考。

描述

1、蝶形算法：

有关蝶形算法的介绍和思想大家百度或者Google一下就很容易找到，这里只是说一下要注意的地方。

蝶形算法中有这样一个有趣的规则：若输入信号的顺序为自然顺序，那么输出信号的顺序就为倒位序（算法参见4）顺序。

2、二维FFT的变换顺序：

首先进行行变换，对变换后的结果再进行列变换。

3、关于二维FFT运算后的结果

按照公式运算出的结果中，能量大部分分集中在四个角，如果我们想要能量集中在中间，我们需要成一个欧拉数，其实也简单，你可以在输入信号时做一个简单的变换，如下描述：

设i，j为输入信号的坐标，那么

输入信号可表示为x(i, j)， 若(i + j) % 2 == 0 则取源信号为输入信号，否则取源信号的相反数为输入信号，即 -x(i, j)。

运算出来的结果中，能量就集中在中间位置了。

4、关于倒位序算法

倒位序：就是将数字的各个尾反过来排序后得到的数字后的顺序，举个例子吧

如我们的输入8个信号，我们只需要三个位就可以描述着写信号的下标，比如1 = 001B， 2 = 010B等等，那么1的倒位后为100B = 4， 010B = 2，依此类推，这就是倒位序，最后生成的新的顺序就是排序后的结果，这个结果有一个特点，那就是把偶数和奇数分开，这也就是FFT的理论基础。

注：但是这个排序需要收到输入信号总量的限制，比如我们输入8个信号，虽然在计算机里我们的取值范围为00000000B ~ 00000111B 但是倒位序的却只针对最后三位，即只在最后三位的范围内倒位排序，如上面1的倒位序是100B而不是10000000B，这点请非常注意！

下面就是倒位序的算法的应用，用C写的，可惜CSDN的粘贴代码不支持C，大家将就看吧，至于原理，谈不上，只是个技巧，大家用手简单模拟算一下就可以完全理解了~

这个片段讲述了如何进行到位需排序，假设输入信号的实部为x，虚部为y，这里j是标示着要与当前位置i交换的元素的下标，如果不交换的话就会出现i >= j 的情况，整个过程还是建议大家拿手算一下，这样来的比较快。

n为输入信号的个数，k只是个中间变量而已。

//  ...
j  =   0 ;
for  ( int  i  =   0 ; i  <  n  -   1 ; i ++ )  ... {
    if (i < j) ...{
        tx = x[i];
        ty = y[i];
        x[i] = x[j];
        y[i] = y[j];
        x[j] = tx;
        y[j] = ty;
    }
    k = n >> 1;
    while (k <= j) ...{
        j -= k;
        k >>= 1;
    }
    j += k;
}
//  ...

 

FFT变换步骤

做好了上面的准备，我们来概览一下FFT的变换步骤：

1、若输入信号的个数不是2的整数次方我们就要补齐信号，用零补齐就可以，详细参见笔者博客中《FFT实践中的一些补充》一文。

2、将输入信号进行倒位序排序

3、将输入信号成一个欧拉数，为的是输出结果中能量集中在中间，不乘也可以，视具体应用而定。

补充：只有一个变换方向可以产生能量及中的频谱图

4、进行傅立叶变换，若是二维FFT，那么先进行行变换，再对变换后的结果进行列变换，体现了二重积分的思想。

5、得到的结果通常不能直接输出谱图，因为值可能会很小，我们可以乘以255后输出，至于谱图，其实其实际意义没有多大，输出方法也不一样，比如matlab中就是取实部并乘以255然后输出，并且没有将能量中心移到屏幕中心，而笔者喜欢输出其乘以255后取模再除以100后的图。

例图：

 

第一幅是原图512 * 512 第二幅是它的谱图

注：这里使用的软件是笔者这个学期自行设计开发的，因为这个学期开设了数字图像处理，我非常感兴趣，并且还跟着系主任的研究生的课，呵呵，学到了不少，如果大家感兴趣，我可以把软件拿出来供大家下载，C + plain Win32 api 的程序，希望大家喜欢~

程序代码

 

    typedef  struct  tagComplex  ... {
        double real;
        double imag;
    }