十分钟搞懂基-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算法思路分析：

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

2. 第二部：蝶形结构分析（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;
        }
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13

// 单个蝶形运算单元对应的旋转因子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;

        }
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

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

4.自己编程实现FFT

1. 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:N
    din_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<j
        t = 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，置0
            j=j-k;
            k=k/2;%变为010，准备比较第二位的1，循环
        end
        if j<k
            j=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<j
        t = 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，置0
            j=j-k;
            k=k/2;%变为010，准备比较第二位的1，循环
        end
        if j<k
            j=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-1
        for j=1:(2^(m-1))  % 注意：j是从1开始，实际用到应该j-1
            len = 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);  % 蝶形山支输出
                        
        end
    end
end

self_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;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

运行结果对比：

2. C语言编程实现

#include "stdio.h"
#include "math.h"
#include "stdlib.h"

#define PI 3.1415926

int  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_EN
   			printf("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_EN
   			printf("w_real : %f\r\n",w_real);
   			printf("w_imag : %f\r\n",w_imag);
   			#endif
   			
   	        temp_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_EN
           	printf("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;	
} 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

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

5.总结

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