通俗细致谈FFT以及FFT算法的基本实现

算法&人工智能
相信网上现在有很多关于FFT的教程，我曾经也参阅了很多网上的教程，感觉都不怎么通俗易懂。在基本上的研究FFT，并且通过编程的形式实现之后。我决定写一篇通俗易懂的关于FFT的讲解。因此我在接下来的叙述中尽量非常通俗细致的讲解。

 

本人最早知道傅里叶变换的时候是沉迷于音乐的频谱跳动无法自拔，当时就很想做一个音乐频谱显示器。搜阅了很多资料之后，才了解到傅里叶变换，和FFT。当然这都是以前的事情了，经过了系统的学习+2个星期的研究，自制了一个FFT的算法，不敢说速度上的优势，但是个人认为是一种通俗易懂的实现方法。经过实际的VC++模拟实验、和STM32跑的也很成功。

 

首先，要会FFT，就必须要对DFT有所了解，因为两者之间本质上是一样的。在此之前，先列出离散傅里叶变换对(DFT):

但是FFT之所以称之为快速傅里叶变换，就利用了以下的几个性质（重中之重！）

         先把这仨公式放到这，接下来会用到。

 

根据这几个特性，就可以将一个长的DFT运算分解为若干短序列的DFT运算的组合，从而减少运算量。在这里，为了方便理解，我就用了一个按时间抽取的快速傅里叶变换（DIT-FFT）的方法。

 

首先，将一个序列x(n)一分为二

 

对于,k=0,1,…N-1   设N是2的整次幂 也就是N=2^M
将x(n)按照奇偶分组

 

x(2r)=x1(r)
x(2r+1)=x2(r)                     r=0,1,…..(N/2)-1
那么将DFT也分为两组来预算

            (第一项是偶  第二项是奇)

 

         这个时候，我们上边提的三个性质中的可约性就起到作用了：
回顾一下：
那么这个式子就可以化为：
  (式1-1)
则X1(k)和X2(k)都是长度为N/2的序列 x1(k)和x2(k)的N/2点的离散傅里叶变换
其中：
  K=0,1,2…N/2-1
至此，一个N点的DFT就被分解为2个N/2的DFT。但是X1(k),和X2(k)只有N/2个点，也就是说式子（1-1）只是DFT前半部分。要求DFT的后半部分可以利用其对称性求出后半部分为：
(式1-2)
那么式（1-1）和（1-2）就可以专用一个蝶形信号流图符号来表示。如图：

以N=8为例，可以用下图表示：

通过这样的分解，每一个N/2点DFT只需(N^2)/4次复数相乘计算，明显的节省了运算量。
以此类推，继续将已经得出的X1(k)和X2（k）按照奇偶分解，还是和上边一样的方法。那么就得出了百度上都可以找到的一大推的这个图片了。  （笑）

对于这张图片我想强调的一点就是第二阶蝶形运算的时候，WN0之后不应该是WN1吗，为什么是W2N了，这个问题之前困扰了我一段时间，但是不要忘了，前者的    W0N的展开是W0N/2  因为此时N已经按照奇偶分开了，所以是N/2  而W2N/2也就是  W2N是根据可约性得出的，在这里不能混淆.。

 

对于运算效率就不用多提了

 

以上就是FFT算法的理论内容了，接下来就是用C语言对这个算法的实现了，对于FFT算法C语言的实现，网上的方法层出不穷，介于本人比较懒（懒得看别人的程序），再加上自给自足丰衣足食的原则，我自己也写了一个个人认为比较通俗易懂的程序，并且为了帮助读者理解，我特意尽量减少了库函数的使用，一些基本的函数都是自己写的（难免有很多BUG），但是作为FFT算法已经够用了，目前这个程序只能处理2^N的数据，理论上来讲如果不够2^N的话可以在后面数列补0这种操作为了简约我也就没有实现，但是主要的功能是具备的，下面是代码：

 

  1. /*
  2.         时间：2018年11月24日
  3.         功能：将input里的数据进行快速傅里叶变换  并且输出
  4. */
  6. #include <stdio.h>
  7. #include <math.h>
  8. #define PI 3.1415926
  9. #define FFT_LENGTH                8
  10. double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};
  11. struct complex1{                //定义一个复数结构体
  12.         double real;        //实部
  13.         double image;        //虚部
  14. };
  15. //将input的实数结果存放为复数
  16. struct complex1 result_dat[8];
  17. /*
  18.         虚数的乘法
  19. */
  20. struct complex1 con_complex(struct complex1 a,struct complex1 b){
  21.         struct complex1 temp;
  22.         temp.real=(a.real*b.real)-(a.image*b.image);
  23.         temp.image=(a.image*b.real)+(a.real*b.image);
  24.         return temp;
  25. }
  27. /*
  28.         简单的a的b次方
  29. */
  30. int mypow(int a,int b){
  31.         int i,sum=a;
  32.         if(b==0)return 1;
  33.         for(i=1;i<b;i++){
  34.                 sum*=a;
  35.         }
  36.         return sum;
  37. }
  38. /*
  39.         简单的求以2为底的正整数
  40. */
  41. int log2(int n){
  42.         unsigned i=1;
  43.         int sum=1;
  44.         for(i;;i++){
  45.                 sum*=2;
  46.                 if(sum>=n)break;
  47.         }
  48.         return i;
  49. }
  50. /*
  51.         简单的交换数据的函数
  52. */
  53. void swap(struct complex1 *a,struct complex1 *b){
  54.         struct complex1 temp;
  55.     temp=*a;
  56.     *a=*b;
  57.     *b=temp;
  58. }
  59. /*
  60.         dat为输入数据的数组
  61.         N为抽样次数  也代表周期  必须是2^N次方
  62. */
  63. void fft(struct complex1 dat[],unsigned char N){
  64.         /*最终  dat_buf计算出 当前蝶形运算奇数项与W  乘积
  65.                         dat_org存放上一个偶数项的值
  66.         */
  67.         struct complex1 dat_buf,dat_org;
  68.         /*        L为几级蝶形运算    也代表了2进制的位数
  69.                 n为当前级蝶形的需要次数  n最初为N/2 每级蝶形运算后都要/2
  70.                 i j为倒位时要用到的自增符号  同时  i也用到了L碟级数   j是计算当前碟级的计算次数
  71.                 re_i i_copy均是倒位时用到的变量
  72.                 k为当前碟级  cos(2*pi/N*k)的  k   也是e^(-j2*pi/N)*k  的  k
  73.         */
  74.         unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;
  75.         //经过观察，发现每级蝶形运算需要N/2次运算，共运算N/2*log2N  次
  76.         unsigned char fft_counter=0;
  77.         //在此要进行补2   N必须是2^n   在此略
  78.         //蝶形级数  (L级)
  79.         L=log2(N);
  80.         //计算每级蝶形计算的次数(这里只是一个初始值)  之后每次要/2
  81.         //n=N/2;
  83.         //对dat的顺序进行倒位
  84.         for(i=1;i<N/2;i++){
  85.                 i_copy=i;
  86.                 re_i=0;
  87.                 for(j=L-1;j>0;j--){
  88.                         //判断i的副本最低位的数字  并且移动到最高位  次高位  ..
  89.                         //re_i为交换的数   每次它的数字是不能移动的 并且循环之后要清0
  90.                         re_i|=((i_copy&0x01)<<j);
  91.                         i_copy>>=1;
  92.                 }
  93.                 swap(&dat[i],&dat[re_i]);
  94.         }
  95.         //进行fft计算
  96.         for(i=0;i<L;i++){
  98.                 fft_flag=1;
  99.                 fft_counter=0;
  100.                 for(j=0;j<N;j++){
  101.                         if(fft_counter==mypow(2,i)){                //控制隔几次，运算几次,
  102.                                 fft_flag=0;
  103.                         }else if(fft_counter==0){                //休止结束，继续运算
  104.                                 fft_flag=1;
  105.                         }
  106.                         //当不判断这个语句的时候  fft_flag保持  这样就可以持续运算了
  107.                         if(fft_flag){
  108.                                 dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));
  109.                                 dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));
  110.                                 dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);
  111.                                                                                                 //计算 当前蝶形运算奇数项与W  乘积
  113.                                 dat_org.real=dat[j].real;
  114.                                 dat_org.image=dat[j].image;                //暂存
  116.                                 dat[j].real=dat_org.real+dat_buf.real;
  117.                                 dat[j].image=dat_org.image+dat_buf.image;
  118.                                                                                                         //实部加实部   虚部加虚部
  120.                                 dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;
  121.                                 dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;
  122.                                                                                                         //实部减实部        虚部减虚部
  124.                                 k++;
  125.                                 fft_counter++;
  126.                         }else{
  127.                                 fft_counter--;                                //运算几次，就休止几次
  128.                                 k=0;
  129.                         }
  130.                 }
  131.         }
  134. }
  135. void main(){
  136.         int i;
  137.         //先将输入信号转换成复数
  138.         for(i=0;i<FFT_LENGTH;i++){
  139.                 result_dat[i].image=0;
  140.                                                                 //输入信号是二维的，暂时不存在复数
  141.                 result_dat[i].real=input[i];
  142.                 //result_dat[i].real=10;
  143.                                                                 //输入信号都为实数
  144.         }
  145.         fft(result_dat,FFT_LENGTH);
  146.         for(i=0;i<FFT_LENGTH;i++){
  147.                 input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);
  148.                 //取模
  149.                 printf("%lf\n",input[i]);
  150.         }
  151.         while(1);
  152. }

复制代码

这个程序中input这个数组是输入信号，在这里只模拟抽样了8次，输出的数据也是input，如果想看其它序列的话，可以改变FFT_LENGTH的值以及 input里的内容，程序输出的是实部和虚部的模，如果单纯想看实部或者虚部的幅度的话，请自行修改程序~

这就是本次浅谈FFT以及FFT算法的基本实现的全部内容了

原创内容转载其注明原作者.
参考书籍：数字信号处理

来源： 通俗细致谈FFT以及FFT算法的基本实现_通信与网络-面包板社区

历史上的今天
1月
19
weinxin
我的微信
一个码农、工程狮、集能量和智慧于一身的、DIY高手、小伙伴er很多的、80后奶爸。
Igor
  • 版权声明： 发表于 2019-01-1911:46:29
  • 转载注明：http://blog.tsingmac.com/prolions/software/manual/1050/
图象处理基本算法[整理] 图形图像

图象处理基本算法[整理]

1)将256*256分辨率的图像变为128*128分辨率可以将源图像划分成2*2的子图像块,然后将2*2的 子图像块的所有像素颜色均按照F(i,j)的颜色值进行设定，达到降低分辨率的目的。 如: F(...
03月02日 6,160 views 评论
匿名

发表评论

匿名网友 填写信息

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: