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