1 /*
2 * This file is part of libfftpack.
3 *
4 * libfftpack is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * libfftpack is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with libfftpack; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
19 /*
20 * libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
21 * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22 * (DLR).
23 */
24
25 /*
26 * Copyright (C) 2005 Max-Planck-Society
27 * \author Martin Reinecke
28 */
29
30 #include <stdlib.h>
31 #include <math.h>
32 #include <string.h>
33 #include "bluestein.h"
34 #include "fftpack.h"
35 #include "ls_fft.h"
36
make_complex_plan(size_t length)37 complex_plan make_complex_plan (size_t length)
38 {
39 complex_plan plan = RALLOC(complex_plan_i,1);
40 size_t pfsum = prime_factor_sum(length);
41 double comp1 = (double)(length*pfsum);
42 double comp2 = 2*3*length*log(3.*length);
43 comp2*=3.; /* fudge factor that appears to give good overall performance */
44 plan->length=length;
45 plan->bluestein = (comp2<comp1);
46 if (plan->bluestein)
47 bluestein_i (length,&(plan->work),&(plan->worksize));
48 else
49 {
50 plan->worksize=4*length+15;
51 plan->work=RALLOC(double,4*length+15);
52 cffti(length, plan->work);
53 }
54 return plan;
55 }
56
copy_complex_plan(complex_plan plan)57 complex_plan copy_complex_plan (complex_plan plan)
58 {
59 if (!plan) return NULL;
60 {
61 complex_plan newplan = RALLOC(complex_plan_i,1);
62 *newplan = *plan;
63 newplan->work=RALLOC(double,newplan->worksize);
64 memcpy(newplan->work,plan->work,sizeof(double)*newplan->worksize);
65 return newplan;
66 }
67 }
68
kill_complex_plan(complex_plan plan)69 void kill_complex_plan (complex_plan plan)
70 {
71 DEALLOC(plan->work);
72 DEALLOC(plan);
73 }
74
complex_plan_forward(complex_plan plan,double * data)75 void complex_plan_forward (complex_plan plan, double *data)
76 {
77 if (plan->bluestein)
78 bluestein (plan->length, data, plan->work, -1);
79 else
80 cfftf (plan->length, data, plan->work);
81 }
82
complex_plan_backward(complex_plan plan,double * data)83 void complex_plan_backward (complex_plan plan, double *data)
84 {
85 if (plan->bluestein)
86 bluestein (plan->length, data, plan->work, 1);
87 else
88 cfftb (plan->length, data, plan->work);
89 }
90
91
make_real_plan(size_t length)92 real_plan make_real_plan (size_t length)
93 {
94 real_plan plan = RALLOC(real_plan_i,1);
95 size_t pfsum = prime_factor_sum(length);
96 double comp1 = .5*length*pfsum;
97 double comp2 = 2*3*length*log(3.*length);
98 comp2*=3; /* fudge factor that appears to give good overall performance */
99 plan->length=length;
100 plan->bluestein = (comp2<comp1);
101 if (plan->bluestein)
102 bluestein_i (length,&(plan->work),&(plan->worksize));
103 else
104 {
105 plan->worksize=2*length+15;
106 plan->work=RALLOC(double,2*length+15);
107 rffti(length, plan->work);
108 }
109 return plan;
110 }
111
copy_real_plan(real_plan plan)112 real_plan copy_real_plan (real_plan plan)
113 {
114 if (!plan) return NULL;
115 {
116 real_plan newplan = RALLOC(real_plan_i,1);
117 *newplan = *plan;
118 newplan->work=RALLOC(double,newplan->worksize);
119 memcpy(newplan->work,plan->work,sizeof(double)*newplan->worksize);
120 return newplan;
121 }
122 }
123
kill_real_plan(real_plan plan)124 void kill_real_plan (real_plan plan)
125 {
126 DEALLOC(plan->work);
127 DEALLOC(plan);
128 }
129
real_plan_forward_fftpack(real_plan plan,double * data)130 void real_plan_forward_fftpack (real_plan plan, double *data)
131 {
132 if (plan->bluestein)
133 {
134 size_t m;
135 size_t n=plan->length;
136 double *tmp = RALLOC(double,2*n);
137 for (m=0; m<n; ++m)
138 {
139 tmp[2*m] = data[m];
140 tmp[2*m+1] = 0.;
141 }
142 bluestein(n,tmp,plan->work,-1);
143 data[0] = tmp[0];
144 memcpy (data+1, tmp+2, (n-1)*sizeof(double));
145 DEALLOC(tmp);
146 }
147 else
148 rfftf (plan->length, data, plan->work);
149 }
150
fftpack2halfcomplex(double * data,size_t n)151 static void fftpack2halfcomplex (double *data, size_t n)
152 {
153 size_t m;
154 double *tmp = RALLOC(double,n);
155 tmp[0]=data[0];
156 for (m=1; m<(n+1)/2; ++m)
157 {
158 tmp[m]=data[2*m-1];
159 tmp[n-m]=data[2*m];
160 }
161 if (!(n&1))
162 tmp[n/2]=data[n-1];
163 memcpy (data,tmp,n*sizeof(double));
164 DEALLOC(tmp);
165 }
166
halfcomplex2fftpack(double * data,size_t n)167 static void halfcomplex2fftpack (double *data, size_t n)
168 {
169 size_t m;
170 double *tmp = RALLOC(double,n);
171 tmp[0]=data[0];
172 for (m=1; m<(n+1)/2; ++m)
173 {
174 tmp[2*m-1]=data[m];
175 tmp[2*m]=data[n-m];
176 }
177 if (!(n&1))
178 tmp[n-1]=data[n/2];
179 memcpy (data,tmp,n*sizeof(double));
180 DEALLOC(tmp);
181 }
182
real_plan_forward_fftw(real_plan plan,double * data)183 void real_plan_forward_fftw (real_plan plan, double *data)
184 {
185 real_plan_forward_fftpack (plan, data);
186 fftpack2halfcomplex (data,plan->length);
187 }
188
real_plan_backward_fftpack(real_plan plan,double * data)189 void real_plan_backward_fftpack (real_plan plan, double *data)
190 {
191 if (plan->bluestein)
192 {
193 size_t m;
194 size_t n=plan->length;
195 double *tmp = RALLOC(double,2*n);
196 tmp[0]=data[0];
197 tmp[1]=0.;
198 memcpy (tmp+2,data+1, (n-1)*sizeof(double));
199 if ((n&1)==0) tmp[n+1]=0.;
200 for (m=2; m<n; m+=2)
201 {
202 tmp[2*n-m]=tmp[m];
203 tmp[2*n-m+1]=-tmp[m+1];
204 }
205 bluestein (n, tmp, plan->work, 1);
206 for (m=0; m<n; ++m)
207 data[m] = tmp[2*m];
208 DEALLOC(tmp);
209 }
210 else
211 rfftb (plan->length, data, plan->work);
212 }
213
real_plan_backward_fftw(real_plan plan,double * data)214 void real_plan_backward_fftw (real_plan plan, double *data)
215 {
216 halfcomplex2fftpack (data,plan->length);
217 real_plan_backward_fftpack (plan, data);
218 }
219
real_plan_forward_c(real_plan plan,double * data)220 void real_plan_forward_c (real_plan plan, double *data)
221 {
222 size_t m;
223 size_t n=plan->length;
224
225 if (plan->bluestein)
226 {
227 for (m=1; m<2*n; m+=2)
228 data[m]=0;
229 bluestein (plan->length, data, plan->work, -1);
230 data[1]=0;
231 for (m=2; m<n; m+=2)
232 {
233 double avg;
234 avg = 0.5*(data[2*n-m]+data[m]);
235 data[2*n-m] = data[m] = avg;
236 avg = 0.5*(data[2*n-m+1]-data[m+1]);
237 data[2*n-m+1] = avg;
238 data[m+1] = -avg;
239 }
240 if ((n&1)==0) data[n+1] = 0.;
241 }
242 else
243 {
244 /* using "m+m" instead of "2*m" to avoid a nasty bug in Intel's compiler */
245 for (m=0; m<n; ++m) data[m+1] = data[m+m];
246 rfftf (n, data+1, plan->work);
247 data[0] = data[1];
248 data[1] = 0;
249 for (m=2; m<n; m+=2)
250 {
251 data[2*n-m] = data[m];
252 data[2*n-m+1] = -data[m+1];
253 }
254 if ((n&1)==0) data[n+1] = 0.;
255 }
256 }
257
real_plan_backward_c(real_plan plan,double * data)258 void real_plan_backward_c (real_plan plan, double *data)
259 {
260 size_t n=plan->length;
261
262 if (plan->bluestein)
263 {
264 size_t m;
265 data[1]=0;
266 for (m=2; m<n; m+=2)
267 {
268 double avg;
269 avg = 0.5*(data[2*n-m]+data[m]);
270 data[2*n-m] = data[m] = avg;
271 avg = 0.5*(data[2*n-m+1]-data[m+1]);
272 data[2*n-m+1] = avg;
273 data[m+1] = -avg;
274 }
275 if ((n&1)==0) data[n+1] = 0.;
276 bluestein (plan->length, data, plan->work, 1);
277 for (m=1; m<2*n; m+=2)
278 data[m]=0;
279 }
280 else
281 {
282 ptrdiff_t m;
283 data[1] = data[0];
284 rfftb (n, data+1, plan->work);
285 for (m=n-1; m>=0; --m)
286 {
287 data[2*m] = data[m+1];
288 data[2*m+1] = 0.;
289 }
290 }
291 }
292