引言
网上关于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;
} COMPLEX, *LPCOMPLEX;
void Powerof2(int num, int * m) ...{
int n = 0, tn = 1;
while (tn * 2 <= num) ...{
n++;
tn *= 2;
}
*m = n;
}
void FFT(int dir, int m, double * x, double * y) ...{
int n = 1, k, j, l1, l2, i2;
double tx, ty, c1 = 1.0, c2 = 0.0, u1, u2, t;
for (int i = 0; i < m; i++) n *= 2;
// 换位
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;
}
// 变换开始
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (int i = 0; i < m; i++) ...{
u1 = 1.0;
u2 = 0.0;
l1 = l2;
l2 <<= 1;
for (int j = 0; j < l1; j++) ...{
for (int p = j; p < n; p += l2) ...{
i2 = p + l1;
tx = u1 * x[i2] - u2 * y[i2];
ty = u1 * y[i2] + u2 * x[i2];
x[i2] = x[p] - tx;
y[i2] = y[p] - ty;
x[p] += tx;
y[p] += ty;
}
t = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = t;
}
c2 = sqrt((1.0 - c1) / 2);
c1 = sqrt((1.0 + c1) / 2);
if (dir == 1) c2 = -c2;
}
// 逆变换系数
if (dir == 1) ...{
for (int i = 0; i < n; i++) ...{
x[i] /= (double)n;
y[i] /= (double)n;
}
}
}
void FFT2D(COMPLEX ** c, int nrow, int ncol, int dir) ...{
int m;
double *real, *imag;
// 行变换
Powerof2(ncol, &m);
real = (double*)malloc(sizeof(double) * ncol);
imag = (double*)malloc(sizeof(double) * ncol);
for (int i = 0; i < nrow; i++) ...{
for (int j = 0; j < ncol; j++) ...{
real[j] = c[i][j].real;
imag[j] = c[i][j].imag;
}
FFT(dir, m, real, imag);
for (int j = 0; j < ncol; j++) ...{
c[i][j].real = real[j];
c[i][j].imag = imag[j];
}
}
free(real);
free(imag);
// 列变换
Powerof2(nrow, &m);
real = (double*)malloc(sizeof(double) * nrow);
imag = (double*)malloc(sizeof(double) * nrow);
for (int i = 0; i < ncol; i++) ...{
for (int j = 0; j < nrow; j++) ...{
real[j] = c[j][i].real;
imag[j] = c[j][i].imag;
}
FFT(dir, m, real, imag);
for (int j = 0; j < nrow; j++) ...{
c[j][i].real = real[j];
c[j][i].imag = imag[j];
}
}
free(real);
free(imag);
}
这是笔者软件中数学包里的一部分程序，希望大家喜欢啊，记得留言~(*^__^*) 嘻嘻……
发表于 @ 2008年05月17日 13:24:10|评论(loading...)|编辑