[开发工具] C语言实现FFT算法

[复制链接]
 楼主| aspoke 发表于 2023-2-23 09:57 | 显示全部楼层 |阅读模式
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4. #define N 8

  5. /*复数结构体类型*/
  6. typedef struct{
  7. double realPart;
  8. double imaginaryPart;
  9. }complexNumber;

  10. complexNumber x[N], W[N/2]; //序列数组以及旋转因子数组
  11. int size_x=N;
  12. double PI=atan(1)*4;        //圆周率
  13. void fastFourierOperation();//快速傅里叶变换算法
  14. void initRotationFactor();  //生成旋转因子数组
  15. void changePosition();      //偶奇变换算法
  16. void outputArray();         //遍历输出数组
  17. void add(complexNumber ,complexNumber ,complexNumber*); //复数加法
  18. void mul(complexNumber ,complexNumber ,complexNumber*); //复数乘法
  19. void sub(complexNumber ,complexNumber ,complexNumber*); //复数减法

  20. /*长度为1024的冲激序列 */
  21. void init1()               
  22. {
  23.         printf("Example:度为1024的冲激序列\n");
  24.         x[0].realPart=1;
  25.         x[0].imaginaryPart=0;
  26.         for(int i=1;i<N;i++)
  27.         {
  28.                 x[i].realPart=0;
  29.                 x[i].imaginaryPart=0;
  30.         }
  31. }

  32. /*长度为1024的矩形序列 */
  33. void init2()               
  34. {               
  35.         printf("Example:长度为1024的矩形序列\n");
  36.         for(int i=0;i<N;i++)
  37.         {
  38.                 x[i].realPart=1;
  39.                 x[i].imaginaryPart=0;
  40.         }
  41. }

  42. /*长度为1024的sin[2πk/1024]序列 */
  43. void init3()               
  44. {            
  45.         printf("Example:长度为1024的sin[2πk/1024]序列\n");
  46.         for(int i=0;i<N;i++)
  47.         {
  48.                 x[i].realPart=sin(2*PI*i/N);
  49.                 x[i].imaginaryPart=0;
  50.         }
  51. }

  52. /*长度为1024的cos[2πk/1024]序列 */
  53. void init4()              
  54. {            
  55.         printf("Example:长度为1024的cos[2πk/1024]序列\n");
  56.         for(int i=0;i<N;i++)
  57.         {
  58.                 x[i].realPart=cos(2*PI*i/N);
  59.                 x[i].imaginaryPart=0;
  60.         }
  61. }

  62. int main()
  63. {
  64.     init3();               //调用不同的初始化函数,使用不同的例程进行测试
  65.         printf("=============================================================================================================================================================================================\n");
  66.         printf("输出原序列:\n");
  67.         outputArray();
  68.         initRotationFactor();  //生成旋转因子数组
  69.         fastFourierOperation();//调用快速傅里叶变换
  70.         printf("\n\n=============================================================================================================================================================================================\n");
  71.         printf("输出FFT后的结果:\n");
  72.         outputArray();         //遍历输出数组
  73.         return 0;
  74. }


  75. /*快速傅里叶变换*/
  76. void fastFourierOperation()
  77. {
  78.         int l=0;
  79.         complexNumber up,down,product;
  80.         changePosition();                            //偶奇变换算法
  81.         for(int i=0;i<log(size_x)/log(2);i++)        //一级蝶形运算
  82.         {   
  83.                 l=1<<i;
  84.                 for(int j=0;j<size_x;j+= 2*l )          //一组蝶形运算
  85.                 {            
  86.                         for(int k=0;k<l;k++)                //一个蝶形运算
  87.                         {      
  88.                                 mul(x[j+k+l],W[size_x*k/2/l],&product);
  89.                                 add(x[j+k],product,&up);
  90.                                 sub(x[j+k],product,&down);
  91.                                 x[j+k]=up;
  92.                                 x[j+k+l]=down;
  93.                         }
  94.                 }
  95.         }
  96. }


  97. /*生成旋转因子数组*/
  98. void initRotationFactor()
  99. {
  100.         for(int i=0;i<size_x/2;i++)
  101.         {
  102.                 W[i].realPart=cos(2*PI/size_x*i);//用欧拉公式计算旋转因子
  103.                 W[i].imaginaryPart=-1*sin(2*PI/size_x*i);
  104.         }
  105. }


  106. /*偶奇变换算法*/
  107. void changePosition()      
  108. {
  109.         int j=0,k;//待会儿进行运算时需要用到的变量
  110.         complexNumber temp;
  111.         for (int i=0;i<size_x-1;i++)
  112.         {
  113.                 if(i<j)
  114.                 {
  115.                         //若倒序数大于顺序数,进行变址(换位);否则不变,防止重复交换,变回原数
  116.                         temp=x[i];
  117.                         x[i]=x[j];
  118.                         x[j]=temp;
  119.                 }
  120.                 k=size_x/2;    //判断j的最高位是否为0
  121.                 while(j>=k)
  122.                 {              //最高位为1
  123.                         j=j-k;
  124.                         k=k/2;
  125.                 }
  126.                 j=j+k;        //最高位为0
  127.         }
  128.         printf("\n\n=============================================================================================================================================================================================\n");
  129.         printf("输出倒序后的结果:\n");
  130.         outputArray();
  131. }


  132. /*遍历输出数组*/
  133. void outputArray()
  134. {
  135.         int num=0;
  136.         for(int i=0;i<size_x;i++)
  137.         {               
  138.        
  139.                 if(x[i].imaginaryPart<0)
  140.                 {
  141.                         printf("%.4f-j%.4f     ",x[i].realPart,fabs(x[i].imaginaryPart));
  142.                 }
  143.                 else
  144.                 {
  145.                         printf("%.4f+j%.4f     ",x[i].realPart,x[i].imaginaryPart);
  146.                 }
  147.                 num++;
  148.                 if(num==10){
  149.                         printf("\n");
  150.                         num=0;
  151.                 }
  152.         }
  153. }


  154. /*复数加法的定义*/
  155. void add(complexNumber a,complexNumber b,complexNumber *c)
  156. {
  157.         c->realPart=a.realPart+b.realPart;
  158.         c->imaginaryPart=a.imaginaryPart+b.imaginaryPart;
  159. }

  160. /*复数乘法的定义*/
  161. void mul(complexNumber a,complexNumber b,complexNumber *c)
  162. {
  163.         c->realPart=a.realPart*b.realPart - a.imaginaryPart*b.imaginaryPart;
  164.         c->imaginaryPart=a.realPart*b.imaginaryPart + a.imaginaryPart*b.realPart;
  165. }

  166. /*复数减法的定义*/
  167. void sub(complexNumber a,complexNumber b,complexNumber *c)
  168. {
  169.         c->realPart=a.realPart-b.realPart;
  170.         c->imaginaryPart=a.imaginaryPart-b.imaginaryPart;
  171. }


hearstnorman323 发表于 2023-3-2 11:12 | 显示全部楼层
FFT的最优算法是什么?              
1988020566 发表于 2023-3-2 11:27 | 显示全部楼层
如何调用FFT32汇编程序实现
tabmone 发表于 2023-3-2 11:39 | 显示全部楼层
怎么利用FFT算法对音频噪声进行处理
lzbf 发表于 2023-3-2 11:59 | 显示全部楼层
为什么FFT算法能够降低DFT的复杂度
usysm 发表于 2023-3-3 22:13 | 显示全部楼层
如何用FFT算法转换成相位差               
pentruman 发表于 2023-3-4 21:06 | 显示全部楼层
为什么用C语言的FFT算法后,结果都溢出了
alvpeg 发表于 2023-3-5 10:30 | 显示全部楼层
为什么用C语言的FFT算法后,结果都溢出了
plsbackup 发表于 2023-3-9 13:17 | 显示全部楼层
为什么用C语言的FFT算法后,结果都溢出了
yorkbarney 发表于 2023-3-9 13:28 | 显示全部楼层
C语言的fft和ifft程序工程文件,能共享一下吗?
modesty3jonah 发表于 2023-3-10 10:11 | 显示全部楼层
求FFT的C语言实现                 
hearstnorman323 发表于 2023-3-10 10:36 | 显示全部楼层
fft函数返回值是什么               
averyleigh 发表于 2023-3-10 10:52 | 显示全部楼层
信号处理有哪些算法可以实现?              
ulystronglll 发表于 2023-3-10 10:56 | 显示全部楼层
FFT运算就是快速傅里叶变换              
mmbs 发表于 2023-3-10 11:19 | 显示全部楼层
二阶滤波器用C语言怎么写              
lihuami 发表于 2023-3-10 11:23 | 显示全部楼层
在学信号系统,还没学到离散傅里叶变
febgxu 发表于 2023-3-10 13:08 | 显示全部楼层
如何用FFT得到谐波幅值频率和相位
louliana 发表于 2023-3-10 13:14 | 显示全部楼层
有谁用C语言实现过correl函数吗
robincotton 发表于 2023-3-10 13:50 | 显示全部楼层
FFT的最优算法是什么?              
houjiakai 发表于 2023-3-10 13:57 | 显示全部楼层
有用C语言实现的互相关算法的源代码吗
您需要登录后才可以回帖 登录 | 注册

本版积分规则

25

主题

2474

帖子

1

粉丝
快速回复 在线客服 返回列表 返回顶部