原标题：初学者看狂野的树叶就知道了。 详细了解FFT算法的原理
我相信网上现在有很多关于FFT的教程。 我也参考了很多网络教程，但我觉得不太容易懂。 基本上是研究FFT，以编程的形式实现后。 我决定对FFT进行通俗易懂的说明。 因此，我会在今后的记述中尽量通俗易懂地仔细说明。
第一个知道傅立叶变换的时候，我沉迷于音乐的频谱跳动而无法自拔，于是我想做一个音乐的频谱显示屏。 在查阅了很多资料之后，我知道了傅立叶变换和FFT。 当然这是很久以前的事了，最近抽出时间调试，自制了FFT的算法。 虽然不能说速度上的好处，但是我个人认为是容易理解的实现方法。 实际的VC模拟，与STM32的行驶也成功了。 所以我在这里和大家分享。
FFT
算法原理
首先，为了能够进行FFT，需要了解DFT。 因为两者之间有本质上相同的东西。 在此之前，先列举离散傅立叶变换对(DFT )。
但是FFT之所以被称为快速傅立叶变换，是利用了以下性质()
先把这三个公式放在这里，接下来用。
这些特性使得可以将长DFT计算分解为几个短序列的DFT计算的组合，从而减少计算量。 为了便于理解，我们使用了按时间提取的快速傅立叶变换(DIT-FFT )方法。
首先，将一个序列x(n )分成两部分。
在以下情况下：
设n为2的整数次方，即N=2^M。
按奇偶校验对x(n )进行分组：
DFT也分为两组进行预算。
此时，在我们上面列举的三个性质中，约定性很有用。
让我们回顾一下：
该公式如下所示。
x1(k )和x2 ) k )都是长度为N/2的序列x1 ) k )和x2 ) k )的N/2点的离散傅立叶变换。
其中：
这样，一个n点的DFT被分解为两个N/2的DFT。 但是，x1(k )和x1(k )只有N/2分。 也就是说，公式(1-1)只是DFT的前半部分。 DFT的后半部分被要求能够利用其对称性求出后半部分。
式(1-1)和式(1-2)可以用蝶形信号流程图符号表示。 图：
以N=8为例，可以用下图表示。
那么，通过这样分解，每1个N/2点DFT，只需要(N^2)/4次的复数乘法，就可以明显地节约运算量。
那么，接下来类推，是继续按奇偶校验分解已经出现的x1(k )和x1(k )，还是与上边相同的方法？ 在百度上也能找到的这张照片出山了。
对于这张图像，我想强调的是二次蝶形运算的时候。 之后不应该变成那样吗？ 为什么会变成这样？ 这个问题之前我烦恼了一会儿，但别忘了。 前者的展开是N/2，因为此时n是奇偶校验分离的，是从根据的可约束性得出的，在此不能混淆。
没有必要提到计算效率。 m级计算总共需要：
FFT
c语言实现
FFT算法c语言的实现，网络方法层出不穷，除了本人懒惰(不擅长看别人的程序)外，还有自给自足的衣食住原则，所以我自己也写了个人容易理解的程序。 而且，为了帮助读者理解，我特意尽量减少了库函数的使用。 一些基本函数是自己写的。 但是，作为FFT算法已经足够了。 目前，该程序只能处理2^N的数据。 理论上，如果2^N不够，可以进行在后面的数列中添加0的操作。 当然，因为简约语句也不能实现，但主要功能如下开源。
左右滑动查看代码
/*
功能：对input内的数据进行快速傅立叶变换并输出
*/
#包含
#包含
#define FFT_LENGTH 8
double input [ FFT _ length ]={ 1，1，1，1，1，1，1，1 }；
定义结构完成1 {//多个结构
双现实； //实部
双图像； //虚部
(；
将input的实数结果保存在复数中
结构完成1 result _ dat [8]；
/*
虚数乘法
*/
结构完成1 con _ complex (结构完成1 a，结构完成1 b ) {
struct complex1 temp;temp.real=(a.real*b.real)-(a.image*b.image);
temp.image=(a.image*b.real)+(a.real*b.image);
return temp;
}
/*
简单的a的b次方
*/
int mypow(int a,int b){
int i,sum=a;
if(b==0)return 1;
for(i=1;i
sum*=a;
}
return sum;
}
/*
简单的求以2为底的正整数
*/
int log2(int n){
unsigned i=1;
int sum=1;
for(i;;i++){
sum*=2;
if(sum>=n)break;
}
return i;
}
/*
简单的交换数据的函数
*/
void swap(struct complex1 *a,struct complex1 *b){
struct complex1 temp;
temp=*a;
*a=*b;
*b=temp;
}
/*
dat为输入数据的数组
N为抽样次数 也代表周期 必须是2^N次方
*/
void fft(struct complex1 dat[],unsigned char N){
/*最终 dat_buf计算出 当前蝶形运算奇数项与W 乘积
dat_org存放上一个偶数项的值
*/
struct complex1 dat_buf,dat_org;
/* L为几级蝶形运算 也代表了2进制的位数
n为当前级蝶形的需要次数 n最初为N/2 每级蝶形运算后都要/2
i j为倒位时要用到的自增符号 同时 i也用到了L碟级数 j是计算当前碟级的计算次数
re_i i_copy均是倒位时用到的变量
k为当前碟级 cos(2*pi/N*k)的 k 也是e^(-j2*pi/N)*k 的 k
*/
unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;
//经过观察，发现每级蝶形运算需要N/2次运算，共运算N/2*log2N 次
unsigned char fft_counter=0;
//在此要进行补2 N必须是2^n 在此略
//蝶形级数 (L级)
L=log2(N);
//计算每级蝶形计算的次数(这里只是一个初始值) 之后每次要/2
//n=N/2;
//对dat的顺序进行倒位
for(i=1;i
i_copy=i;
re_i=0;
for(j=L-1;j>0;j--){
//判断i的副本最低位的数字 并且移动到最高位 次高位 ..
//re_i为交换的数 每次它的数字是不能移动的 并且循环之后要清0
re_i|=((i_copy&0x01)<
i_copy>>=1;
}
swap(&dat[i],&dat[re_i]);
}
//进行fft计算
for(i=0;i
fft_flag=1;
fft_counter=0;
for(j=0;j
if(fft_counter==mypow(2,i)){ //控制隔几次，运算几次
fft_flag=0;
}else if(fft_counter==0){ //休止结束，继续运算
fft_flag=1;
}
//当不判断这个语句的时候 fft_flag保持 这样就可以持续运算了
if(fft_flag){
dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));
dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));
dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);
//计算 当前蝶形运算奇数项与W 乘积
dat_org.real=dat[j].real;
dat_org.image=dat[j].image; //暂存
dat[j].real=dat_org.real+dat_buf.real;
dat[j].image=dat_org.image+dat_buf.image;
//实部加实部 虚部加虚部
dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;
dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;
//实部减实部 虚部减虚部
k++;
fft_counter++;
}else{
fft_counter--; //运算几次，就休止几次
k=0;
}
}
}
}
void main{
int i;
//先将输入信号转换成复数
for(i=0;i
result_dat[i].image=0;
//输入信号是二维的，暂时不存在复数
result_dat[i].real=input[i];
//result_dat[i].real=10;
//输入信号都为实数
}
fft(result_dat,FFT_LENGTH);
for(i=0;i
input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);
//取模
printf("%lfn",input[i]);
}
while(1);
}
这个程序中input这个数组是输入信号，在这里只模拟抽样了8次，输出的数据也是input，如果想看其它序列的话，可以改变FFT_LENGTH 的值以及 input里的内容，程序输出的是实部和虚部的模，如果单纯想看实部或者虚部的幅度的话，请自行修改程序~
↓↓↓↓点击阅读原文，查看更多新闻返回搜狐，查看更多
责任编辑：