十分钟搞懂基-2 FFT原理及编程思想

0.写在最前

写本文的目的一是为了帮人理清FFT算法思路，二是有几个疑问（在5总结部分提到）希望得到解答。看懂本文的基础：至少听说过、简单了解过傅里叶变换、离散傅里叶变换（DFT）、基于时间抽取的基2FFT算法等理论，本文的重点在于介绍编程思想、对于一些理论只做简要介绍，笔者认为：理论部分看下专业的书籍是最靠谱的，推荐的书籍有郑君里的《信号与系统-下》，或者其他一些相关的数字信号处理教材。

1.DFT（离散傅里叶变换）

DFT的由来不再赘述，理解其由来抓住一句话，我们的计算机处理有限长度的数字信号（不管是时域上还是频域上）。从书中截取离散傅里叶变换定义，注意DFT是针对有限长的序列（例如AD采样若干个点构成的序列），且DFT结果是复数（a+bj）的形式。

旋转因子：

2.FFT（快速傅里叶算法）

FFT是计算DFT的一种快速计算方法，主要是利用了旋转因子的周期性、对称性、可约性的特点，能够实现“新点旧算”，通过减少重复计算来减少计算量。DFT和FFT的性能对比在很多书籍、博客都有详细对比，这里不再赘述。下面简单介绍基于时间抽取基-2FFT算法原理：

基2算法的原理就是将N=2^M个点的FFT一步步按照上述原理进行分解，直到分解到2点的FFT。下面是FFT算法思路分析：

第一步：码位倒置
在进行FFT之前需要对数据顺序做一下的调整，这样做的目的是为了保证输出数据是顺序的，否则输入顺序的输出就是乱序的（但是有规律的），这个调整称为“码位倒置”，规律如下：

那么8点FFT8个输入的顺序地址为：000b~111b，调整之后顺序如上表第三列所示。
从上面看，规律是很简单的：1.变换前后的二进制地址最高位和最低位互换、次高位和次低位互换、…；2.变换之后的地址之间的关系：下个地址是在前一个地址的最高位加1并有进位的结果，例如100b就是000b最高位加1得到的，011b是101b最高位加1产生进位得到的。
规律很简单，但是C语言程序实现起来并不简单，因为C语言不支持位操作，这里介绍一种常用的“雷德算法”去实现这个排序，这里也不再赘述，在后面编程是直接使用的。
第二部：蝶形结构分析（8点为例）

1.上图中同一种颜色的结构称之“蝶形运算单元”；

2.上述两个竖的红色虚线分出三块区域分别为：第一级、第二级、第三级，对于N=2^M
点时间抽取基-2 FFT共有M级（8=2^3点FFT就有3级）；
3.每一级共有N/2个蝶形运算单元，例如8点FFT每一级有4个蝶形运算单元，那么M级共有M*N/2个蝶形运算单元。
4.每一个蝶形运算单元：包含两个输入、一个旋转因子、两个输出，也就是要进行一次蝶形运算就得知道两个输入和旋转因子的值，那么就需要找到每个蝶形运算单元输入数据之间的规律、旋转因子的规律，这是编程必须要搞懂的。

3.FFT算法编程思想

对于N = 2^M 点FFT,基于基2算法C语言实现过程:
整体思路:
1)先分级，即上面提到的。
2)再分组,这里组的意思:是指蝶形(直观上)有交叉的为一组,否则就为不同组，例如下图中绿色方框的为一组，第一级有4组、第二级有2组、第三级有1组。
3)最后分蝶形运算单元个数,是指每一组包含的蝶形运算单元的个数，例如第一级每一组包含一个蝶形单元，第二级每组包含2个蝶形单元，第三级每组包含1个蝶形单元。

下面定义几个循环变量：当前级数m（1,2，…，M）,当前组数i，当前某一组的第几个蝶形单元j，那么需要搞清楚下面几个关键信息：
1)N = 2^M 点FFT的总级数:M
2)每一级所包含的蝶形运算单元: N/2=2^(M-1)
3)每一级所包含的组数: N/2^m = 2^(M-m),m表示当前级数(m=1,2,…,M)
4)每一组所包含的蝶形运算单元: 2^(m-1),m表示当前级数(m=1,2,…,M)
5)每一组所包含的蝶形运算点数 = 每一组所包含的蝶形运算单元x2 = 2x2^(m-1),m表示当前级数(m=1,2,…,M)
6)每个蝶形运算单元的两个输入之间的距离: 2^(m-1)
7)每个蝶形运算单元需要2个输入数据 + 1个旋转因子完成计算
8)每一级需要N/2=2^(M-1)
，个蝶形运算单元,即需要N/2=2^(M-1)个旋转因子
9)每一级旋转因子的种类:2^(m-1),m表示当前级数(m=1,2,…,M)
10）还需要了解复数的乘法和加法，这个比较简单这里不赘述

熟知上述参数之后，接下来开始核心编程思路，下面是蝶形运算单元，其实只要搞清每个蝶形单元的两个输入a,b和W，根据复数乘法和加法即可得到输出A、B：

// 单个蝶形运算单元两个输入点的地址
for(m=1;m<=M;m++)  // 以级数作为循环,从1~M
{for(i=0;i<(2^(M-m) );i++)  // 以组数作为循环,从0~(2^(M-m) - 1){for(j=0;j<(2^(m-1) );j++)  // 以每组蝶形数作为循环,从1~((2^(m-1)) - 1){len = 2^(m-1);  // 蝶形距离:蝶形两个输入数据下标的距离第1个输入点地址a = i*当前组蝶形运算点数 + j = i*2^m + j; 第2个输入点地址b = 第1个输入点地址 + 蝶形距离 = 第1个输入点地址 + len;}}
}

// 单个蝶形运算单元对应的旋转因子W,需要知道W下标N和W上标kn
for(m=1;m<=M;m++)  // 以级数作为循环,从1~M
{for(i=0;i<(2^(M-m) );i++)  // 以组数作为循环,从0~(2^(M-m) - 1){for(j=0;j<(2^(m-1) );j++)  // 以每组蝶形数作为循环,从1~((2^(m-1)) - 1){上标 = j;   // kn值下标 = 2^m; // N旋转因子实部 = cos(2π*上标/下标) = cos(2π*j/2^m);旋转因子虚部 = -sin(2π*上标/下标) = -sin(2π*j/2^m);}}
}// 单个蝶形运算单元进行计算,原位计算即计算的两个输出存储到对应的输入地方
for(m=1;m<=M;m++)  // 以级数作为循环,从1~M
{for(i=0;i<(2^(M-m) );i++)  // 以组数作为循环,从0~(2^(M-m) - 1){for(j=0;j<(2^(m-1) );j++)  // 以每组蝶形数作为循环,从1~((2^(m-1)) - 1){tmp.real = x2.real*w.real - x2.imag*w.imag;  // 旋转因子和输入第二个数进行复数乘法,复数乘法得到实部tmp.imag = x2.real*w.imag + x2.imag*w.real;  // 旋转因子和输入第二个数进行复数乘法,复数乘法得到虚部// 计算蝶形下支x2.real  = x1.real - tmp.real;x2.imag  = x1.imag - tmp.imag;// 计算蝶形上支x1.real  = x1.real + tmp.real;x1.imag  = x1.imag + tmp.imag;}}
}

注意：这里分开写的目的是为了帮助大家一步一步的分析，实际编程需要写到同一个for循环里面去的，所以对序列进行时间抽取基-2FFT算法编写思路：
1.顺序变换（初始化）；
2.理清分级、对每一级分组、对每一组分蝶形运算单元个数；
3.根据级数、组数、每一组蝶形个数三层for循环即可搞定。

4.自己编程实现FFT

MATLAB编程实现

%%-----------------------------------FFT main-----------------------------
clc;
clear all;
close all;N = 1024; % FFT计算点数
M = log2(N); % FFT最大级数
fs = 100; %采样率,100Hz%产生输入数据的实部
n=0:N-1;
t=n/fs;
din_real = 1.0*sin(2*pi*20*t); %sin信号,幅值为1.0,频率f=20Hz
din_real_tmp = din_real;%输入数据的虚部全部置0
for i=1:Ndin_imag(i) = 0.0; % 虚部初始化为0
end%-----------------------用matlab自带fft函数进行计算-------------------------
fft_result = fft(din_real_tmp,N);
fft_result_amp = abs(fft_result);
figure(1);
f=n*fs/N;
%plot(f(1:N/2),fft_result_amp(1:N/2) );%画频谱图，幅值和频率关系(直流分量可能有点问题)
plot(fft_result_amp);%频谱和点数图
title('matlab 自带FFT函数计算结果幅度');
%--------------------------------------------------------------------------%----------------------------自己根据基-2算法实现FFT------------------------
%第一步：调整输入数据实部的顺序，其实这里实部和虚部可以放在一起排序，代码更简洁，执行速度也更快，这里为了简单操作分开排序了
data_len = length(din_real);
j=data_len/2;  %数组半长for i=1:data_len/2-1  %这里实现了奇偶前后分开排序 ,if i<jt = din_real(j+1);din_real(j+1) = din_real(i+1);%交换din_real(i+1) = t;end%求下一个倒序数k = data_len/2;%数组半长while(j>=k)%j为下一个倒序数，比较100的最高位1，如果为1，置0j=j-k;k=k/2;%变为010，准备比较第二位的1，循环endif j<kj=j+k;%找到为0 的一位，补成1，j就是下一个倒序数end
end
new_ord_data_real = din_real; % 调试查看,不用于计算%第二步：调整输入数据虚部的顺序
data_len = length(din_imag);
j=data_len/2; %数组半长for i=1:data_len/2-1 %这里实现了奇偶前后分开排序 ,if i<jt = din_imag(j+1);din_imag(j+1) = din_imag(i+1);%交换din_imag(i+1) = t;end%求下一个倒序数k = data_len/2;%数组半长while(j>=k)%j为下一个倒序数，比较100的最高位1，如果为1，置0j=j-k;k=k/2;%变为010，准备比较第二位的1，循环endif j<kj=j+k;%找到为0 的一位，补成1，j就是下一个倒序数end
end
new_ord_data_imag = din_imag; % 调试查看,不用于计算%第三步：蝶形计算
for m=1:M   % 注意：m是从1开始for i=1:(2^(M-m))  % 注意：i是从1开始，实际用到应该i-1for j=1:(2^(m-1))  % 注意：j是从1开始，实际用到应该j-1len = 2^(m-1);  %两个输入数据之间的蝶形距离% 两个输入数据下标索引addr1 = (i-1)*2^m + (j - 1) + 1; %加1是因为matlab数组首地址从1开始寻址addr2 = addr1 + len;% 蝶旋转因子的上下标shang_biao = j - 1;xia_biao   = 2^m;      % 计算旋转因子的实部和虚部w_angle = 2.0*pi*shang_biao/xia_biao;  % 弧度w_real = cos(w_angle);w_imag = -1.0*sin(w_angle);% 计算第二个输入数据与旋转因子的实部和虚部%tmp_real = din_real(addr2)*w_real - din_imag(addr2)*w_imag; % 旋转因子和输入第二个数进行复数乘法,复数乘法得到实部%tmp_imag = din_real(addr2)*w_imag + din_imag(addr2)*w_real; % 旋转因子和输入第二个数进行复数乘法,复数乘法得到虚部      [tmp_real,tmp_imag] = complex_mul(din_real(addr2),din_imag(addr2),w_real,w_imag); % 第二个数与旋转因子的做复数乘法% 原位计算[din_real(addr2),din_imag(addr2)]  = complex_add(din_real(addr1),din_imag(addr1),-1.0*tmp_real,-1.0*tmp_imag);  % 蝶形下支输出[din_real(addr1),din_imag(addr1)]  = complex_add(din_real(addr1),din_imag(addr1),tmp_real,tmp_imag);  % 蝶形山支输出endend
endself_fft_amp = sqrt(din_real.^2 + din_imag.^2);
figure(2);
%plot(f(1:N/2),self_fft_amp(1:N/2));;%画频谱图，幅值和频率关系(直流分量可能有点问题)
plot(self_fft_amp);;%画幅值和点数关系
title('自己编写FFT函数计算结果幅度');%%-----------------------------------complex_mul function-----------------------------
% 两个复数乘法函数,输出计算结果的实部和虚部
function [out_real,out_imag] = complex_mul(d0_real,d0_imag,d1_real,d1_imag)out_real = d0_real*d1_real - d0_imag*d1_imag;out_imag = d0_real*d1_imag + d0_imag*d1_real;%%-----------------------------------complex_mul function-----------------------------
% 两个复数加法函数,输出计算结果的实部和虚部
function [out_real,out_imag] = complex_add(d0_real,d0_imag,d1_real,d1_imag);out_real = d0_real + d1_real;out_imag = d0_imag + d1_imag;

运行结果对比：

C语言编程实现

#include "stdio.h"
#include "math.h"
#include "stdlib.h"#define PI 3.1415926int  N = 256;  // FFT计算点数// 定义复数结构体
typedef struct
{double real;double imag;
}complex;// 输入数据重排函数定义
void Fft_Data_ReOrder(int N,complex Data[N])
{complex temp;unsigned short i=0,j=0,k=0;double t;for(i=0;i<N;i++){k = i;j = 0;t = (log(N)/log(2));while( (t--)>0 ){j  = j<<1;j |= (k & 1);k  = k>>1;}if(j>i){temp = Data[i];Data[i] = Data[j];Data[j] = temp;}}
}// FFT计算函数定义 ,N为FFT计算点数，Data[N]为计算输入和计算结果输出（原位计算）
void Fft_Calculate(int N,complex Data[N])
{ int M = 0;int addr0 = 0;   // 蝶形第一个数据下表索引int addr1 = 0;  // 蝶形第二个数据下表索引int m,i,j,len,shang_biao,xia_biao; float temp_real = 0.0;float temp_imag = 0.0;float w_real = 0.0;  // 旋转因子实部 float w_imag = 0.0;  // 旋转因子虚部    M = log2(N);  // 最大级数 for(m=1;m<=M;m++) % 以级数作为循环{for(i=0;i<( (int)pow(2,M-m) );i++)  // 以组数作为循环，这里的2^(M-m)除了通过pow函数可以通过将1左移M-m也是可以的{for(j=0;j<( ( (int)pow(2,m-1) ) );j++)   // 以每组蝶形个数作为循环{len = (int)pow(2,m-1);addr0 = i*(int)pow(2,m) + j;  // 蝶形运算单元上支输入数据地址 addr1 = addr0 + len;          // 蝶形运算单元下支输入数据地址 shang_biao = j;               // 旋转因子的上标 xia_biao   = (int)pow(2,m);   // 旋转因子的下标 // 调试信息#ifdef DEBUG_PRINTF_ENprintf("m : %d\r\n",m);printf("i : %d\r\n",i);printf("j : %d\r\n",j);printf("len : %d\r\n",len);printf("addr0 : %d\r\n",addr0);printf("addr1 : %d\r\n",addr1);printf("shang_biao : %d\r\n",shang_biao);printf("xia_biao : %d\r\n",xia_biao);#endif              w_real = cos(2.0*PI*shang_biao/xia_biao);       // 计算旋转因子的实部 w_imag = -1.0*sin(2.0*PI*shang_biao/xia_biao);  // 计算旋转因子的实虚部 #ifdef DEBUG_PRINTF_ENprintf("w_real : %f\r\n",w_real);printf("w_imag : %f\r\n",w_imag);#endiftemp_real = Data[addr1].real * w_real - Data[addr1].imag * w_imag;  // 旋转因子和输入第二个数进行复数乘法,复数乘法得到实部temp_imag = Data[addr1].real * w_imag + Data[addr1].imag * w_real;  // 旋转因子和输入第二个数进行复数乘法,复数乘法得到虚部#ifdef DEBUG_PRINTF_ENprintf("temp_real : %f\r\n",temp_real);printf("temp_imag : %f\r\n",temp_imag);#endif// 输出下支计算结果 Data[addr1].real  = Data[addr0].real - temp_real;Data[addr1].imag  = Data[addr0].imag - temp_imag;// 输出上支计算结果               Data[addr0].real  = Data[addr0].real + temp_real;Data[addr0].imag  = Data[addr0].imag + temp_imag;}}}
}// 复数取模函数
float Complex_Mod(complex Data)
{float tmp;tmp =  pow(Data.real,2) + pow(Data.imag,2);tmp = sqrt(tmp);return tmp;
} // FFT主函数
int main(void)
{complex data_complex[N];int i;float amp;// 生成数据    for(i=0;i<N;i++){data_complex[i].real = i;data_complex[i].imag = 0;}Fft_Data_ReOrder(N,data_complex);  // 数据顺序重排 Fft_Calculate(N,data_complex);     // 数据计算 printf("--------FFT Calculate Result is---------\r\n");for(i=0;i<N;i++){    amp = Complex_Mod(data_complex[i]);  // 计算每个复数点模长 printf("%f + j%f   amp:%f\r\n",data_complex[i].real,data_complex[i].imag,amp);}printf("--------FFT Calculator Finishing!--------\r\n");return 0;
}

C语言计算结果，太多了，未完全显示：
MATLAB 自带FFT函数计算结果（部分）：
MATLAB 自己编写FFT函数计算结果（部分），第一列为实部，第二列为虚部：

5.总结

目前笔者感觉思路和代码上没有问题，但从计算结果来看还是存在以下几个关键问题，希望有人能解答以下疑惑
1.自己在MATLAB编写的fft函数和FFT自带的函数运算结果相比，频谱曲线没有那么光滑，主频点的幅值有差异；
2.自己在MATLAB编写的fft函数运行的时间要比FFT自带的函数运行时间长很多，自己编写的fft函数计算65536就要10多秒