引言
网上关于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 的程序,希望大家喜欢~
程序代码
double real;
double imag;
}