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 fftpack.c : A set of FFT routines in C.
27 Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber
28 (Version 4, 1985).
29
30 C port by Martin Reinecke (2010)
31 */
32
33 #ifdef BACKWARD
34 #define PSIGN +
35 #define PMSIGNC(a,b,c,d) { a.r=c.r+d.r; a.i=c.i+d.i; b.r=c.r-d.r; b.i=c.i-d.i; }
36 /* a = b*c */
37 #define MULPMSIGNC(a,b,c) { a.r=b.r*c.r-b.i*c.i; a.i=b.r*c.i+b.i*c.r; }
38 #else
39 #define PSIGN -
40 #define PMSIGNC(a,b,c,d) { a.r=c.r-d.r; a.i=c.i-d.i; b.r=c.r+d.r; b.i=c.i+d.i; }
41 /* a = conj(b)*c */
42 #define MULPMSIGNC(a,b,c) { a.r=b.r*c.r+b.i*c.i; a.i=b.r*c.i-b.i*c.r; }
43 #endif
44
45 static void X(2) (size_t ido, size_t l1, const cmplx *cc, cmplx *ch,
46 const cmplx *wa)
47 {
48 const size_t cdim=2;
49 size_t k,i;
50 cmplx t;
51 if (ido==1)
52 for (k=0;k<l1;++k)
53 PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k))
54 else
55 for (k=0;k<l1;++k)
56 for (i=0;i<ido;++i)
57 {
58 PMC (CH(i,k,0),t,CC(i,0,k),CC(i,1,k))
59 MULPMSIGNC (CH(i,k,1),WA(0,i),t)
60 }
61 }
62
63 static void X(3)(size_t ido, size_t l1, const cmplx *cc, cmplx *ch,
64 const cmplx *wa)
65 {
66 const size_t cdim=3;
67 static const double taur=-0.5, taui= PSIGN 0.86602540378443864676;
68 size_t i, k;
69 cmplx c2, c3, d2, d3, t2;
70
71 if (ido==1)
72 for (k=0; k<l1; ++k)
73 {
74 PMC (t2,c3,CC(0,1,k),CC(0,2,k))
75 ADDC (CH(0,k,0),t2,CC(0,0,k))
76 SCALEC(t2,taur)
77 ADDC(c2,CC(0,0,k),t2)
78 SCALEC(c3,taui)
79 CONJFLIPC(c3)
80 PMC(CH(0,k,1),CH(0,k,2),c2,c3)
81 }
82 else
83 for (k=0; k<l1; ++k)
84 for (i=0; i<ido; ++i)
85 {
86 PMC (t2,c3,CC(i,1,k),CC(i,2,k))
87 ADDC (CH(i,k,0),t2,CC(i,0,k))
88 SCALEC(t2,taur)
89 ADDC(c2,CC(i,0,k),t2)
90 SCALEC(c3,taui)
91 CONJFLIPC(c3)
92 PMC(d2,d3,c2,c3)
93 MULPMSIGNC(CH(i,k,1),WA(0,i),d2)
94 MULPMSIGNC(CH(i,k,2),WA(1,i),d3)
95 }
96 }
97
98 static void X(4)(size_t ido, size_t l1, const cmplx *cc, cmplx *ch,
99 const cmplx *wa)
100 {
101 const size_t cdim=4;
102 size_t i, k;
103 cmplx c2, c3, c4, t1, t2, t3, t4;
104
105 if (ido==1)
106 for (k=0; k<l1; ++k)
107 {
108 PMC(t2,t1,CC(0,0,k),CC(0,2,k))
109 PMC(t3,t4,CC(0,1,k),CC(0,3,k))
110 CONJFLIPC(t4)
111 PMC(CH(0,k,0),CH(0,k,2),t2,t3)
112 PMSIGNC (CH(0,k,1),CH(0,k,3),t1,t4)
113 }
114 else
115 for (k=0; k<l1; ++k)
116 for (i=0; i<ido; ++i)
117 {
118 PMC(t2,t1,CC(i,0,k),CC(i,2,k))
119 PMC(t3,t4,CC(i,1,k),CC(i,3,k))
120 CONJFLIPC(t4)
121 PMC(CH(i,k,0),c3,t2,t3)
122 PMSIGNC (c2,c4,t1,t4)
123 MULPMSIGNC (CH(i,k,1),WA(0,i),c2)
124 MULPMSIGNC (CH(i,k,2),WA(1,i),c3)
125 MULPMSIGNC (CH(i,k,3),WA(2,i),c4)
126 }
127 }
128
129 static void X(5)(size_t ido, size_t l1, const cmplx *cc, cmplx *ch,
130 const cmplx *wa)
131 {
132 const size_t cdim=5;
133 static const double tr11= 0.3090169943749474241,
134 ti11= PSIGN 0.95105651629515357212,
135 tr12=-0.8090169943749474241,
136 ti12= PSIGN 0.58778525229247312917;
137 size_t i, k;
138 cmplx c2, c3, c4, c5, d2, d3, d4, d5, t2, t3, t4, t5;
139
140 if (ido==1)
141 for (k=0; k<l1; ++k)
142 {
143 PMC (t2,t5,CC(0,1,k),CC(0,4,k))
144 PMC (t3,t4,CC(0,2,k),CC(0,3,k))
145 CH(0,k,0).r=CC(0,0,k).r+t2.r+t3.r;
146 CH(0,k,0).i=CC(0,0,k).i+t2.i+t3.i;
147 c2.r=CC(0,0,k).r+tr11*t2.r+tr12*t3.r;
148 c2.i=CC(0,0,k).i+tr11*t2.i+tr12*t3.i;
149 c3.r=CC(0,0,k).r+tr12*t2.r+tr11*t3.r;
150 c3.i=CC(0,0,k).i+tr12*t2.i+tr11*t3.i;
151 c5.r=ti11*t5.r+ti12*t4.r;
152 c5.i=ti11*t5.i+ti12*t4.i;
153 c4.r=ti12*t5.r-ti11*t4.r;
154 c4.i=ti12*t5.i-ti11*t4.i;
155 CONJFLIPC(c5)
156 PMC(CH(0,k,1),CH(0,k,4),c2,c5)
157 CONJFLIPC(c4)
158 PMC(CH(0,k,2),CH(0,k,3),c3,c4)
159 }
160 else
161 for (k=0; k<l1; ++k)
162 for (i=0; i<ido; ++i)
163 {
164 PMC (t2,t5,CC(i,1,k),CC(i,4,k))
165 PMC (t3,t4,CC(i,2,k),CC(i,3,k))
166 CH(i,k,0).r=CC(i,0,k).r+t2.r+t3.r;
167 CH(i,k,0).i=CC(i,0,k).i+t2.i+t3.i;
168 c2.r=CC(i,0,k).r+tr11*t2.r+tr12*t3.r;
169 c2.i=CC(i,0,k).i+tr11*t2.i+tr12*t3.i;
170 c3.r=CC(i,0,k).r+tr12*t2.r+tr11*t3.r;
171 c3.i=CC(i,0,k).i+tr12*t2.i+tr11*t3.i;
172 c5.r=ti11*t5.r+ti12*t4.r;
173 c5.i=ti11*t5.i+ti12*t4.i;
174 c4.r=ti12*t5.r-ti11*t4.r;
175 c4.i=ti12*t5.i-ti11*t4.i;
176 CONJFLIPC(c5)
177 PMC(d2,d5,c2,c5)
178 CONJFLIPC(c4)
179 PMC(d3,d4,c3,c4)
180 MULPMSIGNC (CH(i,k,1),WA(0,i),d2)
181 MULPMSIGNC (CH(i,k,2),WA(1,i),d3)
182 MULPMSIGNC (CH(i,k,3),WA(2,i),d4)
183 MULPMSIGNC (CH(i,k,4),WA(3,i),d5)
184 }
185 }
186
187 static void X(6)(size_t ido, size_t l1, const cmplx *cc, cmplx *ch,
188 const cmplx *wa)
189 {
190 const size_t cdim=6;
191 static const double taui= PSIGN 0.86602540378443864676;
192 cmplx ta1,ta2,ta3,a0,a1,a2,tb1,tb2,tb3,b0,b1,b2,d1,d2,d3,d4,d5;
193 size_t i, k;
194
195 if (ido==1)
196 for (k=0; k<l1; ++k)
197 {
198 PMC(ta1,ta3,CC(0,2,k),CC(0,4,k))
199 ta2.r = CC(0,0,k).r - .5*ta1.r;
200 ta2.i = CC(0,0,k).i - .5*ta1.i;
201 SCALEC(ta3,taui)
202 ADDC(a0,CC(0,0,k),ta1)
203 CONJFLIPC(ta3)
204 PMC(a1,a2,ta2,ta3)
205 PMC(tb1,tb3,CC(0,5,k),CC(0,1,k))
206 tb2.r = CC(0,3,k).r - .5*tb1.r;
207 tb2.i = CC(0,3,k).i - .5*tb1.i;
208 SCALEC(tb3,taui)
209 ADDC(b0,CC(0,3,k),tb1)
210 CONJFLIPC(tb3)
211 PMC(b1,b2,tb2,tb3)
212 PMC(CH(0,k,0),CH(0,k,3),a0,b0)
213 PMC(CH(0,k,4),CH(0,k,1),a1,b1)
214 PMC(CH(0,k,2),CH(0,k,5),a2,b2)
215 }
216 else
217 for (k=0; k<l1; ++k)
218 for (i=0; i<ido; ++i)
219 {
220 PMC(ta1,ta3,CC(i,2,k),CC(i,4,k))
221 ta2.r = CC(i,0,k).r - .5*ta1.r;
222 ta2.i = CC(i,0,k).i - .5*ta1.i;
223 SCALEC(ta3,taui)
224 ADDC(a0,CC(i,0,k),ta1)
225 CONJFLIPC(ta3)
226 PMC(a1,a2,ta2,ta3)
227 PMC(tb1,tb3,CC(i,5,k),CC(i,1,k))
228 tb2.r = CC(i,3,k).r - .5*tb1.r;
229 tb2.i = CC(i,3,k).i - .5*tb1.i;
230 SCALEC(tb3,taui)
231 ADDC(b0,CC(i,3,k),tb1)
232 CONJFLIPC(tb3)
233 PMC(b1,b2,tb2,tb3)
234 PMC(CH(i,k,0),d3,a0,b0)
235 PMC(d4,d1,a1,b1)
236 PMC(d2,d5,a2,b2)
237 MULPMSIGNC (CH(i,k,1),WA(0,i),d1)
238 MULPMSIGNC (CH(i,k,2),WA(1,i),d2)
239 MULPMSIGNC (CH(i,k,3),WA(2,i),d3)
240 MULPMSIGNC (CH(i,k,4),WA(3,i),d4)
241 MULPMSIGNC (CH(i,k,5),WA(4,i),d5)
242 }
243 }
244
X(g)245 static void X(g)(size_t ido, size_t ip, size_t l1, const cmplx *cc, cmplx *ch,
246 const cmplx *wa)
247 {
248 const size_t cdim=ip;
249 cmplx *tarr=RALLOC(cmplx,2*ip);
250 cmplx *ccl=tarr, *wal=tarr+ip;
251 size_t i,j,k,l,jc,lc;
252 size_t ipph = (ip+1)/2;
253
254 for (i=1; i<ip; ++i)
255 wal[i]=wa[ido*(i-1)];
256 for (k=0; k<l1; ++k)
257 for (i=0; i<ido; ++i)
258 {
259 cmplx s=CC(i,0,k);
260 ccl[0] = CC(i,0,k);
261 for(j=1,jc=ip-1; j<ipph; ++j,--jc)
262 {
263 PMC (ccl[j],ccl[jc],CC(i,j,k),CC(i,jc,k))
264 ADDC (s,s,ccl[j])
265 }
266 CH(i,k,0) = s;
267 for (j=1, jc=ip-1; j<=ipph; ++j,--jc)
268 {
269 cmplx abr=ccl[0], abi={0.,0.};
270 size_t iang=0;
271 for (l=1,lc=ip-1; l<ipph; ++l,--lc)
272 {
273 iang+=j;
274 if (iang>ip) iang-=ip;
275 abr.r += ccl[l ].r*wal[iang].r;
276 abr.i += ccl[l ].i*wal[iang].r;
277 abi.r += ccl[lc].r*wal[iang].i;
278 abi.i += ccl[lc].i*wal[iang].i;
279 }
280 #ifndef BACKWARD
281 { abi.i=-abi.i; abi.r=-abi.r; }
282 #endif
283 CONJFLIPC(abi)
284 PMC(CH(i,k,j),CH(i,k,jc),abr,abi)
285 }
286 }
287
288 DEALLOC(tarr);
289
290 if (ido==1) return;
291
292 for (j=1; j<ip; ++j)
293 for (k=0; k<l1; ++k)
294 {
295 size_t idij=(j-1)*ido+1;
296 for(i=1; i<ido; ++i, ++idij)
297 {
298 cmplx t=CH(i,k,j);
299 MULPMSIGNC (CH(i,k,j),wa[idij],t)
300 }
301 }
302 }
303
304 #undef PSIGN
305 #undef PMSIGNC
306 #undef MULPMSIGNC
307