利用C语言实现FFT、IFFT运算

一、FFT算法理论

上述分别为FFT、IFFT公式。下面首先讨论FFT的算法实现。
本文采用输入逆序、输出顺序的FFT计算方法。实质上就是在时域对x(n)进行“奇偶分类”、在频域上对X(k)进行“前后分类”。 值得说明的是，这里的“奇”和“偶”是相对的概念，并不完全是通常我们所理解的“奇”和“偶”。下面将给出一个例子进行说明：

·图中，总共可以分为 “三级” ，每一级包含若干 “蝶形计算单元”。每一个蝶形单元的结构如下：
因此，进行FFT运算只需考虑三件事情。一、将输入序列倒序排列，也就是进行奇偶分类。二、实现最基本的蝶形计算单元。三、根据规律，划分出FFT运算的级数，并且确定出每一级中参与蝶形运算的数据分组。

二、FFT运算的代码实现

1. 输入倒序排列（奇偶分类）。事实上，将输入倒序排列就是进行了奇偶分类。如果FFT的点数为8，那么这8个数的二进制表示分别为000，001，010，011，100，101，110，111。进行二进制层面倒序排列就变为000，100，010，110，001，101，011，111；重新转变为十进制就是0，4，2，6，1，5，3，7。 上面的计算方式对人类来说是很方便的，但是对于电脑来说，效果并不是太好，因此需要另外发现一种规律。
2. 输入奇偶分类（倒序排列）。上面已经提到，倒序排列实际上就是在进行奇偶分类，而进行奇偶分类对电脑来说会简单很多。总结1中提到的规律，三位二进制数可以理解为FFT运算的级数。我们可以发现，在倒序排列完成后，“偶数”总是排在“奇数”前面。 这里需要提出“偶中的偶”，“奇中的奇”的概念。即每次进行“余2”运算，结果是0的都排在前面。现在看来，电脑实现算法的方式就是每次进行“余2”运算，进行奇偶分类；然后进行“除2”运算，再进行“余2”运算，进行奇偶分类……直到最后排序完成。 实现代码如下：

void reverse()		//输入倒序
{ 
	complex temp;
	char* flag = (char*)malloc(N * sizeof(char));	//用于标记是否被换过，若为1，则已经换过
	int i, j, current_n, target_n;
	for (i = 0;i < N;i++) {
		flag[i] = 0;
	}
	for (i = 0;i < N;i++) {
		current_n = i;
		target_n = 0;
		//获取两个互为倒序的标号，换种思路
		for (j = 0;j < bits;j++) {
			target_n = target_n + (int)((current_n % 2) * pow(2, bits - j - 1));
			current_n /= 2;
		}
		current_n = i;
		//对应标号值互换
		if (current_n != target_n && flag[current_n] != 1) {
			temp = x[current_n];
			x[current_n] = x[target_n];
			x[target_n] = temp;
			flag[target_n] = 1;
		}
	}
	free(flag);
}

说明：complex为自己定义的复数数据结构。再次强调，奇偶分类就是倒序排序。二者是一样的。

3. 蝶形运算单元的实现。分析蝶形运算单元，两输入、两输出，因此蝶形运算单元需要注意的就是输入、输出采用同一片内存存储即可，也就是让输出数据覆盖输入数据，这样可以有效节约存储资源。

void butterfly(int x1_point, int x2_point, complex wn)	//蝶形计算单元
{
	complex result1, result2, T;
	T = multiplication(x[x2_point], wn);
	result1 = add(x[x1_point], T);
	result2 = sub(x[x1_point], T);

	x[x1_point] = result1;
	x[x2_point] = result2;
}

说明：上面multiplication、add、sub方法均为自己定义，wn为旋转因子。

4. 进行单级FFT运算。需要注意的就是总结规律。总结出两对偶结点的“结点跨距”、“分组间隔”、“分组数”……有了这些才能够编写将参与蝶形运算的数据分类配对。

void single_fft(int m)			//进行单级的fft运算
{
	//首先进行分组
	int point_distance = (int)pow(2, m - 1);	//结点跨距
	int group_distance = 2 * point_distance;	//分组间隔
	int group_count = N / group_distance;		//分组数

	int group_header = 0;

	int x1, x2;								//参与计算两结点的标号
	complex w;

	int i, j;								//循环控制符
	//分组
	for (i = 0;i < group_count;i++) {
		//获取一组的标号范围
		group_header = i * group_distance;
		//将每组分成两半，前一半均为x1,后一半均为x2
		for (j = 0;j < group_distance / 2;j++) {
			w = w_builder(m, j);
			x1 = j + group_header;
			x2 = x1 + point_distance;
			butterfly(x1, x2, w);
		}
	}
}

5. 有了单级FFT运算以后，将每一级组合就可以实现FFT运算了。

complex* fft(int fft_point, int xn_length, complex* xn)		//fft最终封装
{
	int i;
	//初始化
	fft_init(fft_point, xn_length, xn);
	//输入倒序
	reverse();
	//fft运算
	for (i = 1;i <= bits;i++) {
		single_fft(i);
	}
	return x;
}

三、IFFT算法理论

从上面的等式可以看出，括号里面的内容实际上就是FFT运算。也就是说只需要对输入、输出序列进行一定的处理，就可以利用刚刚已经做好的FFT算法实现IFFT的运算！

四、IFFT运算的代码实现

complex* ifft(int xk_length, complex* xk)		//ifft最终封装
{
	int i;
	//对输入序列进行处理
	ifft_init(xk_length, xk);
	//套用fft算法
	reverse();
	for (i = 1;i <= bits;i++) {
		single_fft(i);
	}
	//对输出序列进行处理
	for (i = 0;i < N;i++) {
		x[i] = conjugate(x[i]);		//求共轭
		x[i].real /= N;
		x[i].imaginary /= N;
	}
	return x;
}

五、总结

本文也程序也考虑到了FFT自动补零问题的处理，大家可以参考完整代码（FFT、IFFT均验证正确）。FFT.zip