快速傅里叶变换算法探幽

学完《信号与系统》之后发现什么还是不会？感觉《数字信号处理》和《信号与系统》差不多？怎么办？期末快到了……是不是感觉《数字信号处理》学了和没学一样？对的，就是这种感觉。大一的时候用FFT算法做过音乐频谱，然而当时对里面的算法并不是很了解，之前面某俱乐部的时候被面试官问到，因此……虽然还是……。学完DSP之后感觉没自己来尝试下FFT算法真是人生一大憾事。于是赶着期末之前当做是复习下（虽然并不在考试范围）。

傅里叶其人以及FFT算法的由来想必在这里不必多讲，在FFT算法广泛应用的今天，其强大的威力也早已被人所熟知。FFT算法是一种可以较高效率的计算离散傅里叶变换的算法，N点的DFT共需要N²次复数乘法和N(N-1)次复数加法，而利用FFT算法，任何一个N为2的整数幂（即N= 2^M）的DFT，都可以通过M次分解，最后成为2点的DFT来计算。M次分解构成了从x(n)到X(k)的M级迭代计算，每级由N/2个蝶形运算组成。完成一个蝶形计算需一次（复数）乘法和两次复数加法。因此，完成N点的时间抽选FFT计算的总运算量可以节省为M*N/2=log₂N*N/2次复数乘法和M*2*N/2= log₂N*N次复数加法。本文将以手算8点FFT变换，以及采用C和C++实现FFT变换算法，最后使用matlab里面的FFT来进行验算的方式带你领略FFT算法的博大精深。

用FFT手算8点DFT

本文采用的例子是余翻译的《数字信号处理》第四版（我们的教材）第165面的例题5.16.其原序列为

原序列

x[0]	x[1]	x[2]	x[3]	x[4]	x[5]	x[6]	x[7]
1	2	2	2	0	1	1	1

首先我们要明白,利用单个N点的DFT计算一个2N点的DFT的基本思想。大概就是利用形成偶序号样本x0[n]=x[2n]和奇序号样本x1[n]=x[2n+1]形成的两个长度为N/2的子序列的两个N/2点的DFT的加权和，来计算原来样本序列的DFT变换。例如，如果先计算N/2点的DFT的话，形式大致是如下所示：

                                                               

先计算两个N/2点的DFT变换，然后计算加权和。然而，每一个N/2点的序列又可以再分为2，即为原来的N/4点，依次类推，我们这里假设N是2^m，那么，最后总是可以分到只剩下2点的DFT变换。对于8点的DFT来说，N/4就已经是剩下2点了。对于N/4点，其计算的流图如下所示：

                                                       

因为两点的DFT已经不可再分解了，并且两点DFT可以很容易计算，其计算的流图如下：

                                                           

其中Wn是旋转因子，

                                                                          

上面的2点DFT运算也是FFT的基本运算模块，由其形状特点被称为蝶形运算，改进后的蝶形运算的基本模块流图如下：

                                                       

所以8点的DFT变换的流图可以进一步变为如下：

                                                         

对于8点的DFT，我们由Wn的计算公式可以得到如下的结果：

                                  

所以8点的整个FFT运算过程如下图所示：

             

由蝶形运算的基本模块公式一步步向后算，便可得到最终的结果，其结果和题目的答案是一致的。（其中小字体的表示旋转因子，横线上的表示该级蝶形运算的结果）我觉得如果我们自己手算过FFT变换的话，对其中的一些来龙去脉会比较清楚，也为后面的算法设计提供了基础。

FFT算法实现分析

我们现在来对FFT算法实现进行分析。FFT算法的实现主要有3点，分别是码位倒置，蝶形运算和旋转因子。

码位倒置

从上面的计算过程我们可以看到，FFT变换最开始的序列的顺序已经不是原始的序列的顺序了，其变换的对应关系如下图所示：

                                        

其对应的二进制序号如上图的有半部分所示，观察得到其码位的序号是反序的，即在进行FFT变换之前（或之后），需要对FFT的码位进行倒置。

假设一个输入序列为N点的，那么它的序号二进制数位数就是t=log2N（我们假设输入序列总是2的幂）。可以利用移位运算实现码位的倒置。假设输入为n2n1n0，则输出为n0n1n2。用一个变量j每次获得原序号的最低位，并不断右移使得低位变成高位，最终可以实现码位倒置的功能，码位倒置的程序可参考下面的代码：

/**
*reverse - 码位倒置
*@data:输入的序列
*/
void FFT::reverse(std::vector<complex>& data)//码位倒置函数
{
	unsigned int i, j, k;
	unsigned int t;
	complex temp;//临时交换变量
	for (i = 0;i<N;i++)//从第0个序号到第N-1个序号
	{//每个大循环实现一次序号交换
		k = i;//当前第i个序号
		j = 0;//存储倒序后的序号，先初始化为0
		for (t = 0;t<log2N;t++)//实现倒置 例如：0101 -> 1010
		{//每个小循环得到一个交换的序号
			j <<= 1;
			j |= (k & 1);//j左移一位然后或上k的最低位
			k >>= 1;    //k右移一位，次低位变为最低位
		}
		if (j>i)//如果倒序后大于原序数，就将两个存储单元进行交换，根据对称性，这样可以减少交换的次数
		{
			temp = data[i];
			data[i] = data[j];
			data[j] = temp;
		}
	}
}

蝶形运算

通过观察上面的蝶形运算图，我们可以总结出下面的一些规律。

对于8点的FFT，一共有3级蝶形运算

第一级有4组蝶形运算，每组有1个蝶形运算，其中每个蝶形运算的两个节点的距离为1（相邻） 

第二级有2组蝶形运算，每组有2个蝶形运算，其中每个蝶形运算的两个节点的距离为2

第三级有1组蝶形运算，每组有4个蝶形运算，其中每个蝶形运算的两个节点的距离为4

大致如下图所示：

   

可以推倒，如果是N点的FFT，则一共有log2N级的蝶形运算，第m级的蝶形运算，每个蝶形两节点的距离是L=2^(m-1)，第m级有N/2L组蝶形，每组有2^(m-1)个蝶形运算（即蝶形节点距离就是该组的蝶形个数）。因此在实现的时候，可以采用3组嵌套循环，第一层为log2N级蝶形运算，第二层为N/2L组蝶形，第三层为2^(m-1)即L个蝶形运算。

旋转因子

旋转因子其实就是所谓的加权。观察上面的旋转因子，对于8点的FFT，有如下的结果(由于旋转因子的可约性，其他地方的形式可能与这里不同，这里所有旋转因子均不约去)

第一级的旋转因子为W8[0](这里0是指数)

第二级的旋转因子为W8[0]，W8[2]

第三级的旋转因子为W8[0]，W8[1]，W8[2]，W8[3]。

通过推导，对于N点的FFT，有

第一级的旋转因子为WN[0]

第二级的旋转因子为WN[0],WN[N/4]，可以表示成WN[(N/2^2)*k],其中0<=k<2

第三级的旋转因子为WN[0],WN[N/8],WN[2N/8],WN[3N/8]，可以表示成WN[(N/2^3)*k],其中0<=k<4

第m级的旋转因子为WN[0],WN[1],……,WN[N/2-1]，可以表示成WN[(N/2^m)*k],其中0<=k<2^(m-1)

因此，对于不进行约操作的旋转因子，我们可以得到，第m级的第k个旋转因子为W[(N/2^m)*k]，其中0<=k<2^(m-1)，即第m级共有2^(m-1)个旋转因子。

因为我这里的C++代码主要跑在PC上，因此这里的旋转因子我采用运行计算的方式生成，利用欧拉公式将旋转因子的式子分解为三角函数的形式，可以利用下面的代码生成旋转因子表。

for (k = 0; k < N / 2; k++)//计算旋转因子表
	{
		temp.real = (float)cos(2 * PI / N * k);
		temp.imag = -(float)sin(2 * PI / N * k);
		WN_array.push_back(temp);
	}

另外，可以发现，我们上面整个过程中都只是用到了一半的旋转因子，因此上面的代码也只是生成前半部分的旋转因子。如果是在MCU上面跑的FFT的话，为了节省CPU的计算，建议先生成旋转因子表，然后复制到代码中即可。

蝶形运算的程序参考如下：

/**
 *fft - fft变换的外部接口
 *@data:原始的数据
 *@n:点数
 */
void FFT::fft(std::vector<complex>& data, unsigned int n)
{
	unsigned int i, j, k, l;
	complex top, bottom;//蝶形运算的上下两个部分
	complex temp;
	init(n);//初始化
	print_something();//查看一下
	reverse(data); //码位倒序
	for (i = 0;i < log2N;i++)//共log2N级
	{ //一级蝶形运算
		l = 1 << i;//l等于2的i次方
		for (j = 0;j < N;j += 2 * l)//每l个蝶形是一组，每级有N/2l组，j += 2 * l得到下一组蝶形运算的下标
		{   //一组蝶形运算
			for (k = 0;k < l;k++)//每组有l个
			{  //课本里改进后的蝶形运算，只做一次复数乘积
				temp = complex_mul(data[j + k + l], WN_array[N / (2 * l)*k]);//碟间距为l//每组的第k个蝶形	
				top = complex_add(data[j + k], temp);//上面
				bottom = complex_sub(data[j + k], temp);//下面
				data[j + k] = top;
				data[j + k + l] = bottom;
			}
		}
	}
}

FFT算法C++实现

C++实现的FFT算法仿照matlab中fft的一种接口，输入运算的数组以及点数N，从而进行FFT变换。采用一个vector来存放输入和输出序列，定义一个FFT类，该类的头文件如下所示：

/**
 *From zero to greatness
 *@author:LinChuangwei
 *@Email:979951191@qq.com
 *@date:11/27/2015 
 *@brief:FFT类头文件
 */
#ifndef _FFT_H_
#define _FFT_H_

#include <iostream>
#include <vector>

#define PI 3.1415926
//复数结构体，用于复数类型
typedef struct complex 
{
	float real;//实部
	float imag;//虚部
}complex;

//FFT 类
class FFT
{
public:
	void fft(std::vector<complex>& data,unsigned int n);//fft变换接口
	void print_something();//打印点数以及旋转因子表
private:
	unsigned int N;//N点的FFT
	unsigned log2N;
	std::vector<complex> WN_array;//旋转因子数组
	//内部接口
	void init(unsigned int n);//初始化函数
	complex complex_add(complex a, complex b);//复数加法
	complex complex_sub(complex a, complex b);//复数减法
	complex complex_mul(complex a, complex b);//复数乘法
	void reverse(std::vector<complex>& data);//码位倒置函数
};

#endif /*FFT.h*/

实现文件如下所示：

/**
 *From zero to greatness
 *@author:LinChuangwei
 *@Email:979951191@qq.com
 *@date:11/27/2015 
 *@brief:FFT实现文件
 */
#include "stdafx.h"
#include "fft.h"

 /**
 *fft - fft变换的外部接口
 *@data:原始的数据
 *@n:点数
 */
void FFT::fft(std::vector<complex>& data, unsigned int n)
{
	unsigned int i, j, k, l;
	complex top, bottom;//蝶形运算的上下两个部分
	complex temp;
	init(n);//初始化
	print_something();//查看一下
	reverse(data); //码位倒序
	for (i = 0;i < log2N;i++)//共log2N级
	{ //一级蝶形运算
		l = 1 << i;//l等于2的i次方
		for (j = 0;j < N;j += 2 * l)//每l个蝶形是一组，每级有N/2l组，j += 2 * l得到下一组蝶形运算的下标
		{   //一组蝶形运算
			for (k = 0;k < l;k++)//每组有l个
			{  //课本里改进后的蝶形运算，只做一次复数乘积
				temp = complex_mul(data[j + k + l], WN_array[N / (2 * l)*k]);//碟间距为l//每组的第k个蝶形	
				top = complex_add(data[j + k], temp);//上面
				bottom = complex_sub(data[j + k], temp);//下面
				data[j + k] = top;
				data[j + k + l] = bottom;
			}
		}
	}
}

/**
 *print_something - 打印旋转因子表等
 */
void FFT::print_something()
{
	std::vector<complex>::iterator it;
	std::cout << "进行FFT变换的点数："<< N << std::endl;
	std::cout << "旋转因子表如下(前面一半)：" << std::endl;
	for (it = WN_array.begin();it != WN_array.end();it++)
	{
		std::cout << "real:" << it->real << "	" << "imag:" << it->imag << std::endl;
	}
}

/**
*init - 初始化函数
*@n:点数
*主要完成旋转因子表的计算
*/
void FFT::init(unsigned int n)
{
	unsigned int k;
	complex temp;
	N = n;//初始化两个数
	log2N = (unsigned int)log2(N);
	for (k = 0; k < N / 2; k++)//计算旋转因子表
	{
		temp.real = (float)cos(2 * PI / N * k);
		temp.imag = -(float)sin(2 * PI / N * k);
		WN_array.push_back(temp);
	}
}

/**
*complex_add - 复数加法
*@a:输入的复数
*@b:输入的复数
*返回运算的复数结果
*/
complex FFT::complex_add(complex a, complex b)//复数加法
{//复数加法：实部+实部，虚部+虚部
	complex result;//用于保存一些临时的变量
	result.real = a.real + b.real;
	result.imag = a.imag + b.imag;
	return result;
}

/**
*complex_sub - 复数减法
*@a:输入的复数
*@b:输入的复数
*返回运算的复数结果
*/
complex FFT::complex_sub(complex a, complex b)//复数减法
{//复数减法：实部-实部，虚部-虚部
	complex result;//用于保存一些临时的变量
	result.real = a.real - b.real;
	result.imag = a.imag - b.imag;
	return result;
}

/**
*complex_mul - 复数乘法
*@a:输入的复数
*@b:输入的复数
*返回运算的复数结果
*/
complex FFT::complex_mul(complex a, complex b)//复数乘法
{//按复数乘法规则
	complex result;//用于保存一些临时的变量
	result.real = a.real*b.real - a.imag*b.imag;
	result.imag = a.real*b.imag + a.imag*b.real;
	return result;
}

/**
*reverse - 码位倒置
*@data:输入的序列
*/
void FFT::reverse(std::vector<complex>& data)//码位倒置函数
{
	unsigned int i, j, k;
	unsigned int t;
	complex temp;//临时交换变量
	for (i = 0;i<N;i++)//从第0个序号到第N-1个序号
	{//每个大循环实现一次序号交换
		k = i;//当前第i个序号
		j = 0;//存储倒序后的序号，先初始化为0
		for (t = 0;t<log2N;t++)//实现倒置 例如：0101 -> 1010
		{//每个小循环得到一个交换的序号
			j <<= 1;
			j |= (k & 1);//j左移一位然后或上k的最低位
			k >>= 1;    //k右移一位，次低位变为最低位
		}
		if (j>i)//如果倒序后大于原序数，就将两个存储单元进行交换，根据对称性，这样可以减少交换的次数
		{
			temp = data[i];
			data[i] = data[j];
			data[j] = temp;
		}
	}
}

测试的代码如下：

/**
 *From zero to greatness
 *@author:LinChuangwei
 *@Email:979951191@qq.com
 *@date:11/27/2015 
 *@brief:FFT算法测试程序
 */
#include "stdafx.h"
#include "fft.h"
int main()
{
	std::vector<complex> data = { {1,0}, {2,0}, {2,0}, {2,0} ,{0,0},{1,0},{1,0},{1,0} };
	std::vector<complex>::iterator it;
	FFT lcw_fft;//定义一个FFT类
	std::cout << "FFT变换前的序列：" << std::endl;
	for (it = data.begin();it != data.end();it++)
	{
		std::cout <<"real:" << it->real << "	" <<"imag:" << it->imag << std::endl;
	}
	lcw_fft.fft(data,8);//进行FFT变换
	std::cout << "FFT变换后的序列：" << std::endl;
	for (it = data.begin();it != data.end();it++)
	{
		std::cout << "real:" << it->real << "	" << "imag:" << it->imag << std::endl;
	}
    return 0;
}

改代码的输入测试数据也是上面的那道例题，输出采用实部和虚部分开输出的方式，其运行结果如下：

                                 
可以看到，其运算的结果是正确的。

FFT算法C实现

因为我和MCU的关系比较密切，最后可能很多的算法是要在MCU上面跑的，为了增加算法的可移植性，将上面的C++程序改为了C语言的版本。

C语言的实现这里没有实现多种点数的FFT变换（这里实现的是固定点数的FFT变换），因为建议是要先算好某个点数的旋转因子表，然后到时查表即可，这里假设我们已经学会了FFT的算法并且学会了如何改代码，因为这里实现的也是上面例子的代码（8点的FFT）

（移植时，一般只需修改旋转因子表，以及几个宏定义）

C语言实现的头文件如下：

/**
 *From zero to greatness
 *@author:LinChuangwei
 *@Email:979951191@qq.com
 *@date:11/27/2015 
 *@brief:fft算法C实现头文件
 */
#ifndef _FFT_H_
#define _FFT_H_

#define PI 3.141592
#define N 8 //点数
#define log2N 3

 //复数结构体，用于复数类型
typedef struct complex
{
	float real;//实部
	float imag;//虚部
}complex;

void fft(complex data[]);//FFT变换函数
void reverse(complex data[]);//码位倒置函数
complex complex_add(complex a, complex b);//复数加法
complex complex_sub(complex a, complex b);//复数减法
complex complex_mul(complex a, complex b);//复数乘法

#endif/*FFT.h*/

实现文件如下：

/**
 *From zero to greatness
 *@author:LinChuangwei
 *@Email:979951191@qq.com
 *@date:11/27/2015 
 *@brief:FFT变换C实现头文件
 */

#include "FFT.h"
 //旋转因子表
complex WN_array[N / 2] =
{//在嵌入式的处理器中，为节省CPU的计算，一般应先把旋转因子计算好
 //可用如下的语句计算,因为只用到前面一半，可计算一半的表即可
 /*	for (k = 0; k < N / 2; k++)//计算旋转因子表
 {
 WN_array[k].real = cos(2 * PI / N * k);
 WN_array[k].imag = -sin(2 * PI / N * k);
 }
 */
	{ 1,0 },{ 0.707107,-0.707107 },{ 2.67949e-08,-1 },{ -0.707107,-0.707107 }
};
 /**
 *fft - fft变换的外部接口
 *@data:原始的数据
 */
void fft(complex data[])
{
	unsigned int i, j, k, l;
	complex top, bottom;//蝶形运算的上下两个部分
	complex temp;
	reverse(data); //码位倒序
	for (i = 0;i < log2N;i++)//共log2N级
	{ //一级蝶形运算
		l = 1 << i;//l等于2的i次方
		for (j = 0;j < N;j += 2 * l)//每l个蝶形是一组，每级有N/2l组
		{   //一组蝶形运算
			for (k = 0;k < l;k++)//每组有L个
			{  //课本里改进后的蝶形运算，只做一次复数乘积
				temp = complex_mul(data[j + k + l], WN_array[N / (2 * l)*k]);//碟间距为l//每组的第k个蝶形			
				top = complex_add(data[j + k], temp);//上面
				bottom = complex_sub(data[j + k], temp);//下面
				data[j + k] = top;
				data[j + k + l] = bottom;
			}
		}
	}
}

/**
*complex_add - 复数加法
*@a:输入的复数
*@b:输入的复数
*返回运算的复数结果
*/
complex complex_add(complex a, complex b)//复数加法
{//复数加法：实部+实部，虚部+虚部
	complex result;//用于保存一些临时的变量
	result.real = a.real + b.real;
	result.imag = a.imag + b.imag;
	return result;
}

/**
*complex_sub - 复数减法
*@a:输入的复数
*@b:输入的复数
*返回运算的复数结果
*/
complex complex_sub(complex a, complex b)//复数减法
{//复数减法：实部-实部，虚部-虚部
	complex result;//用于保存一些临时的变量
	result.real = a.real - b.real;
	result.imag = a.imag - b.imag;
	return result;
}

/**
*complex_mul - 复数乘法
*@a:输入的复数
*@b:输入的复数
*返回运算的复数结果
*/
complex complex_mul(complex a, complex b)//复数乘法
{//按复数乘法规则
	complex result;//用于保存一些临时的变量
	result.real = a.real*b.real - a.imag*b.imag;
	result.imag = a.real*b.imag + a.imag*b.real;
	return result;
}

/**
*complex_add - 复数乘法
*@data[]:输入的序列
*/
void reverse(complex data[])//码位倒置函数
{
	unsigned int i, j, k;
	unsigned int t;
	complex temp;//临时交换变量
	for (i = 0;i<N;i++)//从第0个序号到第N-1个序号
	{//每个大循环实现一次序号交换
		k = i;//当前第i个序号
		j = 0;//存储倒序后的序号，先初始化为0
		for (t = 0;t<log2N;t++)//实现倒置 例如：0101 -> 1010
		{//每个小循环得到一个交换的序号
			j <<= 1;
			j |= (k & 1);//j左移一位然后或上k的最低位
			k >>= 1;    //k右移一位，次低位变为最低位
		}
		if (j>i)//如果倒序后大于原序数，就将两个存储单元进行交换，根据对称性，这样可以减少交换的次数
		{
			temp = data[i];
			data[i] = data[j];
			data[j] = temp;
		}
	}
}

测试的代码如下：

/**
 *From zero to greatness
 *@author:LinChuangwei
 *@Email:979951191@qq.com
 *@date:11/27/2015 
 *@brief:FFT测试程序
 */
#include "FFT.h"
#include <stdio.h>

int main()
{
	complex data[] = { { 1,0 },{ 2,0 },{ 2,0 },{ 2,0 } ,{ 0,0 },{ 1,0 },{ 1,0 },{ 1,0 } };
	int i = 0;
	printf("FFT变换前：
");
	for (i = 0;i < N;i++)
	{
		printf("real:%f	imag:%f
", data[i].real, data[i].imag);
	}
	fft(data);
	printf("FFT变换后：
");
	for (i = 0;i < N;i++)
	{
		printf("real:%f	imag:%f
", data[i].real, data[i].imag);
	}
	return 0;
}

运行的结果和上图是类似的。

matlab验算

可以使用一个简短的matlab代码验算一下我们手算，以及程序的正确性。

function lcw_fft
close all;
data=[1 2 2 2 0 1 1 1];
x = fft(data,8);
real_ = real(x)
imag_ = imag(x)

运算的结果如下：

                             
对比一下，可以知道我们的结果是正确的。之后将会使用FFT算法做一些有趣的东西，敬请期待。

由于时间仓促，本文的推导及程序可能存在错误，希望各位提出批评和建议。