1 /* fft/c_radix2.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 int
FUNCTION(gsl_fft_complex,radix2_forward)21 FUNCTION(gsl_fft_complex,radix2_forward) (TYPE(gsl_complex_packed_array) data,
22 const size_t stride, const size_t n)
23 {
24 gsl_fft_direction sign = gsl_fft_forward;
25 int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
26 return status;
27 }
28
29 int
FUNCTION(gsl_fft_complex,radix2_backward)30 FUNCTION(gsl_fft_complex,radix2_backward) (TYPE(gsl_complex_packed_array) data,
31 const size_t stride, const size_t n)
32 {
33 gsl_fft_direction sign = gsl_fft_backward;
34 int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
35 return status;
36 }
37
38 int
FUNCTION(gsl_fft_complex,radix2_inverse)39 FUNCTION(gsl_fft_complex,radix2_inverse) (TYPE(gsl_complex_packed_array) data,
40 const size_t stride, const size_t n)
41 {
42 gsl_fft_direction sign = gsl_fft_backward;
43 int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
44
45 if (status)
46 {
47 return status;
48 }
49
50 /* normalize inverse fft with 1/n */
51
52 {
53 const ATOMIC norm = 1.0 / n;
54 size_t i;
55 for (i = 0; i < n; i++)
56 {
57 REAL(data,stride,i) *= norm;
58 IMAG(data,stride,i) *= norm;
59 }
60 }
61
62 return status;
63 }
64
65
66
67 int
FUNCTION(gsl_fft_complex,radix2_transform)68 FUNCTION(gsl_fft_complex,radix2_transform) (TYPE(gsl_complex_packed_array) data,
69 const size_t stride,
70 const size_t n,
71 const gsl_fft_direction sign)
72 {
73 int result ;
74 size_t dual;
75 size_t bit;
76 size_t logn = 0;
77
78 if (n == 1) /* identity operation */
79 {
80 return 0 ;
81 }
82
83 /* make sure that n is a power of 2 */
84
85 result = fft_binary_logn(n) ;
86
87 if (result == -1)
88 {
89 GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
90 }
91 else
92 {
93 logn = result ;
94 }
95
96 /* bit reverse the ordering of input data for decimation in time algorithm */
97
98 (void) FUNCTION(fft_complex,bitreverse_order) (data, stride, n, logn) ;
99
100 /* apply fft recursion */
101
102 dual = 1;
103
104 for (bit = 0; bit < logn; bit++)
105 {
106 ATOMIC w_real = 1.0;
107 ATOMIC w_imag = 0.0;
108
109 const double theta = 2.0 * ((int) sign) * M_PI / (2.0 * (double) dual);
110
111 const ATOMIC s = sin (theta);
112 const ATOMIC t = sin (theta / 2.0);
113 const ATOMIC s2 = 2.0 * t * t;
114
115 size_t a, b;
116
117 /* a = 0 */
118
119 for (b = 0; b < n; b += 2 * dual)
120 {
121 const size_t i = b ;
122 const size_t j = b + dual;
123
124 const ATOMIC z1_real = REAL(data,stride,j) ;
125 const ATOMIC z1_imag = IMAG(data,stride,j) ;
126
127 const ATOMIC wd_real = z1_real ;
128 const ATOMIC wd_imag = z1_imag ;
129
130 REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
131 IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
132 REAL(data,stride,i) += wd_real;
133 IMAG(data,stride,i) += wd_imag;
134 }
135
136 /* a = 1 .. (dual-1) */
137
138 for (a = 1; a < dual; a++)
139 {
140
141 /* trignometric recurrence for w-> exp(i theta) w */
142
143 {
144 const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
145 const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
146 w_real = tmp_real;
147 w_imag = tmp_imag;
148 }
149
150 for (b = 0; b < n; b += 2 * dual)
151 {
152 const size_t i = b + a;
153 const size_t j = b + a + dual;
154
155 const ATOMIC z1_real = REAL(data,stride,j) ;
156 const ATOMIC z1_imag = IMAG(data,stride,j) ;
157
158 const ATOMIC wd_real = w_real * z1_real - w_imag * z1_imag;
159 const ATOMIC wd_imag = w_real * z1_imag + w_imag * z1_real;
160
161 REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
162 IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
163 REAL(data,stride,i) += wd_real;
164 IMAG(data,stride,i) += wd_imag;
165 }
166 }
167 dual *= 2;
168 }
169
170 return 0;
171
172 }
173
174
175 int
FUNCTION(gsl_fft_complex,radix2_dif_forward)176 FUNCTION(gsl_fft_complex,radix2_dif_forward) (TYPE(gsl_complex_packed_array) data,
177 const size_t stride,
178 const size_t n)
179 {
180 gsl_fft_direction sign = gsl_fft_forward;
181 int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
182 return status;
183 }
184
185 int
FUNCTION(gsl_fft_complex,radix2_dif_backward)186 FUNCTION(gsl_fft_complex,radix2_dif_backward) (TYPE(gsl_complex_packed_array) data,
187 const size_t stride,
188 const size_t n)
189 {
190 gsl_fft_direction sign = gsl_fft_backward;
191 int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
192 return status;
193 }
194
195 int
FUNCTION(gsl_fft_complex,radix2_dif_inverse)196 FUNCTION(gsl_fft_complex,radix2_dif_inverse) (TYPE(gsl_complex_packed_array) data,
197 const size_t stride,
198 const size_t n)
199 {
200 gsl_fft_direction sign = gsl_fft_backward;
201 int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
202
203 if (status)
204 {
205 return status;
206 }
207
208 /* normalize inverse fft with 1/n */
209
210 {
211 const ATOMIC norm = 1.0 / n;
212 size_t i;
213 for (i = 0; i < n; i++)
214 {
215 REAL(data,stride,i) *= norm;
216 IMAG(data,stride,i) *= norm;
217 }
218 }
219
220 return status;
221 }
222
223 int
FUNCTION(gsl_fft_complex,radix2_dif_transform)224 FUNCTION(gsl_fft_complex,radix2_dif_transform) (TYPE(gsl_complex_packed_array) data,
225 const size_t stride,
226 const size_t n,
227 const gsl_fft_direction sign)
228 {
229 int result ;
230 size_t dual;
231 size_t bit;
232 size_t logn = 0;
233
234 if (n == 1) /* identity operation */
235 {
236 return 0 ;
237 }
238
239 /* make sure that n is a power of 2 */
240
241 result = fft_binary_logn(n) ;
242
243 if (result == -1)
244 {
245 GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
246 }
247 else
248 {
249 logn = result ;
250 }
251
252 /* apply fft recursion */
253
254 dual = n / 2;
255
256 for (bit = 0; bit < logn; bit++)
257 {
258 ATOMIC w_real = 1.0;
259 ATOMIC w_imag = 0.0;
260
261 const double theta = 2.0 * ((int) sign) * M_PI / ((double) (2 * dual));
262
263 const ATOMIC s = sin (theta);
264 const ATOMIC t = sin (theta / 2.0);
265 const ATOMIC s2 = 2.0 * t * t;
266
267 size_t a, b;
268
269 for (b = 0; b < dual; b++)
270 {
271 for (a = 0; a < n; a+= 2 * dual)
272 {
273 const size_t i = b + a;
274 const size_t j = b + a + dual;
275
276 const ATOMIC t1_real = REAL(data,stride,i) + REAL(data,stride,j);
277 const ATOMIC t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
278 const ATOMIC t2_real = REAL(data,stride,i) - REAL(data,stride,j);
279 const ATOMIC t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
280
281 REAL(data,stride,i) = t1_real;
282 IMAG(data,stride,i) = t1_imag;
283 REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
284 IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
285 }
286
287 /* trignometric recurrence for w-> exp(i theta) w */
288
289 {
290 const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
291 const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
292 w_real = tmp_real;
293 w_imag = tmp_imag;
294 }
295 }
296 dual /= 2;
297 }
298
299 /* bit reverse the ordering of output data for decimation in
300 frequency algorithm */
301
302 (void) FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
303
304 return 0;
305
306 }
307
308
309
310
311
312
313
314
315