快速傅里叶变换算法探幽
学完《信号与系统》之后发现什么还是不会？
感觉《数字信号处理》和《信号与系统》差不多？
怎么办？期末快到了……是不是感觉《数字信号处理》学了和没学一样？
对的，就是这种感觉。大一的时候用FFT算法做过音乐频谱，然而当时对里面的算法并不是很了解，之前面某俱乐部的时候被面试官问到，因此……虽然还是……。学完DSP之后感觉没自己来尝试下FFT算法真是人生一大憾事。于是赶着期末之前当做是复习下（虽然并不在考试范围）。
傅里叶其人以及FFT算法的由来想必在这里不必多讲，在FFT算法广泛应用的今天，其强大的威力也早已被人所熟知。FFT算法是一种可以较高效率的计算离散傅里叶变换的算法，N点的DFT共需要N2次复数乘法和N(N-1)次复数加法，而利用FFT算法，任何一个N为2的整数幂（即N= 2M）的DFT，都可以通过M次分解，最后成为2点的DFT来计算。M次分解构成了从x(n)到X(k)的M级迭代计算，每级由N/2个蝶形运算组成。完成一个蝶形计算需一次（复数）乘法和两次复数加法。因此，完成N点的时间抽选FFT计算的总运算量可以节省为M*N/2=log2N*N/2次复数乘法和M*2*N/2= log2N*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);
}蝶形运算的程序参考如下：
/**
*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算法做一些有趣的东西，敬请期待。
由于时间仓促，本文的推导及程序可能存在错误，希望各位提出批评和建议。
