用FFT手算8点DFT

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

首先我们要明白,利用单个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的幂）。可以利用移位运算实现码位的倒置。假设输入为n2n1 n0 ，则输出为 n0 n1n2。用一个变量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 << "\t" << "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 << "\t" <<"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 << "\t" << "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变换前：\n");
	for (i = 0;i < N;i++)
		printf("real:%f\timag:%f\n", data[i].real, data[i].imag);
	fft(data);
	printf("FFT变换后：\n");
	for (i = 0;i < N;i++)
		printf("real:%f\timag:%f\n", 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算法做一些有趣的东西，敬请期待。

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