傅立叶变换与傅立叶反变换的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. }