1 /* fft/hc_main.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 #include "gsl__config.h"
21 #include <stdlib.h>
22 #include <math.h>
23
24 #include "gsl_errno.h"
25 #include "gsl_complex.h"
26 #include "gsl_fft_halfcomplex.h"
27
28 #include "gsl_fft__hc_pass.h"
29
30 int
FUNCTION(gsl_fft_halfcomplex,backward)31 FUNCTION(gsl_fft_halfcomplex,backward) (BASE data[], const size_t stride,
32 const size_t n,
33 const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
34 TYPE(gsl_fft_real_workspace) * work)
35 {
36 int status = FUNCTION(gsl_fft_halfcomplex,transform) (data, stride, n, wavetable, work) ;
37 return status ;
38 }
39
40 int
FUNCTION(gsl_fft_halfcomplex,inverse)41 FUNCTION(gsl_fft_halfcomplex,inverse) (BASE data[], const size_t stride,
42 const size_t n,
43 const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
44 TYPE(gsl_fft_real_workspace) * work)
45 {
46 int status = FUNCTION(gsl_fft_halfcomplex,transform) (data, stride, n, wavetable, work);
47
48 if (status)
49 {
50 return status;
51 }
52
53 /* normalize inverse fft with 1/n */
54
55 {
56 const double norm = 1.0 / n;
57 size_t i;
58 for (i = 0; i < n; i++)
59 {
60 data[stride*i] *= norm;
61 }
62 }
63 return status;
64 }
65
66 int
FUNCTION(gsl_fft_halfcomplex,transform)67 FUNCTION(gsl_fft_halfcomplex,transform) (BASE data[], const size_t stride, const size_t n,
68 const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
69 TYPE(gsl_fft_real_workspace) * work)
70 {
71 BASE * const scratch = work->scratch;
72
73 BASE * in;
74 BASE * out;
75 size_t istride, ostride ;
76
77
78 size_t factor, product, q;
79 size_t i;
80 size_t nf;
81 int state;
82 int tskip;
83 TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4;
84
85 if (n == 0)
86 {
87 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
88 }
89
90 if (n == 1)
91 { /* FFT of one data point is the identity */
92 return 0;
93 }
94
95 if (n != wavetable->n)
96 {
97 GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
98 }
99
100 if (n != work->n)
101 {
102 GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
103 }
104
105 nf = wavetable->nf;
106 product = 1;
107 state = 0;
108
109 for (i = 0; i < nf; i++)
110 {
111 factor = wavetable->factor[i];
112 product *= factor;
113 q = n / product;
114
115 tskip = (q + 1) / 2 - 1;
116
117 if (state == 0)
118 {
119 in = data;
120 istride = stride;
121 out = scratch;
122 ostride = 1;
123 state = 1;
124 }
125 else
126 {
127 in = scratch;
128 istride = 1;
129 out = data;
130 ostride = stride;
131 state = 0;
132 }
133
134 if (factor == 2)
135 {
136 twiddle1 = wavetable->twiddle[i];
137 FUNCTION(fft_halfcomplex,pass_2) (in, istride, out, ostride,
138 product, n, twiddle1);
139 }
140 else if (factor == 3)
141 {
142 twiddle1 = wavetable->twiddle[i];
143 twiddle2 = twiddle1 + tskip;
144 FUNCTION(fft_halfcomplex,pass_3) (in, istride, out, ostride,
145 product, n, twiddle1, twiddle2);
146 }
147 else if (factor == 4)
148 {
149 twiddle1 = wavetable->twiddle[i];
150 twiddle2 = twiddle1 + tskip;
151 twiddle3 = twiddle2 + tskip;
152 FUNCTION(fft_halfcomplex,pass_4) (in, istride, out, ostride,
153 product, n, twiddle1, twiddle2,
154 twiddle3);
155 }
156 else if (factor == 5)
157 {
158 twiddle1 = wavetable->twiddle[i];
159 twiddle2 = twiddle1 + tskip;
160 twiddle3 = twiddle2 + tskip;
161 twiddle4 = twiddle3 + tskip;
162 FUNCTION(fft_halfcomplex,pass_5) (in, istride, out, ostride,
163 product, n, twiddle1, twiddle2,
164 twiddle3, twiddle4);
165 }
166 else
167 {
168 twiddle1 = wavetable->twiddle[i];
169 FUNCTION(fft_halfcomplex,pass_n) (in, istride, out, ostride,
170 factor, product, n, twiddle1);
171 }
172 }
173
174 if (state == 1) /* copy results back from scratch to data */
175 {
176 for (i = 0; i < n; i++)
177 {
178 data[stride*i] = scratch[i] ;
179 }
180 }
181
182 return 0;
183
184 }
185
186
187