1 #include <math.h>
2 #include <string.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5
6 typedef struct complex { float r , i ; } complex ;
7
8 /*--------------------------------------------------------------------
9 This is the older FFT routine by AJ and RWC,
10 from which the routines in csfft.c were derived.
11 ----------------------------------------------------------------------*/
12
13 static complex * csplus = NULL , * csminus = NULL ; /* trig consts */
14 static int nold = -666 ;
15
16 #undef PI
17 #define PI (3.141592653589793238462643) /* should be close enough */
18
csfft_trigconsts(int idim)19 static void csfft_trigconsts( int idim )
20 {
21 register unsigned int m, n, i0, i1, i2, i3, k;
22 register complex *r0, *csp;
23 register float co, si, f0, f1, f2, f3, f4;
24 double al;
25
26 if( idim == nold ) return ;
27
28 if( idim > nold ){
29 if( csplus != NULL ){ free(csplus) ; free(csminus) ; }
30 csplus = (complex *) malloc( sizeof(complex) * idim ) ;
31 csminus = (complex *) malloc( sizeof(complex) * idim ) ;
32 if( csplus == NULL || csminus == NULL ){
33 fprintf(stderr,"\n*** fft cannot malloc space! ***\n"); exit(-1) ;
34 }
35 }
36 nold = n = idim ;
37
38 f1 = 1.0 ; /* csplus init */
39 m = 1;
40 k = 0;
41 while (n > m) {
42 i3 = m << 1;
43 f2 = m;
44 al = f1*PI/f2;
45 co = cos(al);
46 si = sin(al);
47 (csplus + k)->r = 1.;
48 (csplus + k)->i = 0.;
49 for (i0=0; i0 < m; i0++) {
50 k++;
51 csp = csplus + k;
52 r0 = csp - 1;
53 csp->r = r0->r * co - r0->i * si;
54 csp->i = r0->i * co + r0->r * si;
55 }
56 m = i3;
57 }
58
59 f1 = -1.0 ; /* csminus init */
60 m = 1;
61 k = 0;
62 while (n > m) {
63 i3 = m << 1;
64 f2 = m;
65 al = f1*PI/f2;
66 co = cos(al);
67 si = sin(al);
68 (csminus + k)->r = 1.;
69 (csminus + k)->i = 0.;
70 for (i0=0; i0 < m; i0++) {
71 k++;
72 csp = csminus + k;
73 r0 = csp - 1;
74 csp->r = r0->r * co - r0->i * si;
75 csp->i = r0->i * co + r0->r * si;
76 }
77 m = i3;
78 }
79 return ;
80 }
81
82 /*--------------------------------------------------------------------
83 Complex-to-complex FFT in place:
84 mode = -1 or +1 (NO SCALING ON INVERSE!)
85 idim = dimension (power of 2)
86 xc = input/output array
87 Automagically re-initializes itself when idim changes from
88 previous call. By AJJesmanowicz, modified by RWCox.
89 ----------------------------------------------------------------------*/
90
csfft(int mode,int idim,complex * xc)91 void csfft( int mode , int idim , complex * xc )
92 {
93 register unsigned int m, n, i0, i1, i2, i3, k;
94 register complex *r0, *r1, *csp;
95 register float co, si, f0, f1, f2, f3, f4;
96 double al;
97
98 /** perhaps initialize **/
99
100 if( nold != idim ) csfft_trigconsts( idim ) ;
101
102 /** Main loop starts here **/
103
104 n = idim;
105 i2 = idim >> 1;
106 i1 = 0;
107 csp = (mode > 0) ? csplus : csminus ; /* choose const array */
108
109 for (i0=0; i0 < n; i0 ++) {
110 if ( i1 > i0 ) {
111 r0 = xc + i0;
112 r1 = xc + i1;
113 f1 = r0->r;
114 f2 = r0->i;
115 r0->r = r1->r;
116 r0->i = r1->i;
117 r1->r = f1;
118 r1->i = f2;
119 }
120 m = i2;
121 while ( m && !(i1 < m) ) {
122 i1 -= m;
123 m >>= 1;
124 }
125 i1 += m;
126 }
127
128 m = 1;
129 k = 0;
130 while (n > m) {
131 i3 = m << 1;
132 for (i0=0; i0 < m; i0 ++) {
133 co = (csp + k)->r;
134 si = (csp + k)->i;
135 for (i1=i0; i1 < n; i1 += i3) {
136 r0 = xc + i1;
137 r1 = r0 + m;
138 f1 = r1->r * co;
139 f2 = r1->i * si;
140 f3 = r1->r * si;
141 f4 = r1->i * co;
142 f1 -= f2;
143 f3 += f4;
144 f2 = r0->r;
145 f4 = r0->i;
146 r1->r = f2 - f1;
147 r1->i = f4 - f3;
148 r0->r = f2 + f1;
149 r0->i = f4 + f3;
150 }
151 k++;
152 }
153 m = i3;
154 }
155
156 #ifdef SCALE_INVERSE
157 if (mode > 0) {
158 f0 = 1.0 / idim ;
159 i0 = 0;
160 i1 = 1;
161 while (i0 < n) {
162 r0 = xc + i0;
163 r1 = xc + i1;
164 f1 = r0->r;
165 f2 = r0->i;
166 f3 = r1->r;
167 f4 = r1->i;
168 f1 *= f0;
169 f2 *= f0;
170 f3 *= f0;
171 f4 *= f0;
172 r0->r = f1;
173 r0->i = f2;
174 r1->r = f3;
175 r1->i = f4;
176 i0 += 2;
177 i1 += 2;
178 }
179 }
180 #endif
181
182 return ;
183 }
184