1 /*
2 * Copyright (c) 2007 - 2015 Joseph Gaeddert
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23 //
24 // fft_dft.c : definitions for regular, slow DFTs
25 //
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30 #include "liquid.internal.h"
31
32 // create FFT plan for regular DFT
33 // _nfft : FFT size
34 // _x : input array [size: _nfft x 1]
35 // _y : output array [size: _nfft x 1]
36 // _dir : fft direction: {LIQUID_FFT_FORWARD, LIQUID_FFT_BACKWARD}
37 // _method : fft method
FFT(_create_plan_dft)38 FFT(plan) FFT(_create_plan_dft)(unsigned int _nfft,
39 TC * _x,
40 TC * _y,
41 int _dir,
42 int _flags)
43 {
44 // allocate plan and initialize all internal arrays to NULL
45 FFT(plan) q = (FFT(plan)) malloc(sizeof(struct FFT(plan_s)));
46
47 q->nfft = _nfft;
48 q->x = _x;
49 q->y = _y;
50 q->flags = _flags;
51 q->type = (_dir == LIQUID_FFT_FORWARD) ? LIQUID_FFT_FORWARD : LIQUID_FFT_BACKWARD;
52 q->direction = (_dir == LIQUID_FFT_FORWARD) ? LIQUID_FFT_FORWARD : LIQUID_FFT_BACKWARD;
53 q->method = LIQUID_FFT_METHOD_DFT;
54
55 q->data.dft.twiddle = NULL;
56 q->data.dft.dotprod = NULL;
57
58 // check size, use specific codelet for small DFTs
59 if (q->nfft == 2) q->execute = FFT(_execute_dft_2);
60 else if (q->nfft == 3) q->execute = FFT(_execute_dft_3);
61 else if (q->nfft == 4) q->execute = FFT(_execute_dft_4);
62 else if (q->nfft == 5) q->execute = FFT(_execute_dft_5);
63 else if (q->nfft == 6) q->execute = FFT(_execute_dft_6);
64 else if (q->nfft == 7) q->execute = FFT(_execute_dft_7);
65 else if (q->nfft == 8) q->execute = FFT(_execute_dft_8);
66 else if (q->nfft == 16) q->execute = FFT(_execute_dft_16);
67 else {
68 q->execute = FFT(_execute_dft);
69
70 // initialize twiddle factors
71 q->data.dft.twiddle = (TC *) malloc(q->nfft * sizeof(TC));
72
73 // create dotprod objects
74 q->data.dft.dotprod = (DOTPROD()*) malloc(q->nfft * sizeof(DOTPROD()));
75
76 // create dotprod objects
77 // twiddles: exp(-j*2*pi*W/n), W=
78 // 0 0 0 0 0...
79 // 0 1 2 3 4...
80 // 0 2 4 6 8...
81 // 0 3 6 9 12...
82 // ...
83 // Note that first row/column is zero, no multiplication necessary.
84 // Create dotprod for first row anyway because it's still faster...
85 unsigned int i;
86 unsigned int k;
87 T d = (q->direction == LIQUID_FFT_FORWARD) ? -1.0 : 1.0;
88 for (i=0; i<q->nfft; i++) {
89 // initialize twiddle factors
90 // NOTE: no need to compute first twiddle because exp(-j*2*pi*0) = 1
91 for (k=1; k<q->nfft; k++)
92 q->data.dft.twiddle[k-1] = cexpf(_Complex_I*d*2*M_PI*(T)(k*i) / (T)(q->nfft));
93
94 // create dotprod object
95 q->data.dft.dotprod[i] = DOTPROD(_create)(q->data.dft.twiddle, q->nfft-1);
96 }
97 }
98
99 return q;
100 }
101
102 // destroy FFT plan
FFT(_destroy_plan_dft)103 void FFT(_destroy_plan_dft)(FFT(plan) _q)
104 {
105 // free twiddle factors
106 if (_q->data.dft.twiddle != NULL)
107 free(_q->data.dft.twiddle);
108
109 // free dotprod objects
110 if (_q->data.dft.dotprod != NULL) {
111 unsigned int i;
112 for (i=0; i<_q->nfft; i++)
113 DOTPROD(_destroy)(_q->data.dft.dotprod[i]);
114
115 // free dotprod array
116 free(_q->data.dft.dotprod);
117 }
118
119 // free main object memory
120 free(_q);
121 }
122
123 // execute DFT (slow but functionally correct)
FFT(_execute_dft)124 void FFT(_execute_dft)(FFT(plan) _q)
125 {
126 unsigned int i;
127 unsigned int nfft = _q->nfft;
128
129 #if 0
130 // DC value is sum of input
131 _q->y[0] = _q->x[0];
132 for (i=1; i<nfft; i++)
133 _q->y[0] += _q->x[i];
134
135 // compute remaining DFT values
136 unsigned int k;
137 for (i=1; i<nfft; i++) {
138 _q->y[i] = _q->x[0];
139 for (k=1; k<nfft; k++) {
140 _q->y[i] += _q->x[k] * _q->data.dft.twiddle[(i*k)%_q->nfft];
141 }
142 }
143 #else
144 // use vector dot products
145 // NOTE: no need to compute first multiplication because exp(-j*2*pi*0) = 1
146 for (i=0; i<nfft; i++) {
147 DOTPROD(_execute)(_q->data.dft.dotprod[i], &_q->x[1], &_q->y[i]);
148 _q->y[i] += _q->x[0];
149 }
150 #endif
151 }
152
153 //
154 // codelets for small DFTs
155 //
156
157 //
FFT(_execute_dft_2)158 void FFT(_execute_dft_2)(FFT(plan) _q)
159 {
160 _q->y[0] = _q->x[0] + _q->x[1];
161 _q->y[1] = _q->x[0] - _q->x[1];
162 }
163
164 //
FFT(_execute_dft_3)165 void FFT(_execute_dft_3)(FFT(plan) _q)
166 {
167 #if 0
168 // NOTE: not as fast as other method, but perhaps useful for
169 // fixed-point algorithm
170 // x = a + jb
171 // y = c + jd
172 // We want to compute both x*y and x*conj(y) with as few
173 // multiplications as possible. If we define
174 // k1 = a*(c+d);
175 // k2 = d*(a+b);
176 // k3 = c*(b-a);
177 // k4 = b*(c+d);
178 // then
179 // x*y = (k1-k2) + j(k1+k3)
180 // x*conj(y) = (k4-k3) + j(k4-k2)
181 T a, b; // c=real(g)=-0.5, d=imag(g)=-sqrt(3)/2
182 T k1, k2, k3, k4;
183
184 // compute both _q->x[1]*g and _q->x[1]*conj(g) with only 4 real multiplications
185 a = crealf(_q->x[1]);
186 b = cimagf(_q->x[1]);
187 //k1 = a*(-0.5f + -0.866025403784439f);
188 k1 = -1.36602540378444f*a;
189 k2 = -0.866025403784439f*( a + b);
190 k3 = -0.5f*( b - a);
191 //k4 = b*(-0.5f + -0.866025403784439f);
192 k4 = -1.36602540378444f*b;
193
194 TC ta1 = (k1-k2) + _Complex_I*(k1+k3); //
195 TC tb1 = (k4-k3) + _Complex_I*(k4-k2); //
196
197 // compute both _q->x[2]*g and _q->x[2]*conj(g) with only 4 real multiplications
198 a = crealf(_q->x[2]);
199 b = cimagf(_q->x[2]);
200 //k1 = a*(-0.5f + -0.866025403784439f);
201 k1 = -1.36602540378444f*a;
202 k2 = -0.866025403784439f*( a + b);
203 k3 = -0.5f*( b - a);
204 //k4 = b*(-0.5f + -0.866025403784439f);
205 k4 = -1.36602540378444f*b;
206
207 TC ta2 = (k1-k2) + _Complex_I*(k1+k3); //
208 TC tb2 = (k4-k3) + _Complex_I*(k4-k2); //
209
210 // set return values
211 _q->y[0] = _q->x[0] + _q->x[1] + _q->x[2];
212 if (_q->direction == LIQUID_FFT_FORWARD) {
213 _q->y[1] = _q->x[0] + ta1 + tb2;
214 _q->y[2] = _q->x[0] + tb1 + ta2;
215 } else {
216 _q->y[1] = _q->x[0] + tb1 + ta2;
217 _q->y[2] = _q->x[0] + ta1 + tb2;
218 }
219 #else
220 TC g = -0.5f - _Complex_I*0.866025403784439; // sqrt(3)/2
221
222 _q->y[0] = _q->x[0] + _q->x[1] + _q->x[2];
223 TC ta = _q->x[0] + _q->x[1]*g + _q->x[2]*conjf(g);
224 TC tb = _q->x[0] + _q->x[1]*conjf(g) + _q->x[2]*g;
225
226 // set return values
227 if (_q->direction == LIQUID_FFT_FORWARD) {
228 _q->y[1] = ta;
229 _q->y[2] = tb;
230 } else {
231 _q->y[1] = tb;
232 _q->y[2] = ta;
233 }
234 #endif
235 }
236
237 //
FFT(_execute_dft_4)238 void FFT(_execute_dft_4)(FFT(plan) _q)
239 {
240 TC yp;
241 TC * x = _q->x;
242 TC * y = _q->y;
243
244 // index reversal
245 y[0] = x[0];
246 y[1] = x[2];
247 y[2] = x[1];
248 y[3] = x[3];
249
250 // k0 = 0, k1=1
251 yp = y[1];
252 y[1] = y[0] - yp;
253 y[0] = y[0] + yp;
254
255 // k0 = 2, k1=3
256 yp = y[3];
257 y[3] = y[2] - yp;
258 y[2] = y[2] + yp;
259
260 // k0 = 0, k1=2
261 yp = y[2];
262 y[2] = y[0] - yp;
263 y[0] = y[0] + yp;
264
265 // k0 = 1, k1=3
266 yp = cimagf(y[3]) - _Complex_I*crealf(y[3]);
267 if (_q->direction == LIQUID_FFT_BACKWARD)
268 yp = -yp;
269 y[3] = y[1] - yp;
270 y[1] = y[1] + yp;
271 }
272
273 //
FFT(_execute_dft_5)274 void FFT(_execute_dft_5)(FFT(plan) _q)
275 {
276 TC * x = _q->x;
277 TC * y = _q->y;
278
279 // DC value is sum of inputs
280 y[0] = x[0] + x[1] + x[2] + x[3] + x[4];
281
282 // exp(-j*2*pi*1/5)
283 TC g0 = 0.309016994374947 - 0.951056516295154*_Complex_I;
284
285 // exp(-j*2*pi*2/5)
286 TC g1 = -0.809016994374947 - 0.587785252292473*_Complex_I;
287
288 if (_q->direction == LIQUID_FFT_BACKWARD) {
289 g0 = conjf(g0);
290 g1 = conjf(g1);
291 }
292 TC g0_conj = conjf(g0);
293 TC g1_conj = conjf(g1);
294
295 y[1] = x[0] + x[1]*g0 + x[2]*g1 + x[3]*g1_conj + x[4]*g0_conj;
296 y[2] = x[0] + x[1]*g1 + x[2]*g0_conj + x[3]*g0 + x[4]*g1_conj;
297 y[3] = x[0] + x[1]*g1_conj + x[2]*g0 + x[3]*g0_conj + x[4]*g1;
298 y[4] = x[0] + x[1]*g0_conj + x[2]*g1_conj + x[3]*g1 + x[4]*g0;
299 }
300
301 //
FFT(_execute_dft_6)302 void FFT(_execute_dft_6)(FFT(plan) _q)
303 {
304 TC * x = _q->x;
305 TC * y = _q->y;
306
307 // DC value is sum of inputs
308 y[0] = x[0] + x[1] + x[2] + x[3] + x[4] + x[5];
309
310 // exp(-j*2*pi*1/6) = 1/2 - j*sqrt(3)/2
311 TC g = 0.5 - 0.866025403784439*_Complex_I;
312
313 TC g1, g2, g4, g5;
314
315 if (_q->direction == LIQUID_FFT_FORWARD) {
316 g1 = g; // exp(-j*2*pi*1/6)
317 g2 = -conjf(g); // exp(-j*2*pi*2/6)
318 g4 = -g; // exp(-j*2*pi*4/6)
319 g5 = conjf(g); // exp(-j*2*pi*5/6)
320 } else {
321 g1 = conjf(g); // exp( j*2*pi*1/6)
322 g2 = -g; // exp( j*2*pi*2/6)
323 g4 = -conjf(g); // exp( j*2*pi*4/6)
324 g5 = g; // exp( j*2*pi*5/6)
325 }
326
327 y[1] = x[0] + x[1]*g1 + x[2]*g2 - x[3] + x[4]*g4 + x[5]*g5;
328 y[2] = x[0] + x[1]*g2 + x[2]*g4 + x[3] + x[4]*g2 + x[5]*g4;
329 y[3] = x[0] - x[1] + x[2] - x[3] + x[4] - x[5];
330 y[4] = x[0] + x[1]*g4 + x[2]*g2 + x[3] + x[4]*g4 + x[5]*g2;
331 y[5] = x[0] + x[1]*g5 + x[2]*g4 - x[3] + x[4]*g2 + x[5]*g1;
332 }
333
334 //
FFT(_execute_dft_7)335 void FFT(_execute_dft_7)(FFT(plan) _q)
336 {
337 TC * x = _q->x;
338 TC * y = _q->y;
339
340 // DC value is sum of inputs
341 y[0] = x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6];
342
343 // initialize twiddle factors
344 TC g1 = 0.623489801858734 - 0.781831482468030 * _Complex_I; // exp(-j*2*pi*1/7)
345 TC g2 = -0.222520933956314 - 0.974927912181824 * _Complex_I; // exp(-j*2*pi*2/7)
346 TC g3 = -0.900968867902419 - 0.433883739117558 * _Complex_I; // exp(-j*2*pi*3/7)
347
348 if (_q->direction == LIQUID_FFT_FORWARD) {
349 } else {
350 g1 = conjf(g1); // exp(+j*2*pi*1/7)
351 g2 = conjf(g2); // exp(+j*2*pi*2/7)
352 g3 = conjf(g3); // exp(+j*2*pi*3/7)
353 }
354
355 TC g4 = conjf(g3);
356 TC g5 = conjf(g2);
357 TC g6 = conjf(g1);
358
359 y[1] = x[0] + x[1]*g1 + x[2]*g2 + x[3]*g3 + x[4]*g4 + x[5]*g5 + x[6]*g6;
360 y[2] = x[0] + x[1]*g2 + x[2]*g4 + x[3]*g6 + x[4]*g1 + x[5]*g3 + x[6]*g5;
361 y[3] = x[0] + x[1]*g3 + x[2]*g6 + x[3]*g2 + x[4]*g5 + x[5]*g1 + x[6]*g4;
362 y[4] = x[0] + x[1]*g4 + x[2]*g1 + x[3]*g5 + x[4]*g2 + x[5]*g6 + x[6]*g3;
363 y[5] = x[0] + x[1]*g5 + x[2]*g3 + x[3]*g1 + x[4]*g6 + x[5]*g4 + x[6]*g2;
364 y[6] = x[0] + x[1]*g6 + x[2]*g5 + x[3]*g4 + x[4]*g3 + x[5]*g2 + x[6]*g1;
365 }
366
367 //
FFT(_execute_dft_8)368 void FFT(_execute_dft_8)(FFT(plan) _q)
369 {
370 TC yp;
371 TC * x = _q->x;
372 TC * y = _q->y;
373
374 // fft or ifft?
375 int fft = _q->direction == LIQUID_FFT_FORWARD ? 1 : 0;
376
377 // index reversal
378 y[0] = x[0];
379 y[1] = x[4];
380 y[2] = x[2];
381 y[3] = x[6];
382 y[4] = x[1];
383 y[5] = x[5];
384 y[6] = x[3];
385 y[7] = x[7];
386
387 // i=0
388 yp = y[1]; y[1] = y[0]-yp; y[0] += yp;
389 yp = y[3]; y[3] = y[2]-yp; y[2] += yp;
390 yp = y[5]; y[5] = y[4]-yp; y[4] += yp;
391 yp = y[7]; y[7] = y[6]-yp; y[6] += yp;
392
393
394 // i=1
395 yp = y[2]; y[2] = y[0]-yp; y[0] += yp;
396 yp = y[6]; y[6] = y[4]-yp; y[4] += yp;
397
398 if (fft) yp = cimagf(y[3]) - crealf(y[3])*_Complex_I;
399 else yp = -cimagf(y[3]) + crealf(y[3])*_Complex_I;
400 y[3] = y[1]-yp;
401 y[1] += yp;
402
403 if (fft) yp = cimagf(y[7]) - crealf(y[7])*_Complex_I;
404 else yp = -cimagf(y[7]) + crealf(y[7])*_Complex_I;
405 y[7] = y[5]-yp;
406 y[5] += yp;
407
408
409 // i=2
410 yp = y[4]; y[4] = y[0]-yp; y[0] += yp;
411
412 if (fft) yp = y[5]*(M_SQRT1_2 - M_SQRT1_2*_Complex_I);
413 else yp = y[5]*(M_SQRT1_2 + M_SQRT1_2*_Complex_I);
414 y[5] = y[1]-yp;
415 y[1] += yp;
416
417 if (fft) yp = cimagf(y[6]) - crealf(y[6])*_Complex_I;
418 else yp = -cimagf(y[6]) + crealf(y[6])*_Complex_I;
419 y[6] = y[2]-yp;
420 y[2] += yp;
421
422 if (fft) yp = y[7]*(-M_SQRT1_2 - M_SQRT1_2*_Complex_I);
423 else yp = y[7]*(-M_SQRT1_2 + M_SQRT1_2*_Complex_I);
424 y[7] = y[3]-yp;
425 y[3] += yp;
426 }
427
428
429 //
FFT(_execute_dft_16)430 void FFT(_execute_dft_16)(FFT(plan) _q)
431 {
432 TC yp;
433 TC * x = _q->x;
434 TC * y = _q->y;
435
436 // fft or ifft?
437 int fft = _q->direction == LIQUID_FFT_FORWARD ? 1 : 0;
438
439 // index reversal
440 y[ 0] = x[ 0];
441 y[ 1] = x[ 8];
442 y[ 2] = x[ 4];
443 y[ 3] = x[12];
444 y[ 4] = x[ 2];
445 y[ 5] = x[10];
446 y[ 6] = x[ 6];
447 y[ 7] = x[14];
448 y[ 8] = x[ 1];
449 y[ 9] = x[ 9];
450 y[10] = x[ 5];
451 y[11] = x[13];
452 y[12] = x[ 3];
453 y[13] = x[11];
454 y[14] = x[ 7];
455 y[15] = x[15];
456
457 // i=0
458 yp = y[ 1]; y[ 1] = y[ 0] - yp; y[ 0] += yp;
459 yp = y[ 3]; y[ 3] = y[ 2] - yp; y[ 2] += yp;
460 yp = y[ 5]; y[ 5] = y[ 4] - yp; y[ 4] += yp;
461 yp = y[ 7]; y[ 7] = y[ 6] - yp; y[ 6] += yp;
462 yp = y[ 9]; y[ 9] = y[ 8] - yp; y[ 8] += yp;
463 yp = y[11]; y[11] = y[10] - yp; y[10] += yp;
464 yp = y[13]; y[13] = y[12] - yp; y[12] += yp;
465 yp = y[15]; y[15] = y[14] - yp; y[14] += yp;
466
467 // i=1
468 yp = y[ 2]; y[ 2] = y[ 0] - yp; y[ 0] += yp;
469 yp = y[ 6]; y[ 6] = y[ 4] - yp; y[ 4] += yp;
470 yp = y[10]; y[10] = y[ 8] - yp; y[ 8] += yp;
471 yp = y[14]; y[14] = y[12] - yp; y[12] += yp;
472 if (fft) {
473 yp = -y[ 3]*_Complex_I; y[ 3] = y[ 1] - yp; y[ 1] += yp;
474 yp = -y[ 7]*_Complex_I; y[ 7] = y[ 5] - yp; y[ 5] += yp;
475 yp = -y[11]*_Complex_I; y[11] = y[ 9] - yp; y[ 9] += yp;
476 yp = -y[15]*_Complex_I; y[15] = y[13] - yp; y[13] += yp;
477 } else {
478 yp = y[ 3]*_Complex_I; y[ 3] = y[ 1] - yp; y[ 1] += yp;
479 yp = y[ 7]*_Complex_I; y[ 7] = y[ 5] - yp; y[ 5] += yp;
480 yp = y[11]*_Complex_I; y[11] = y[ 9] - yp; y[ 9] += yp;
481 yp = y[15]*_Complex_I; y[15] = y[13] - yp; y[13] += yp;
482 }
483
484 // i=2
485 yp = y[ 4]; y[ 4] = y[ 0] - yp; y[ 0] += yp;
486 yp = y[12]; y[12] = y[ 8] - yp; y[ 8] += yp;
487 if (fft) {
488 yp = y[ 5]*( 0.70710677 + _Complex_I* -0.70710677); y[ 5] = y[ 1] - yp; y[ 1] += yp;
489 yp = y[13]*( 0.70710677 + _Complex_I* -0.70710677); y[13] = y[ 9] - yp; y[ 9] += yp;
490 yp = -y[ 6]*_Complex_I; y[ 6] = y[ 2] - yp; y[ 2] += yp;
491 yp = -y[14]*_Complex_I; y[14] = y[10] - yp; y[10] += yp;
492 yp = y[ 7]*( -0.70710677 + _Complex_I* -0.70710677); y[ 7] = y[ 3] - yp; y[ 3] += yp;
493 yp = y[15]*( -0.70710677 + _Complex_I* -0.70710677); y[15] = y[11] - yp; y[11] += yp;
494 } else {
495 yp = y[ 5]*( 0.70710677 - _Complex_I* -0.70710677); y[ 5] = y[ 1] - yp; y[ 1] += yp;
496 yp = y[13]*( 0.70710677 - _Complex_I* -0.70710677); y[13] = y[ 9] - yp; y[ 9] += yp;
497 yp = y[ 6]*_Complex_I; y[ 6] = y[ 2] - yp; y[ 2] += yp;
498 yp = y[14]*_Complex_I; y[14] = y[10] - yp; y[10] += yp;
499 yp = y[ 7]*( -0.70710677 - _Complex_I* -0.70710677); y[ 7] = y[ 3] - yp; y[ 3] += yp;
500 yp = y[15]*( -0.70710677 - _Complex_I* -0.70710677); y[15] = y[11] - yp; y[11] += yp;
501 }
502
503 // i=3
504 yp = y[ 8]; y[ 8] = y[ 0] - yp; y[ 0] += yp;
505 if (fft) {
506 yp = y[ 9]*( 0.92387950 + _Complex_I* -0.38268346); y[ 9] = y[ 1] - yp; y[ 1] += yp;
507 yp = y[10]*( 0.70710677 + _Complex_I* -0.70710677); y[10] = y[ 2] - yp; y[ 2] += yp;
508 yp = y[11]*( 0.38268343 + _Complex_I* -0.92387950); y[11] = y[ 3] - yp; y[ 3] += yp;
509 yp = -y[12]*_Complex_I; y[12] = y[ 4] - yp; y[ 4] += yp;
510 yp = y[13]*( -0.38268340 + _Complex_I* -0.92387956); y[13] = y[ 5] - yp; y[ 5] += yp;
511 yp = y[14]*( -0.70710677 + _Complex_I* -0.70710677); y[14] = y[ 6] - yp; y[ 6] += yp;
512 yp = y[15]*( -0.92387950 + _Complex_I* -0.38268349); y[15] = y[ 7] - yp; y[ 7] += yp;
513 } else {
514 yp = y[ 9]*( 0.92387950 - _Complex_I* -0.38268346); y[ 9] = y[ 1] - yp; y[ 1] += yp;
515 yp = y[10]*( 0.70710677 - _Complex_I* -0.70710677); y[10] = y[ 2] - yp; y[ 2] += yp;
516 yp = y[11]*( 0.38268343 - _Complex_I* -0.92387950); y[11] = y[ 3] - yp; y[ 3] += yp;
517 yp = y[12]*_Complex_I; y[12] = y[ 4] - yp; y[ 4] += yp;
518 yp = y[13]*( -0.38268340 - _Complex_I* -0.92387956); y[13] = y[ 5] - yp; y[ 5] += yp;
519 yp = y[14]*( -0.70710677 - _Complex_I* -0.70710677); y[14] = y[ 6] - yp; y[ 6] += yp;
520 yp = y[15]*( -0.92387950 - _Complex_I* -0.38268349); y[15] = y[ 7] - yp; y[ 7] += yp;
521 }
522 }
523
524