返回列表 发新帖本帖赏金 10.00元(功能说明)

浅谈FFT以及FFT算法的基本实现

[复制链接]
2236|21
 楼主 | 2018-12-14 22:34 | 显示全部楼层 |阅读模式
本帖最后由 enderman1 于 2018-12-20 13:05 编辑

相信网上现在有很多关于FFT的教程,我曾经也参阅了很多网上的教程,感觉都不怎么通俗易懂。在基本上的研究FFT,并且通过编程的形式实现之后。我决定写一篇通俗易懂的关于FFT的讲解。因此我在接下来的叙述中尽量非常通俗细致的讲解。
本人最早知道傅里叶变换的时候是沉迷于音乐的频谱跳动无法自拔,当时就很想做一个音乐频谱显示器。搜阅了很多资料之后,才了解到傅里叶变换,和FFT。当然这都是以前的事情了,经过了系统的学习+2个星期的研究,自制了一个FFT的算法,不敢说速度上的优势,但是个人认为是一种通俗易懂的实现方法。经过实际的VC++模拟实验、和STM32跑的也很成功。

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

         ,k=0,1,…N-1

           n=0,1…N-1



      其中:



但是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/2DFT只需(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. */

  5. #include <stdio.h>
  6. #include <math.h>
  7. #define PI 3.1415926
  8. #define FFT_LENGTH                8
  9. double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};
  10. struct complex1{                //定义一个复数结构体
  11.         double real;        //实部
  12.         double image;        //虚部
  13. };        
  14. //将input的实数结果存放为复数
  15. struct complex1 result_dat[8];
  16. /*
  17.         虚数的乘法
  18. */
  19. struct complex1 con_complex(struct complex1 a,struct complex1 b){
  20.         struct complex1 temp;
  21.         temp.real=(a.real*b.real)-(a.image*b.image);
  22.         temp.image=(a.image*b.real)+(a.real*b.image);
  23.         return temp;
  24. }

  25. /*
  26.         简单的a的b次方
  27. */
  28. int mypow(int a,int b){
  29.         int i,sum=a;
  30.         if(b==0)return 1;
  31.         for(i=1;i<b;i++){
  32.                 sum*=a;        
  33.         }
  34.         return sum;
  35. }
  36. /*
  37.         简单的求以2为底的正整数
  38. */
  39. int log2(int n){
  40.         unsigned i=1;
  41.         int sum=1;
  42.         for(i;;i++){
  43.                 sum*=2;
  44.                 if(sum>=n)break;
  45.         }
  46.         return i;
  47. }
  48. /*
  49.         简单的交换数据的函数
  50. */
  51. void swap(struct complex1 *a,struct complex1 *b){
  52.         struct complex1 temp;
  53.     temp=*a;
  54.     *a=*b;
  55.     *b=temp;
  56. }
  57. /*
  58.         dat为输入数据的数组
  59.         N为抽样次数  也代表周期  必须是2^N次方
  60. */
  61. void fft(struct complex1 dat[],unsigned char N){        
  62.         /*最终  dat_buf计算出 当前蝶形运算奇数项与W  乘积
  63.                         dat_org存放上一个偶数项的值
  64.         */
  65.         struct complex1 dat_buf,dat_org;
  66.         /*        L为几级蝶形运算    也代表了2进制的位数
  67.                 n为当前级蝶形的需要次数  n最初为N/2 每级蝶形运算后都要/2
  68.                 i j为倒位时要用到的自增符号  同时  i也用到了L碟级数   j是计算当前碟级的计算次数
  69.                 re_i i_copy均是倒位时用到的变量
  70.                 k为当前碟级  cos(2*pi/N*k)的  k   也是e^(-j2*pi/N)*k  的  k
  71.         */
  72.         unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;
  73.         //经过观察,发现每级蝶形运算需要N/2次运算,共运算N/2*log2N  次  
  74.         unsigned char fft_counter=0;
  75.         //在此要进行补2   N必须是2^n   在此略
  76.         //蝶形级数  (L级)
  77.         L=log2(N);        
  78.         //计算每级蝶形计算的次数(这里只是一个初始值)  之后每次要/2
  79.         //n=N/2;

  80.         //对dat的顺序进行倒位
  81.         for(i=1;i<N/2;i++){
  82.                 i_copy=i;
  83.                 re_i=0;
  84.                 for(j=L-1;j>0;j--){
  85.                         //判断i的副本最低位的数字  并且移动到最高位  次高位  ..
  86.                         //re_i为交换的数   每次它的数字是不能移动的 并且循环之后要清0
  87.                         re_i|=((i_copy&0x01)<<j);               
  88.                         i_copy>>=1;
  89.                 }
  90.                 swap(&dat[i],&dat[re_i]);
  91.         }
  92.         //进行fft计算
  93.         for(i=0;i<L;i++){
  94.                
  95.                 fft_flag=1;
  96.                 fft_counter=0;
  97.                 for(j=0;j<N;j++){
  98.                         if(fft_counter==mypow(2,i)){                //控制隔几次,运算几次,
  99.                                 fft_flag=0;
  100.                         }else if(fft_counter==0){                //休止结束,继续运算
  101.                                 fft_flag=1;
  102.                         }
  103.                         //当不判断这个语句的时候  fft_flag保持  这样就可以持续运算了
  104.                         if(fft_flag){
  105.                                 dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));
  106.                                 dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));
  107.                                 dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);
  108.                                                                                                 //计算 当前蝶形运算奇数项与W  乘积

  109.                                 dat_org.real=dat[j].real;
  110.                                 dat_org.image=dat[j].image;                //暂存

  111.                                 dat[j].real=dat_org.real+dat_buf.real;
  112.                                 dat[j].image=dat_org.image+dat_buf.image;               
  113.                                                                                                         //实部加实部   虚部加虚部

  114.                                 dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;
  115.                                 dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;
  116.                                                                                                         //实部减实部        虚部减虚部

  117.                                 k++;
  118.                                 fft_counter++;
  119.                         }else{
  120.                                 fft_counter--;                                //运算几次,就休止几次
  121.                                 k=0;
  122.                         }
  123.                 }
  124.         }

  125.         
  126. }
  127. void main(){
  128.         int i;
  129.         //先将输入信号转换成复数
  130.         for(i=0;i<FFT_LENGTH;i++){
  131.                 result_dat[i].image=0;        
  132.                                                                 //输入信号是二维的,暂时不存在复数
  133.                 result_dat[i].real=input[i];
  134.                 //result_dat[i].real=10;
  135.                                                                 //输入信号都为实数
  136.         }
  137.         fft(result_dat,FFT_LENGTH);
  138.         for(i=0;i<FFT_LENGTH;i++){
  139.                 input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);
  140.                 //取模
  141.                 printf("%lf\n",input[i]);
  142.         }
  143.         while(1);
  144. }
复制代码

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

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


原创内容转载其注明原作者.

参考书籍:数字信号处理

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

打赏榜单

21ic小喇叭 打赏了 10.00 元 2018-12-28
理由:欢迎继续在21ic进行分享

| 2018-12-17 18:55 | 显示全部楼层
这个格式调整一下更好

评论

enderman1 2018-12-20 13:08 回复TA
谢谢提醒:-) 
| 2018-12-18 08:53 | 显示全部楼层
收藏收藏,FFT,确实很有用,感谢楼主给我们普及相关知识。
| 2018-12-18 14:29 | 显示全部楼层
mark
| 2018-12-18 16:17 | 显示全部楼层
mark
| 2018-12-19 09:45 | 显示全部楼层
把 file:///D:/temp_data 什么的重新上传一下。

评论

enderman1 2018-12-20 12:46 回复TA
非常抱歉,其实那个是公式,但是不知道为什么无法上传了,我尽快补上 
| 2018-12-19 18:45 | 显示全部楼层
烧脑,但还是给楼主点赞
| 2018-12-19 21:18 | 显示全部楼层
学习一下.
 楼主 | 2018-12-20 13:08 | 显示全部楼层
非常感谢大家的支持,由于楼主是第一次在这发文章,难免会有很多错误,请谅解;还有就是文章里的什么什么乱码的错误已经改正过来了,以后可能会发表跟多的文章  以供参考学习。谢谢~
| 2018-12-20 15:31 | 显示全部楼层
enderman1 发表于 2018-12-20 13:08
非常感谢大家的支持,由于楼主是第一次在这发文章,难免会有很多错误,请谅解;还有就是文章里的什么什么乱 ...

加油!
| 2018-12-25 11:25 | 显示全部楼层
你这程序也是和网上一样的,那3层叠形运算,不用现成的处理方法 ,自己写很烧脑子,还是抄网上的好
 楼主 | 2019-2-2 15:54 | 显示全部楼层
今天在做数字信号处理中滤波这一方面的工作的时候,突然想起了当时写的算法还有一些瑕疵。因为傅里叶变换不能只变换过去而不变回来;因此我更新了一下算法,其中ifft这个变量相当于一个标志位吧,当它等于1时是fft的逆变换,当它是0的时候,还是fft变换。这样补充之后算是整个傅里叶变换对了,因此根据公式我对之前的做了一些小小的修改,也算是对我过去工作的一个完善吧~
以下是源代码:
  1. /*
  2.         dat为输入数据的数组
  3.         N为抽样次数  也代表周期  必须是2^N次方
  4.         ifft 为0 时启动fft  为1时 fft逆变换
  5. */
  6. void fft(struct complex1 dat[],unsigned char N,unsigned char ifft){        
  7.         /*最终  dat_buf计算出 当前蝶形运算奇数项与W  乘积
  8.                         dat_org存放上一个偶数项的值
  9.         */
  10.         struct complex1 dat_buf,dat_org;
  11.         /*        L为几级蝶形运算    也代表了2进制的位数
  12.                 n为当前级蝶形的需要次数  n最初为N/2 每级蝶形运算后都要/2
  13.                 i j为倒位时要用到的自增符号  同时  i也用到了L碟级数   j是计算当前碟级的计算次数
  14.                 re_i i_copy均是倒位时用到的变量
  15.                 k为当前碟级  cos(2*pi/N*k)的  k   也是e^(-j2*pi/N)*k  的  k
  16.         */
  17.         unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;
  18.         //经过观察,发现每级蝶形运算需要N/2次运算,共运算N/2*log2N  次  
  19.         unsigned char fft_counter=0;
  20.         //在此要进行补2   N必须是2^n   在此略

  21.         //double iw_real;
  22.         //蝶形级数  (L级)
  23.         L=log2(N);        
  24.         //计算每级蝶形计算的次数(这里只是一个初始值)  之后每次要/2
  25.         //n=N/2;

  26.         //对dat的顺序进行倒位
  27.         for(i=1;i<N/2;i++){
  28.                 i_copy=i;
  29.                 re_i=0;
  30.                 for(j=L-1;j>0;j--){
  31.                         //判断i的副本最低位的数字  并且移动到最高位  次高位  ..
  32.                         //re_i为交换的数   每次它的数字是不能移动的 并且循环之后要清0
  33.                         re_i|=((i_copy&0x01)<<j);               
  34.                         i_copy>>=1;
  35.                 }
  36.                 swap(&dat[i],&dat[re_i]);
  37.         }
  38.         //进行fft计算
  39.         for(i=0;i<L;i++){
  40.                
  41.                 fft_flag=1;
  42.                 fft_counter=0;
  43.                 for(j=0;j<N;j++){
  44.                         if(fft_counter==mypow(2,i)){                //控制隔几次,运算几次,
  45.                                 fft_flag=0;
  46.                         }else if(fft_counter==0){                //休止结束,继续运算
  47.                                 fft_flag=1;
  48.                         }
  49.                         //当不判断这个语句的时候  fft_flag保持  这样就可以持续运算了
  50.                         if(fft_flag){
  51.                                 /*
  52.                                 dat_buf=dat[j];
  53.                                 dat[j]=dat[j]+dat[j+mypow(2,i)]*cos((2*PI*k)/(N/mypow(2,L-i-1)));
  54.                                 dat[j+mypow(2,i)]=dat_buf-dat[j+mypow(2,i)]*cos((2*PI*k)/(N/mypow(2,L-i-1)));        
  55.                                 */
  56.                                 dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));
  57.                                 if(ifft){
  58.                                         dat_buf.image=sin((2*PI*k)/(N/mypow(2,L-i-1)));
  59.                                 }else{
  60.                                         dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));
  61.                                 }

  62.                                 dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);
  63.                                                                                                 //计算 当前蝶形运算奇数项与W  乘积

  64.                                 dat_org.real=dat[j].real;
  65.                                 dat_org.image=dat[j].image;                //暂存

  66.                                 dat[j].real=dat_org.real+dat_buf.real;
  67.                                 dat[j].image=dat_org.image+dat_buf.image;               
  68.                                                                                                         //实部加实部   虚部加虚部

  69.                                 dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;
  70.                                 dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;
  71.                                                                                                         //实部减实部        虚部减虚部

  72.                                 k++;
  73.                                 fft_counter++;
  74.                         }else{
  75.                                 fft_counter--;                                //运算几次,就休止几次
  76.                                 k=0;
  77.                         }
  78.                 }
  79.         }
  80.         if(ifft){
  81.                 for(i=0;i<N;i++){
  82.                         dat[i].image/=N;
  83.                         dat[i].real/=N;
  84.                 }
  85.         }
  86.         
  87. }
复制代码


评论

enderman1 2019-2-2 15:57 回复TA
对于这个算法的准确性本人还是亲自和matlab进行对比试验过的,准确度还是很不错的,另外输入数组支持复数。请大家参考吧~ 
| 2019-2-3 09:21 | 显示全部楼层
从DFT到FFT是自然而然的。从离散周期信号的傅立叶级数到离散傅立叶变换,不同的教科书有不同的过渡方法,令人信服的推导方法不多见。
| 2019-2-3 09:21 | 显示全部楼层
从DFT到FFT是自然而然的。从离散周期信号的傅立叶级数到离散傅立叶变换,不同的教科书有不同的过渡方法,令人信服的推导方法不多见。
| 2019-2-3 09:31 | 显示全部楼层
采样频率、频谱分辨率,截断效应,比纯粹的算法难理解。
| 2019-2-3 15:00 | 显示全部楼层
干货,收存
| 2019-2-5 16:09 | 显示全部楼层
收藏
| 2019-2-9 14:22 | 显示全部楼层
纯粹的算法随便可以找得到,傅里叶变换思想的理解才是最重要的。

评论

enderman1 2019-2-9 22:24 回复TA
是的 
扫描二维码,随时随地手机跟帖
返回列表 发新帖 本帖赏金 10.00元(功能说明)
*滑动验证:
您需要登录后才可以回帖 登录 | 注册

本版积分规则

我要发帖 投诉建议 创建版块 申请版主

快速回复

您需要登录后才可以回帖
登录 | 注册
高级模式

论坛热帖

关闭

热门推荐上一条 /4 下一条

快速回复 返回顶部 返回列表