devbox
笔触狂放9
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

快速傅里叶变换(FFT)c语言实现_fft代码

作者：笔触狂放9 | 2024-02-18 05:03:35

赞

踩

fft代码

基本原理在这里就不多讲了，可以看看其他高浏览量的博文，这篇文章针对c语言的实现

复数运算算子

        我们都知道C语言本身是没有复数运算的，很多DSP、单片机要用到也没有开源库可以使用复数运算，针对FFT在硬件上运行只能手动从底层开始

定义复数类型

        这里用最简单高效的方法——结构体

  1. struct complex
  2. {
  3.     double real;
  4.     double image;
  5. };

复数加法

  1. struct complex complex_add(struct complex c1,struct complex c2)     //复数加法
  2. {
  3.     struct complex p;
  4.  
  5.     p.real = c1.real + c2.real;
  6.     p.image = c1.image + c2.image;
  7.  
  8.     return p;
  9.  
  10. }

复数减法

  1. struct complex complex_sub(struct complex c1,struct complex c2)     //复数减
  2. {
  3.     struct complex p;
  4.  
  5.     p.real = c1.real - c2.real;
  6.     p.image = c1.image - c2.image;
  7.  
  8.     return p;
  9.  
  10. }

复数乘法

  1. struct complex complex_multi(struct complex c1,struct complex c2)  //复数乘法
  2. {
  3.     struct complex c3;
  4.     c3.real=c1.real*c2.real - c1.image *c2.image;
  5.     c3.image = c2.real* c1.image + c1.real*c2.image;
  6.  
  7.     return c3;
  8. };

复数取模

  1. double mold_length(struct complex c)        //幅度
  2. {
  3.     return sqrt(c.real * c.real + c.image * c.image);
  4. };

旋转因子

    教科书上的旋转因子一般长成这样W_{N}^{nk}把它变成可以运算的形式   

W_{N}^{nk} = e^{-j\frac{2\pi}{N}nk}=cos(\frac{2\pi}{N}nk)-jsin(\frac{2\pi}{N}nk)

  1. struct complex rotation_factor(int N,int n,int k)       //旋转因子
  2. {
  3.     struct complex w;
  4.  
  5.     w.real = cos(2* PI * n * k /N);
  6.     w.image = - sin(2* PI * n * k /N);        //这里的PI是宏定义
  7.  
  8.     return w;
  9. }

反位序操作

        我用的方法是基于时间抽取的基2FFT （DIT-FFT）所以开始进入FFT之前要先把序列顺序取反位序

  1. int reverse_num(int l,int oringin_num)          //反位序
  2. {
  3.     int q=0,m=0;
  4.  
  5.     for(int k=l-1;k>=0;k--)
  6.     {
  7.         q = oringin_num % 2;
  8.         m += q*(1<<k);
  9.         oringin_num = oringin_num / 2;
  10.     }
  11.  
  12.     return m;
  13. }

解释一下这个代码含义，用取余操作取出最低位，把最低位放在m中最高位置，取出后除2操作把每一位右移，直到为0

FFT

我们首先拿出祖传的蝶形运算信号流图，

 在其中抽取一个蝶形详细分析，注意这其中上半部分是相加【+】下半部分是相减【-】，而且下半部分都会带有一个旋转因子W_{N}^{r}，先记住这个特征，接下来编程会利用这个特征

 为了方便编程，把它表达成数学式子，这其中前一级与后一级的关系用等式表达

X_{m}(k)=X_{m-1}(k)+X_{m-1}(j)W_{N}^{r}

X_{m}(j)=X_{m-1}(k)-X_{m-1}(j)W_{N}^{r}

现在问题就是k和j的确定，重新观察完整的蝶形图在第【1】级两两运算的节点序号是【0、1】【2、3】【4、5】【6、7】；在第【2】级两两运算的节点序号是【0、2】【1、3】【4、6】【5、7】；第【3】级运算的节点是【0、4】【1、5】【2、6】【3、7】。注意这里不要用图上的反位置序，而是反位后的顺序，也就是下图红色序号

首先【j】与【k】的间隔可以容易观察出规律，每增加一级，间隔就会扩大两倍，换句话说【j】与【k】的间搁可以表达为下面的式子

dist=j-k=2^{m-1} (m=1,2,3,...)

好了，这里我们破解了第一条规律，急需确定的是这个旋转因子W_{N}^{r}当中的r，这里我没有观察到规律，只好直接引用教材上的解释

 我们可以简单验证一下：

比如第【3】级 【k=2】的情况

r = 2\times2^{3-3}=2   W_{N}^{r} = W_{8}^{2}

第【2】级【k=5】的情况

r = 5\times2^{3-2} = 10   W_{N}^{r} = W_{8}^{10}=W_{8}^{8}W_{8}^{2}=W_{8}^{2}

两次验证都符合实际FFT蝶形图的结果。

接下来看似没有急需破解的规律了，但是真正编程时候还是存在一个问题，【k】的取值不是从上到下按顺序的例如第一级的【k】取值为【0，2，4，6】第二级的取值为【0，1，4，5】第三级的取值为【0，1，2，3】。

这样如果按找从头到尾的顺序确实找不到什么规律，但是整体没规律，部分总会有规律吧，注意看我笔迹画的线为啥颜色是这样分的呢，没有发现很有规律吗？我们按找相同颜色分组，第一级的第一组两个蝶形运算的【k】之间相差距离刚好是【dist】的两倍,第二级的第一组两个蝶形运算的【k】之间相差的距离也是【dist】的两倍，但是到了第三级， 【dist】的两倍是【8】如果【k=0】 那么【k+8】就会超出最大序号【7】，那么考虑第三级分成【4】组，每一组只有一个蝶形运算，有别于其他两个，可以用for循环这样写

for(k=0;k<len;k+=(dist<<1))    //在二进制左移一位相当于乘2 ，这里len是FFT的点数N

同时考虑每一级内又有很多组，可以用for循环嵌套

  1. for(j=0;j<dist;j++)
  2.     for(k=0;k<len;k+=(dist<<1))

最终把分析的结果整合到一起，写成单独的函数 ，注意一下，程序中的一些标号和前面分析的会有点不同

  1. void fft(int len, struct complex in_x[],struct complex out_y[])
  2. {
  3.  
  4.     /*
  5.         param len 序列长度，目前只能是2的指数
  6.         param in_x输入的序列
  7.         param out_y输出的序列
  9.     */
  10.  
  11.  
  12.     int l,k,r,z,dist,q,j;       //l是级
  13.     struct complex w,tmp;
  14.     struct complex in_x_mem[len];
  15.  
  16.     l = log2(len);
  17.  
  18.     for(k=0;k<len;k++)
  19.     {
  20.         in_x_mem[k] = in_x[reverse_num(l,k)];       //反位序号操作
  21.     }
  22.  
  23.     for(r = 0;r<l;r++)      //遍历每一级蝶形运算
  24.     {
  25.  
  26.         dist = 1<<r;        //提前计算每一级的间隔距离
  27.  
  28.         for( j=0;j<dist;j++ )      //计算策略是拆成上下两组，先上计算，后下计算，j是计算的起始序号
  29.         {
  30.  
  31.             for(k=j;k<len;k+=(dist<<1))     //不好解释，得看画图理解
  32.             {
  33.                 q = k+dist; //q同一组蝶形运算第二个序号
  34.                 z = k << (l - r -1);    //确定旋转因子的指数
  35.  
  36.  
  37.                 w = rotation_factor(len,1,z);
  38.                 //由于不是并行计算，必须先缓存
  39.                 tmp = in_x_mem[k];
  40.  
  41.                 in_x_mem[k] = complex_add( in_x_mem[k] , complex_multi(w, in_x_mem[q]) );
  42.                 in_x_mem[q] = complex_sub(tmp , complex_multi(w, in_x_mem[q]) );
  43.  
  44.                 
  45.  
  46.             }
  47.  
  48.  
  49.         }
  50.     }
  51.     
  52.     memcpy(out_y,in_x_mem,len*sizeof(struct complex));    //拷贝到输出的数组
  53. }

结果测试

这里我用序列x[n]=cos(\pi n)=(-1)^n输入到FFT，输出应该是X[k]=8\delta[k-4]

 对照结果，和计算一样，程序成功

完整程序

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5.  
  6.  
  7. #define PI 3.1415926
  8. struct complex
  9. {
  10.     double real;
  11.     double image;
  12. };
  13. struct complex complex_add(struct complex c1,struct complex c2);
  14. struct complex complex_sub(struct complex c1,struct complex c2);
  15. struct complex complex_multi(struct complex c1,struct complex c2);
  16. struct complex rotation_factor(int N,int n,int k);
  17. double mold_length(struct complex c);
  18. void fft(int len, struct complex in_x[],struct complex out_y[]);
  19.  
  20.  
  21.  
  22.  
  23. int main()
  24. {
  25.     int sam[8] = {1,-1,1,-1,1,-1,1,-1};
  26.     int jhg[8] = {0,0,0,0,0,0,0,0};
  27.     struct complex x[8];
  28.     struct complex y[8];
  29.  
  30.     for(int i=0;i<8;i++)
  31.     {
  32.         x[i].real = sam[i];
  33.         x[i].image = jhg[i];
  34.     }
  35.     printf("时域序列\n");
  36.  
  37.     for(int i=0;i<8;i++)
  38.     {
  39.         printf("(%.2f, %.2fi) \n",x[i].real,x[i].image);
  40.     }
  41.  
  42.     fft(8,x,y);
  43.  
  44.     printf("频域序列\n");
  45.  
  46.     for(int i=0;i<8;i++)
  47.     {
  48.         printf("(%.2f, %.2fi)\n",y[i].real,y[i].image);
  49.     }
  50.  
  51.     
  52.  
  53.     return 0;
  54. }
  55.  
  56. struct complex complex_add(struct complex c1,struct complex c2)     //复数加法
  57. {
  58.     struct complex p;
  59.  
  60.     p.real = c1.real + c2.real;
  61.     p.image = c1.image + c2.image;
  62.  
  63.     return p;
  64.  
  65. }
  66.  
  67. struct complex complex_sub(struct complex c1,struct complex c2)     //复数减
  68. {
  69.     struct complex p;
  70.  
  71.     p.real = c1.real - c2.real;
  72.     p.image = c1.image - c2.image;
  73.  
  74.     return p;
  75.  
  76. }
  77.  
  78. struct complex complex_multi(struct complex c1,struct complex c2)  //复数乘法
  79. {
  80.     struct complex c3;
  81.     c3.real=c1.real*c2.real - c1.image *c2.image;
  82.     c3.image = c2.real* c1.image + c1.real*c2.image;
  83.  
  84.     return c3;
  85. };
  86.  
  87.  
  88.  
  89. struct complex rotation_factor(int N,int n,int k)       //旋转因子
  90. {
  91.     struct complex w;
  92.  
  93.     w.real = cos(2* PI * n * k /N);
  94.     w.image = - sin(2* PI * n * k /N);
  95.  
  96.     return w;
  97. }
  98.  
  99. double mold_length(struct complex c)        //幅度
  100. {
  101.     return sqrt(c.real * c.real + c.image * c.image);
  102. };
  103.  
  104.  
  105. int reverse_num(int l,int oringin_num)          //反位序
  106. {
  107.     int q=0,m=0;
  108.  
  109.     for(int k=l-1;k>=0;k--)
  110.     {
  111.         q = oringin_num % 2;
  112.         m += q*(1<<k);
  113.         oringin_num = oringin_num / 2;
  114.     }
  115.  
  116.     return m;
  117. }
  118.  
  119.  
  120. void fft(int len, struct complex in_x[],struct complex out_y[])
  121. {
  122.     /*
  123.         param len 序列长度，目前只能是2的指数
  124.         param in_x输入的序列
  125.         param out_y输出的序列
  127.     */
  128.  
  129.     int l,k,r,z,dist,q,j;       //l是级
  130.     struct complex w,tmp;
  131.     struct complex in_x_mem[len];
  132.  
  133.     l = log2(len);
  134.  
  135.     for(k=0;k<len;k++)
  136.     {
  137.         in_x_mem[k] = in_x[reverse_num(l,k)];       //反位序号操作
  138.     }
  139.  
  140.     for(r = 0;r<l;r++)      //遍历每一级蝶形运算
  141.     {
  142.  
  143.         dist = 1<<r;        //提前计算每一级的间隔距离
  144.  
  145.         for( j=0;j<dist;j++ )      //计算策略是拆成上下两组，先上计算，后下计算，j是计算的起始序号
  146.         {
  147.             for(k=j;k<len;k+=(dist<<1))     //不好解释，得画图理解
  148.             {
  149.                 q = k+dist; //q同一组蝶形运算第二个序号
  150.                 z = k << (l - r -1);    //确定旋转因子的指数
  151.  
  152.                 w = rotation_factor(len,1,z);
  153.                 //由于不是并行计算，必须先缓存
  154.                 tmp = in_x_mem[k];
  155.  
  156.                 in_x_mem[k] = complex_add( in_x_mem[k] , complex_multi(w, in_x_mem[q]) );
  157.                 in_x_mem[q] = complex_sub(tmp , complex_multi(w, in_x_mem[q]) );
  158.  
  159.             }
  160.         }
  161.     }
  162.     memcpy(out_y,in_x_mem,len*sizeof(struct complex));
  163. }
  164.  
  165.  
  166.
本文内容由网友自发贡献，转载请注明出处：https://www.wpsshop.cn/w/笔触狂放9/article/detail/104994?site
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号