傅立叶变换与傅立叶反变换的C语言实现

  1. #include <math.h>  
  2. #include <malloc.h>  
  3. #define pi (double) 3.14159265359  
  4. /*复数的定义*/  
  5. typedef struct  
  6. {  
  7.     double re;  
  8.     double im;  
  9. }COMPLEX;  
  10. /*复数的加运算*/  
  11. COMPLEX Add(COMPLEX c1, COMPLEX c2)  
  12. {  
  13.     COMPLEX c;  
  14.     c.re = c1.re + c2.re;  
  15.     c.im = c1.im + c2.im;  
  16.     return c;  
  17. }  
  18. /*负数的减运算*/  
  19. COMPLEX Sub(COMPLEX c1, COMPLEX c2)  
  20. {  
  21.     COMPLEX c;  
  22.     c.re = c1.re - c2.re;  
  23.     c.im = c1.im - c2.im;  
  24.     return c;  
  25. }  
  26. /*复数的乘运算*/  
  27. COMPLEX Mul(COMPLEX c1, COMPLEX c2)  
  28. {  
  29.     COMPLEX c;  
  30.     c.re = c1.re*c2.re - c1.im*c2.im;  
  31.     c.im = c1.re*c2.im + c1.im*c2.re;  
  32.     return c;  
  33. }  
  34. /*快速傅立叶变换 
  35. TD为时域值,FD为频域值,power为2的幂数*/  
  36. void FFT(COMPLEX *TD, COMPLEX *FD, int power)  
  37. {  
  38.     int count;  
  39.     int i,j,k,bfsize,p;  
  40.     double angle;  
  41.     COMPLEX *W,*X1,*X2,*X;  
  42.     /*计算傅立叶变换点数*/  
  43.     count=1<<power;  
  44.     /*分配运算器所需存储器*/  
  45.     W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);  
  46.     X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);  
  47.     X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);  
  48.     /*计算加权系数*/  
  49.     for(i=0;i<count/2;i++)  
  50.     {  
  51.         angle=-i*pi*2/count;  
  52.         W[i].re=cos(angle);  
  53.         W[i].im=sin(angle);  
  54.     }  
  55.     /*将时域点写入存储器*/  
  56.     memcpy(X1, TD, sizeof(COMPLEX)*count);  
  57.     /*蝶形运算*/  
  58.     for(k=0; k<power; k++)  
  59.     {  
  60.         for(j=0;j<1<<k;j++)  
  61.         {  
  62.             bfsize=1<<power-k;  
  63.             for(i=0;i<bfsize/2;i++)  
  64.             {  
  65.                 p=j*bfsize;  
  66.                 X2[i+p]=Add(X1[i+p], X1[i+p+bfsize/2]);  
  67.                 X2[i+p+bfsize/2]=Mul(Sub(X1[i+p], X1[i+p+bfsize/2]),W[i*(1<<k)]);  
  68.             }  
  69.         }  
  70.         X=X1;  
  71.         X1=X2;  
  72.         X2=X;  
  73.     }  
  74.     /*重新排序*/  
  75.     for(j=0;j<count;j++)  
  76.     {  
  77.         p=0;  
  78.         for(i=0;i<power;i++)  
  79.         {  
  80.             if(j&(1<<i))  
  81.             p+=1<<power-i-1;  
  82.         }  
  83.         FD[j]=X1[p];  
  84.     }  
  85.     /*释放存储器*/  
  86.     free(W);  
  87.     free(X1);  
  88.     free(X2);  
  89. }  
  90. /*快速傅立叶反变换,利用快速傅立叶变换 
  91. FD为频域值,TD为时域值,power为2的幂数*/  
  92. void IFFT(COMPLEX *FD, COMPLEX *TD, int power)  
  93. {  
  94.     int i,count;  
  95.     COMPLEX *x;  
  96.     /*计算傅立叶反变换点数*/  
  97.     count=1<<power;  
  98.     /*分配运算所需存储器*/  
  99.     x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);  
  100.     /*将频域点写入存储器*/  
  101.     memcpy(x,FD,sizeof(COMPLEX)*count);  
  102.     /*求频域点的共轭*/  
  103.     for(i=0;i<count;i++)  
  104.     {  
  105.         x[i].im=-x[i].im;  
  106.     }  
  107.     /*调用快速傅立叶变换*/  
  108.     FFT(x,TD,power);  
  109.     /*求时域点的共轭*/  
  110.     for(i=0;i<count;i++)  
  111.     {  
  112.         TD[i].re/=count;  
  113.         TD[i].im=-TD[i].im/count;  
  114.     }  
  115.     /*释放存储器*/  
  116.     free(x);  
  117. }  
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值