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 <config.h>
21 #include <stdlib.h>
22 #include <math.h>
23
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_complex.h>
26 #include <gsl/gsl_fft_halfcomplex.h>
27
28 #include "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 product_1;
83 int tskip;
84 TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4;
85
86 if (n == 0)
87 {
88 GSL_ERROR ("length n must be positive integer", GSL_EDOM);
89 }
90
91 if (n == 1)
92 { /* FFT of one data point is the identity */
93 return 0;
94 }
95
96 if (n != wavetable->n)
97 {
98 GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
99 }
100
101 if (n != work->n)
102 {
103 GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
104 }
105
106 nf = wavetable->nf;
107 product = 1;
108 state = 0;
109
110 for (i = 0; i < nf; i++)
111 {
112 factor = wavetable->factor[i];
113 product_1 = product;
114 product *= factor;
115 q = n / product;
116
117 tskip = (q + 1) / 2 - 1;
118
119 if (state == 0)
120 {
121 in = data;
122 istride = stride;
123 out = scratch;
124 ostride = 1;
125 state = 1;
126 }
127 else
128 {
129 in = scratch;
130 istride = 1;
131 out = data;
132 ostride = stride;
133 state = 0;
134 }
135
136 if (factor == 2)
137 {
138 twiddle1 = wavetable->twiddle[i];
139 FUNCTION(fft_halfcomplex,pass_2) (in, istride, out, ostride,
140 product, n, twiddle1);
141 }
142 else if (factor == 3)
143 {
144 twiddle1 = wavetable->twiddle[i];
145 twiddle2 = twiddle1 + tskip;
146 FUNCTION(fft_halfcomplex,pass_3) (in, istride, out, ostride,
147 product, n, twiddle1, twiddle2);
148 }
149 else if (factor == 4)
150 {
151 twiddle1 = wavetable->twiddle[i];
152 twiddle2 = twiddle1 + tskip;
153 twiddle3 = twiddle2 + tskip;
154 FUNCTION(fft_halfcomplex,pass_4) (in, istride, out, ostride,
155 product, n, twiddle1, twiddle2,
156 twiddle3);
157 }
158 else if (factor == 5)
159 {
160 twiddle1 = wavetable->twiddle[i];
161 twiddle2 = twiddle1 + tskip;
162 twiddle3 = twiddle2 + tskip;
163 twiddle4 = twiddle3 + tskip;
164 FUNCTION(fft_halfcomplex,pass_5) (in, istride, out, ostride,
165 product, n, twiddle1, twiddle2,
166 twiddle3, twiddle4);
167 }
168 else
169 {
170 twiddle1 = wavetable->twiddle[i];
171 FUNCTION(fft_halfcomplex,pass_n) (in, istride, out, ostride,
172 factor, product, n, twiddle1);
173 }
174 }
175
176 if (state == 1) /* copy results back from scratch to data */
177 {
178 for (i = 0; i < n; i++)
179 {
180 data[stride*i] = scratch[i] ;
181 }
182 }
183
184 return 0;
185
186 }
187
188
189