1 /********************************************************************
2 * *
3 * THIS FILE IS PART OF THE OggSQUISH SOFTWARE CODEC SOURCE CODE. *
4 * *
5 ********************************************************************
6
7 function: Fast discrete Fourier and cosine transforms and inverses
8 author: Monty <xiphmont@mit.edu>
9 modifications by: Monty
10 last modification date: Jul 1 1996
11
12 djmw 20030630 Adapted for praat (replaced 'int' declarations with 'long').
13 djmw 20040511 Made all local variables type double to increase numerical precision.
14 djmw 20171003 Replaced `long` declarations with `integer`).
15
16 ********************************************************************/
17
18 /* These Fourier routines were originally based on the Fourier routines of
19 the same names from the NETLIB bihar and fftpack fortran libraries
20 developed by Paul N. Swarztrauber at the National Center for Atmospheric
21 Research in Boulder, CO USA. They have been reimplemented in C and
22 optimized in a few ways for OggSquish. */
23
24 /* As the original fortran libraries are public domain, the C Fourier
25 routines in this file are hereby released to the public domain as well.
26 The C routines here produce output exactly equivalent to the original
27 fortran routines. Of particular interest are the facts that (like the
28 original fortran), these routines can work on arbitrary length vectors
29 that need not be powers of two in length. */
30
31 #include "melder.h" /* for integer */
32
drfti1(integer n,FFT_DATA_TYPE * wa,integer * ifac)33 static void drfti1 (integer n, FFT_DATA_TYPE * wa, integer *ifac)
34 {
35 static constexpr integer ntryh [4] = { 4, 2, 3, 5 };
36 static constexpr double tpi = 6.28318530717958647692528676655900577;
37 integer ntry = 0, j = -1;
38 integer nl = n;
39 integer nf = 0;
40
41 L101:
42 j++;
43 if (j < 4)
44 ntry = ntryh [j];
45 else
46 ntry += 2;
47
48 L104:
49 const integer nq = nl / ntry;
50 const integer nr = nl - ntry * nq;
51 if (nr != 0)
52 goto L101;
53
54 nf++;
55 ifac [nf + 1] = ntry;
56 nl = nq;
57 if (ntry != 2)
58 goto L107;
59 if (nf == 1)
60 goto L107;
61
62 for (integer i = 1; i < nf; i++)
63 {
64 const integer ib = nf - i + 1;
65 ifac [ib + 1] = ifac [ib];
66 }
67 ifac [2] = 2;
68
69 L107:
70 if (nl != 1)
71 goto L104;
72 ifac [0] = n;
73 ifac [1] = nf;
74 const double argh = tpi / n;
75 integer is = 0;
76 const integer nfm1 = nf - 1;
77 integer l1 = 1;
78
79 if (nfm1 == 0)
80 return;
81 for (integer k1 = 0; k1 < nfm1; k1++)
82 {
83 const integer ip = ifac [k1 + 2];
84 integer ld = 0;
85 const integer l2 = l1 * ip;
86 const integer ido = n / l2;
87 const integer ipm = ip - 1;
88
89 for (j = 0; j < ipm; j++)
90 {
91 ld += l1;
92 integer i = is;
93 const double argld = (double) ld *argh;
94
95 double fi = 0.0;
96 for (integer ii = 2; ii < ido; ii += 2)
97 {
98 fi += 1.0;
99 const double arg = fi * argld;
100 wa [i++] = cos (arg);
101 wa [i++] = sin (arg);
102 }
103 is += ido;
104 }
105 l1 = l2;
106 }
107 }
108
NUMrffti(integer n,FFT_DATA_TYPE * wsave,integer * ifac)109 static void NUMrffti (integer n, FFT_DATA_TYPE * wsave, integer *ifac)
110 {
111 if (n == 1)
112 return;
113 drfti1 (n, wsave + n, ifac);
114 }
115
116 /* void NUMcosqi(integer n, FFT_DATA_TYPE *wsave, integer *ifac){ static
117 double pih = 1.57079632679489661923132169163975; static integer k;
118 static double fk, dt;
119
120 dt=pih/n; fk=0.; for(k=0;k<n;k++){ fk+=1.; wsave [k] = cos(fk*dt); }
121
122 NUMrffti(n, wsave+n,ifac); } */
123
dradf2(integer ido,integer l1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa1)124 static void dradf2 (integer ido, integer l1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa1)
125 {
126 integer t1 = 0;
127 integer t2, t0 = (t2 = l1 * ido);
128 integer t3 = ido << 1;
129 for (integer k = 0; k < l1; k++)
130 {
131 ch [t1 << 1] = cc [t1] + cc [t2];
132 ch [(t1 << 1) + t3 - 1] = cc [t1] - cc [t2];
133 t1 += ido;
134 t2 += ido;
135 }
136
137 if (ido < 2)
138 return;
139 if (ido == 2)
140 goto L105;
141
142 t1 = 0;
143 t2 = t0;
144 for (integer k = 0; k < l1; k++)
145 {
146 t3 = t2;
147 integer t4 = (t1 << 1) + (ido << 1);
148 integer t5 = t1;
149 integer t6 = t1 + t1;
150 for (integer i = 2; i < ido; i += 2)
151 {
152 t3 += 2;
153 t4 -= 2;
154 t5 += 2;
155 t6 += 2;
156 const double tr2 = wa1 [i - 2] * cc [t3 - 1] + wa1 [i - 1] * cc [t3];
157 const double ti2 = wa1 [i - 2] * cc [t3] - wa1 [i - 1] * cc [t3 - 1];
158 ch [t6] = cc [t5] + ti2;
159 ch [t4] = ti2 - cc [t5];
160 ch [t6 - 1] = cc [t5 - 1] + tr2;
161 ch [t4 - 1] = cc [t5 - 1] - tr2;
162 }
163 t1 += ido;
164 t2 += ido;
165 }
166
167 if (ido % 2 == 1)
168 return;
169
170 L105:
171 t3 = (t2 = (t1 = ido) - 1);
172 t2 += t0;
173 for (integer k = 0; k < l1; k++)
174 {
175 ch [t1] = -cc [t2];
176 ch [t1 - 1] = cc [t3];
177 t1 += ido << 1;
178 t2 += ido;
179 t3 += ido;
180 }
181 }
182
dradf4(integer ido,integer l1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa1,FFT_DATA_TYPE * wa2,FFT_DATA_TYPE * wa3)183 static void dradf4 (integer ido, integer l1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa1,
184 FFT_DATA_TYPE * wa2, FFT_DATA_TYPE * wa3)
185 {
186 static constexpr double hsqt2 = .70710678118654752440084436210485;
187 integer t5, t6;
188
189 integer t0 = l1 * ido;
190
191 integer t1 = t0;
192 integer t4 = t1 << 1;
193 integer t2 = t1 + (t1 << 1);
194 integer t3 = 0;
195
196 for (integer k = 0; k < l1; k++)
197 {
198 const double tr1 = cc [t1] + cc [t2];
199 const double tr2 = cc [t3] + cc [t4];
200 ch [t5 = t3 << 2] = tr1 + tr2;
201 ch [(ido << 2) + t5 - 1] = tr2 - tr1;
202 ch [(t5 += (ido << 1)) - 1] = cc [t3] - cc [t4];
203 ch [t5] = cc [t2] - cc [t1];
204
205 t1 += ido;
206 t2 += ido;
207 t3 += ido;
208 t4 += ido;
209 }
210
211 if (ido < 2)
212 return;
213 if (ido == 2)
214 goto L105;
215
216 t1 = 0;
217 for (integer k = 0; k < l1; k++)
218 {
219 t2 = t1;
220 t4 = t1 << 2;
221 t5 = (t6 = ido << 1) + t4;
222 for (integer i = 2; i < ido; i += 2)
223 {
224 t3 = (t2 += 2);
225 t4 += 2;
226 t5 -= 2;
227
228 t3 += t0;
229 const double cr2 = wa1 [i - 2] * cc [t3 - 1] + wa1 [i - 1] * cc [t3];
230 const double ci2 = wa1 [i - 2] * cc [t3] - wa1 [i - 1] * cc [t3 - 1];
231 t3 += t0;
232 const double cr3 = wa2 [i - 2] * cc [t3 - 1] + wa2 [i - 1] * cc [t3];
233 const double ci3 = wa2 [i - 2] * cc [t3] - wa2 [i - 1] * cc [t3 - 1];
234 t3 += t0;
235 const double cr4 = wa3 [i - 2] * cc [t3 - 1] + wa3 [i - 1] * cc [t3];
236 const double ci4 = wa3 [i - 2] * cc [t3] - wa3 [i - 1] * cc [t3 - 1];
237
238 const double tr1 = cr2 + cr4;
239 const double tr4 = cr4 - cr2;
240 const double ti1 = ci2 + ci4;
241 const double ti4 = ci2 - ci4;
242 const double ti2 = cc [t2] + ci3;
243 const double ti3 = cc [t2] - ci3;
244 const double tr2 = cc [t2 - 1] + cr3;
245 const double tr3 = cc [t2 - 1] - cr3;
246
247 ch [t4 - 1] = tr1 + tr2;
248 ch [t4] = ti1 + ti2;
249
250 ch [t5 - 1] = tr3 - ti4;
251 ch [t5] = tr4 - ti3;
252
253 ch [t4 + t6 - 1] = ti4 + tr3;
254 ch [t4 + t6] = tr4 + ti3;
255
256 ch [t5 + t6 - 1] = tr2 - tr1;
257 ch [t5 + t6] = ti1 - ti2;
258 }
259 t1 += ido;
260 }
261 if (ido % 2 == 1)
262 return;
263
264 L105:
265
266 t2 = (t1 = t0 + ido - 1) + (t0 << 1);
267 t3 = ido << 2;
268 t4 = ido;
269 t5 = ido << 1;
270 t6 = ido;
271
272 for (integer k = 0; k < l1; k++)
273 {
274 const double ti1 = -hsqt2 * (cc [t1] + cc [t2]);
275 const double tr1 = hsqt2 * (cc [t1] - cc [t2]);
276 ch [t4 - 1] = tr1 + cc [t6 - 1];
277 ch [t4 + t5 - 1] = cc [t6 - 1] - tr1;
278 ch [t4] = ti1 - cc [t1 + t0];
279 ch [t4 + t5] = ti1 + cc [t1 + t0];
280 t1 += ido;
281 t2 += ido;
282 t4 += t3;
283 t6 += ido;
284 }
285 }
286
dradfg(integer ido,integer ip,integer l1,integer idl1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * c1,FFT_DATA_TYPE * c2,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * ch2,FFT_DATA_TYPE * wa)287 static void dradfg (integer ido, integer ip, integer l1, integer idl1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * c1,
288 FFT_DATA_TYPE * c2, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * ch2, FFT_DATA_TYPE * wa)
289 {
290
291 static constexpr double tpi = 6.28318530717958647692528676655900577;
292 integer is;
293 integer t1, t2, t3, t4, t5, t6, t7, t8, t9;
294
295 const double arg = tpi / (double) ip;
296 const double dcp = cos (arg);
297 const double dsp = sin (arg);
298 const integer ipph = (ip + 1) >> 1;
299 const integer ipp2 = ip;
300 const integer idp2 = ido;
301 const integer nbd = (ido - 1) >> 1;
302 const integer t0 = l1 * ido;
303 const integer t10 = ip * ido;
304
305 if (ido == 1)
306 goto L119;
307 for (integer ik = 0; ik < idl1; ik++)
308 ch2 [ik] = c2 [ik];
309
310 t1 = 0;
311 for (integer j = 1; j < ip; j++)
312 {
313 t1 += t0;
314 t2 = t1;
315 for (integer k = 0; k < l1; k++)
316 {
317 ch [t2] = c1 [t2];
318 t2 += ido;
319 }
320 }
321
322 is = -ido;
323 t1 = 0;
324 if (nbd > l1)
325 {
326 for (integer j = 1; j < ip; j++)
327 {
328 t1 += t0;
329 is += ido;
330 t2 = -ido + t1;
331 for (integer k = 0; k < l1; k++)
332 {
333 integer idij = is - 1;
334 t2 += ido;
335 t3 = t2;
336 for (integer i = 2; i < ido; i += 2)
337 {
338 idij += 2;
339 t3 += 2;
340 ch [t3 - 1] = wa [idij - 1] * c1 [t3 - 1] + wa [idij] * c1 [t3];
341 ch [t3] = wa [idij - 1] * c1 [t3] - wa [idij] * c1 [t3 - 1];
342 }
343 }
344 }
345 }
346 else
347 {
348
349 for (integer j = 1; j < ip; j++)
350 {
351 is += ido;
352 integer idij = is - 1;
353 t1 += t0;
354 t2 = t1;
355 for (integer i = 2; i < ido; i += 2)
356 {
357 idij += 2;
358 t2 += 2;
359 t3 = t2;
360 for (integer k = 0; k < l1; k++)
361 {
362 ch [t3 - 1] = wa [idij - 1] * c1 [t3 - 1] + wa [idij] * c1 [t3];
363 ch [t3] = wa [idij - 1] * c1 [t3] - wa [idij] * c1 [t3 - 1];
364 t3 += ido;
365 }
366 }
367 }
368 }
369
370 t1 = 0;
371 t2 = ipp2 * t0;
372 if (nbd < l1)
373 {
374 for (integer j = 1; j < ipph; j++)
375 {
376 t1 += t0;
377 t2 -= t0;
378 t3 = t1;
379 t4 = t2;
380 for (integer i = 2; i < ido; i += 2)
381 {
382 t3 += 2;
383 t4 += 2;
384 t5 = t3 - ido;
385 t6 = t4 - ido;
386 for (integer k = 0; k < l1; k++)
387 {
388 t5 += ido;
389 t6 += ido;
390 c1 [t5 - 1] = ch [t5 - 1] + ch [t6 - 1];
391 c1 [t6 - 1] = ch [t5] - ch [t6];
392 c1 [t5] = ch [t5] + ch [t6];
393 c1 [t6] = ch [t6 - 1] - ch [t5 - 1];
394 }
395 }
396 }
397 }
398 else
399 {
400 for (integer j = 1; j < ipph; j++)
401 {
402 t1 += t0;
403 t2 -= t0;
404 t3 = t1;
405 t4 = t2;
406 for (integer k = 0; k < l1; k++)
407 {
408 t5 = t3;
409 t6 = t4;
410 for (integer i = 2; i < ido; i += 2)
411 {
412 t5 += 2;
413 t6 += 2;
414 c1 [t5 - 1] = ch [t5 - 1] + ch [t6 - 1];
415 c1 [t6 - 1] = ch [t5] - ch [t6];
416 c1 [t5] = ch [t5] + ch [t6];
417 c1 [t6] = ch [t6 - 1] - ch [t5 - 1];
418 }
419 t3 += ido;
420 t4 += ido;
421 }
422 }
423 }
424
425 L119:
426 for (integer ik = 0; ik < idl1; ik++)
427 c2 [ik] = ch2 [ik];
428
429 t1 = 0;
430 t2 = ipp2 * idl1;
431 for (integer j = 1; j < ipph; j++)
432 {
433 t1 += t0;
434 t2 -= t0;
435 t3 = t1 - ido;
436 t4 = t2 - ido;
437 for (integer k = 0; k < l1; k++)
438 {
439 t3 += ido;
440 t4 += ido;
441 c1 [t3] = ch [t3] + ch [t4];
442 c1 [t4] = ch [t4] - ch [t3];
443 }
444 }
445
446 double ar1 = 1.;
447 double ai1 = 0.;
448 t1 = 0;
449 t2 = ipp2 * idl1;
450 t3 = (ip - 1) * idl1;
451 for (integer l = 1; l < ipph; l++)
452 {
453 t1 += idl1;
454 t2 -= idl1;
455 const double ar1h = dcp * ar1 - dsp * ai1;
456 ai1 = dcp * ai1 + dsp * ar1;
457 ar1 = ar1h;
458 t4 = t1;
459 t5 = t2;
460 t6 = t3;
461 t7 = idl1;
462
463 for (integer ik = 0; ik < idl1; ik++)
464 {
465 ch2 [t4++] = c2 [ik] + ar1 * c2 [t7++];
466 ch2 [t5++] = ai1 * c2 [t6++];
467 }
468
469 const double dc2 = ar1;
470 const double ds2 = ai1;
471 double ar2 = ar1;
472 double ai2 = ai1;
473
474 t4 = idl1;
475 t5 = (ipp2 - 1) * idl1;
476 for (integer j = 2; j < ipph; j++)
477 {
478 t4 += idl1;
479 t5 -= idl1;
480
481 const double ar2h = dc2 * ar2 - ds2 * ai2;
482 ai2 = dc2 * ai2 + ds2 * ar2;
483 ar2 = ar2h;
484
485 t6 = t1;
486 t7 = t2;
487 t8 = t4;
488 t9 = t5;
489 for (integer ik = 0; ik < idl1; ik++)
490 {
491 ch2 [t6++] += ar2 * c2 [t8++];
492 ch2 [t7++] += ai2 * c2 [t9++];
493 }
494 }
495 }
496
497 t1 = 0;
498 for (integer j = 1; j < ipph; j++)
499 {
500 t1 += idl1;
501 t2 = t1;
502 for (integer ik = 0; ik < idl1; ik++)
503 ch2 [ik] += c2 [t2++];
504 }
505
506 if (ido < l1)
507 goto L132;
508
509 t1 = 0;
510 t2 = 0;
511 for (integer k = 0; k < l1; k++)
512 {
513 t3 = t1;
514 t4 = t2;
515 for (integer i = 0; i < ido; i++)
516 cc [t4++] = ch [t3++];
517 t1 += ido;
518 t2 += t10;
519 }
520
521 goto L135;
522
523 L132:
524 for (integer i = 0; i < ido; i++)
525 {
526 t1 = i;
527 t2 = i;
528 for (integer k = 0; k < l1; k++)
529 {
530 cc [t2] = ch [t1];
531 t1 += ido;
532 t2 += t10;
533 }
534 }
535
536 L135:
537 t1 = 0;
538 t2 = ido << 1;
539 t3 = 0;
540 t4 = ipp2 * t0;
541 for (integer j = 1; j < ipph; j++)
542 {
543
544 t1 += t2;
545 t3 += t0;
546 t4 -= t0;
547
548 t5 = t1;
549 t6 = t3;
550 t7 = t4;
551
552 for (integer k = 0; k < l1; k++)
553 {
554 cc [t5 - 1] = ch [t6];
555 cc [t5] = ch [t7];
556 t5 += t10;
557 t6 += ido;
558 t7 += ido;
559 }
560 }
561
562 if (ido == 1)
563 return;
564 if (nbd < l1)
565 goto L141;
566
567 t1 = -ido;
568 t3 = 0;
569 t4 = 0;
570 t5 = ipp2 * t0;
571 for (integer j = 1; j < ipph; j++)
572 {
573 t1 += t2;
574 t3 += t2;
575 t4 += t0;
576 t5 -= t0;
577 t6 = t1;
578 t7 = t3;
579 t8 = t4;
580 t9 = t5;
581 for (integer k = 0; k < l1; k++)
582 {
583 for (integer i = 2; i < ido; i += 2)
584 {
585 const integer ic = idp2 - i;
586 cc [i + t7 - 1] = ch [i + t8 - 1] + ch [i + t9 - 1];
587 cc [ic + t6 - 1] = ch [i + t8 - 1] - ch [i + t9 - 1];
588 cc [i + t7] = ch [i + t8] + ch [i + t9];
589 cc [ic + t6] = ch [i + t9] - ch [i + t8];
590 }
591 t6 += t10;
592 t7 += t10;
593 t8 += ido;
594 t9 += ido;
595 }
596 }
597 return;
598
599 L141:
600
601 t1 = -ido;
602 t3 = 0;
603 t4 = 0;
604 t5 = ipp2 * t0;
605 for (integer j = 1; j < ipph; j++)
606 {
607 t1 += t2;
608 t3 += t2;
609 t4 += t0;
610 t5 -= t0;
611 for (integer i = 2; i < ido; i += 2)
612 {
613 t6 = idp2 + t1 - i;
614 t7 = i + t3;
615 t8 = i + t4;
616 t9 = i + t5;
617 for (integer k = 0; k < l1; k++)
618 {
619 cc [t7 - 1] = ch [t8 - 1] + ch [t9 - 1];
620 cc [t6 - 1] = ch [t8 - 1] - ch [t9 - 1];
621 cc [t7] = ch [t8] + ch [t9];
622 cc [t6] = ch [t9] - ch [t8];
623 t6 += t10;
624 t7 += t10;
625 t8 += ido;
626 t9 += ido;
627 }
628 }
629 }
630 }
631
drftf1(integer n,FFT_DATA_TYPE * c,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa,integer * ifac)632 static void drftf1 (integer n, FFT_DATA_TYPE * c, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa, integer *ifac)
633 {
634 const integer nf = ifac [1];
635 integer na = 1;
636 integer l2 = n;
637 integer iw = n;
638
639 for (integer k1 = 0; k1 < nf; k1++)
640 {
641 const integer kh = nf - k1;
642 const integer ip = ifac [kh + 1];
643 const integer l1 = l2 / ip;
644 const integer ido = n / l2;
645 const integer idl1 = ido * l1;
646 integer ix2, ix3;
647 iw -= (ip - 1) * ido;
648 na = 1 - na;
649
650 if (ip != 4)
651 goto L102;
652
653 ix2 = iw + ido;
654 ix3 = ix2 + ido;
655 if (na != 0)
656 dradf4 (ido, l1, ch, c, wa + iw - 1, wa + ix2 - 1, wa + ix3 - 1);
657 else
658 dradf4 (ido, l1, c, ch, wa + iw - 1, wa + ix2 - 1, wa + ix3 - 1);
659 goto L110;
660
661 L102:
662 if (ip != 2)
663 goto L104;
664 if (na != 0)
665 goto L103;
666
667 dradf2 (ido, l1, c, ch, wa + iw - 1);
668 goto L110;
669
670 L103:
671 dradf2 (ido, l1, ch, c, wa + iw - 1);
672 goto L110;
673
674 L104:
675 if (ido == 1)
676 na = 1 - na;
677 if (na != 0)
678 goto L109;
679
680 dradfg (ido, ip, l1, idl1, c, c, c, ch, ch, wa + iw - 1);
681 na = 1;
682 goto L110;
683
684 L109:
685 dradfg (ido, ip, l1, idl1, ch, ch, ch, c, c, wa + iw - 1);
686 na = 0;
687
688 L110:
689 l2 = l1;
690 }
691
692 if (na == 1)
693 return;
694
695 for (integer i = 0; i < n; i++)
696 c [i] = ch [i];
697 }
698
dradb2(integer ido,integer l1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa1)699 static void dradb2 (integer ido, integer l1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa1)
700 {
701 const integer t0 = l1 * ido;
702
703 integer t1 = 0;
704 integer t2 = 0;
705 integer t3 = (ido << 1) - 1;
706 for (integer k = 0; k < l1; k++)
707 {
708 ch [t1] = cc [t2] + cc [t3 + t2];
709 ch [t1 + t0] = cc [t2] - cc [t3 + t2];
710 t2 = (t1 += ido) << 1;
711 }
712
713 if (ido < 2)
714 return;
715 if (ido == 2)
716 goto L105;
717
718 t1 = 0;
719 t2 = 0;
720 for (integer k = 0; k < l1; k++)
721 {
722 t3 = t1;
723 integer t4, t5 = (t4 = t2) + (ido << 1);
724 integer t6 = t0 + t1;
725 for (integer i = 2; i < ido; i += 2)
726 {
727 t3 += 2;
728 t4 += 2;
729 t5 -= 2;
730 t6 += 2;
731 ch [t3 - 1] = cc [t4 - 1] + cc [t5 - 1];
732 const double tr2 = cc [t4 - 1] - cc [t5 - 1];
733 ch [t3] = cc [t4] - cc [t5];
734 const double ti2 = cc [t4] + cc [t5];
735 ch [t6 - 1] = wa1 [i - 2] * tr2 - wa1 [i - 1] * ti2;
736 ch [t6] = wa1 [i - 2] * ti2 + wa1 [i - 1] * tr2;
737 }
738 t2 = (t1 += ido) << 1;
739 }
740
741 if (ido % 2 == 1)
742 return;
743
744 L105:
745 t1 = ido - 1;
746 t2 = ido - 1;
747 for (integer k = 0; k < l1; k++)
748 {
749 ch [t1] = cc [t2] + cc [t2];
750 ch [t1 + t0] = -(cc [t2 + 1] + cc [t2 + 1]);
751 t1 += ido;
752 t2 += ido << 1;
753 }
754 }
755
dradb3(integer ido,integer l1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa1,FFT_DATA_TYPE * wa2)756 static void dradb3 (integer ido, integer l1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa1,
757 FFT_DATA_TYPE * wa2)
758 {
759 static constexpr double taur = -.5;
760 static constexpr double taui = .86602540378443864676372317075293618;
761
762 integer t0 = l1 * ido;
763
764 integer t1 = 0;
765 integer t2 = t0 << 1;
766 integer t3 = ido << 1;
767 integer t4 = ido + (ido << 1);
768 integer t5 = 0;
769 for (integer k = 0; k < l1; k++)
770 {
771 const double tr2 = cc [t3 - 1] + cc [t3 - 1];
772 const double cr2 = cc [t5] + (taur * tr2);
773 ch [t1] = cc [t5] + tr2;
774 const double ci3 = taui * (cc [t3] + cc [t3]);
775 ch [t1 + t0] = cr2 - ci3;
776 ch [t1 + t2] = cr2 + ci3;
777 t1 += ido;
778 t3 += t4;
779 t5 += t4;
780 }
781
782 if (ido == 1)
783 return;
784
785 t1 = 0;
786 t3 = ido << 1;
787 for (integer k = 0; k < l1; k++)
788 {
789 integer t7 = t1 + (t1 << 1);
790 integer t6 = (t5 = t7 + t3);
791 integer t8 = t1;
792 integer t9, t10 = (t9 = t1 + t0) + t0;
793
794 for (integer i = 2; i < ido; i += 2)
795 {
796 t5 += 2;
797 t6 -= 2;
798 t7 += 2;
799 t8 += 2;
800 t9 += 2;
801 t10 += 2;
802 const double tr2 = cc [t5 - 1] + cc [t6 - 1];
803 const double cr2 = cc [t7 - 1] + (taur * tr2);
804 ch [t8 - 1] = cc [t7 - 1] + tr2;
805 const double ti2 = cc [t5] - cc [t6];
806 const double ci2 = cc [t7] + (taur * ti2);
807 ch [t8] = cc [t7] + ti2;
808 const double cr3 = taui * (cc [t5 - 1] - cc [t6 - 1]);
809 const double ci3 = taui * (cc [t5] + cc [t6]);
810 const double dr2 = cr2 - ci3;
811 const double dr3 = cr2 + ci3;
812 const double di2 = ci2 + cr3;
813 const double di3 = ci2 - cr3;
814 ch [t9 - 1] = wa1 [i - 2] * dr2 - wa1 [i - 1] * di2;
815 ch [t9] = wa1 [i - 2] * di2 + wa1 [i - 1] * dr2;
816 ch [t10 - 1] = wa2 [i - 2] * dr3 - wa2 [i - 1] * di3;
817 ch [t10] = wa2 [i - 2] * di3 + wa2 [i - 1] * dr3;
818 }
819 t1 += ido;
820 }
821 }
822
dradb4(integer ido,integer l1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa1,FFT_DATA_TYPE * wa2,FFT_DATA_TYPE * wa3)823 static void dradb4 (integer ido, integer l1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa1,
824 FFT_DATA_TYPE * wa2, FFT_DATA_TYPE * wa3)
825 {
826 static constexpr double sqrt2 = 1.4142135623730950488016887242097;
827
828 const integer t0 = l1 * ido;
829
830 integer t1 = 0;
831 integer t2 = ido << 2;
832 integer t3 = 0, t4, t5, t8;
833 const integer t6 = ido << 1;
834 for (integer k = 0; k < l1; k++)
835 {
836 t4 = t3 + t6;
837 t5 = t1;
838 const double tr3 = cc [t4 - 1] + cc [t4 - 1];
839 const double tr4 = cc [t4] + cc [t4];
840 const double tr1 = cc [t3] - cc [(t4 += t6) - 1];
841 const double tr2 = cc [t3] + cc [t4 - 1];
842 ch [t5] = tr2 + tr3;
843 ch [t5 += t0] = tr1 - tr4;
844 ch [t5 += t0] = tr2 - tr3;
845 ch [t5 += t0] = tr1 + tr4;
846 t1 += ido;
847 t3 += t2;
848 }
849
850 if (ido < 2)
851 return;
852 if (ido == 2)
853 goto L105;
854
855 t1 = 0;
856 for (integer k = 0; k < l1; k++)
857 {
858 integer t5 = (t4 = (t3 = (t2 = t1 << 2) + t6)) + t6;
859 integer t7 = t1;
860 for (integer i = 2; i < ido; i += 2)
861 {
862 t2 += 2;
863 t3 += 2;
864 t4 -= 2;
865 t5 -= 2;
866 t7 += 2;
867 const double ti1 = cc [t2] + cc [t5];
868 const double ti2 = cc [t2] - cc [t5];
869 const double ti3 = cc [t3] - cc [t4];
870 const double tr4 = cc [t3] + cc [t4];
871 const double tr1 = cc [t2 - 1] - cc [t5 - 1];
872 const double tr2 = cc [t2 - 1] + cc [t5 - 1];
873 const double ti4 = cc [t3 - 1] - cc [t4 - 1];
874 const double tr3 = cc [t3 - 1] + cc [t4 - 1];
875 ch [t7 - 1] = tr2 + tr3;
876 const double cr3 = tr2 - tr3;
877 ch [t7] = ti2 + ti3;
878 const double ci3 = ti2 - ti3;
879 const double cr2 = tr1 - tr4;
880 const double cr4 = tr1 + tr4;
881 const double ci2 = ti1 + ti4;
882 const double ci4 = ti1 - ti4;
883
884 ch [(t8 = t7 + t0) - 1] = wa1 [i - 2] * cr2 - wa1 [i - 1] * ci2;
885 ch [t8] = wa1 [i - 2] * ci2 + wa1 [i - 1] * cr2;
886 ch [(t8 += t0) - 1] = wa2 [i - 2] * cr3 - wa2 [i - 1] * ci3;
887 ch [t8] = wa2 [i - 2] * ci3 + wa2 [i - 1] * cr3;
888 ch [(t8 += t0) - 1] = wa3 [i - 2] * cr4 - wa3 [i - 1] * ci4;
889 ch [t8] = wa3 [i - 2] * ci4 + wa3 [i - 1] * cr4;
890 }
891 t1 += ido;
892 }
893
894 if (ido % 2 == 1)
895 return;
896
897 L105:
898
899 t1 = ido;
900 t2 = ido << 2;
901 t3 = ido - 1;
902 t4 = ido + (ido << 1);
903 for (integer k = 0; k < l1; k++)
904 {
905 t5 = t3;
906 const double ti1 = cc [t1] + cc [t4];
907 const double ti2 = cc [t4] - cc [t1];
908 const double tr1 = cc [t1 - 1] - cc [t4 - 1];
909 const double tr2 = cc [t1 - 1] + cc [t4 - 1];
910 ch [t5] = tr2 + tr2;
911 ch [t5 += t0] = sqrt2 * (tr1 - ti1);
912 ch [t5 += t0] = ti2 + ti2;
913 ch [t5 += t0] = -sqrt2 * (tr1 + ti1);
914
915 t3 += ido;
916 t1 += t2;
917 t4 += t2;
918 }
919 }
920
dradbg(integer ido,integer ip,integer l1,integer idl1,FFT_DATA_TYPE * cc,FFT_DATA_TYPE * c1,FFT_DATA_TYPE * c2,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * ch2,FFT_DATA_TYPE * wa)921 static void dradbg (integer ido, integer ip, integer l1, integer idl1, FFT_DATA_TYPE * cc, FFT_DATA_TYPE * c1,
922 FFT_DATA_TYPE * c2, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * ch2, FFT_DATA_TYPE * wa)
923 {
924 static constexpr double tpi = 6.28318530717958647692528676655900577;
925 integer is, t1, t2, t3, t4, t5, t6, t7, t8, t9, t11, t12;
926
927 const integer t10 = ip * ido;
928 const integer t0 = l1 * ido;
929 const double arg = tpi / (double) ip;
930 const double dcp = cos (arg);
931 const double dsp = sin (arg);
932 const integer nbd = (ido - 1) >> 1;
933 const integer ipp2 = ip;
934 const integer ipph = (ip + 1) >> 1;
935 if (ido < l1)
936 goto L103;
937
938 t1 = 0;
939 t2 = 0;
940 for (integer k = 0; k < l1; k++)
941 {
942 t3 = t1;
943 t4 = t2;
944 for (integer i = 0; i < ido; i++)
945 {
946 ch [t3] = cc [t4];
947 t3++;
948 t4++;
949 }
950 t1 += ido;
951 t2 += t10;
952 }
953 goto L106;
954
955 L103:
956 t1 = 0;
957 for (integer i = 0; i < ido; i++)
958 {
959 t2 = t1;
960 t3 = t1;
961 for (integer k = 0; k < l1; k++)
962 {
963 ch [t2] = cc [t3];
964 t2 += ido;
965 t3 += t10;
966 }
967 t1++;
968 }
969
970 L106:
971 t1 = 0;
972 t2 = ipp2 * t0;
973 t7 = (t5 = ido << 1);
974 for (integer j = 1; j < ipph; j++)
975 {
976 t1 += t0;
977 t2 -= t0;
978 t3 = t1;
979 t4 = t2;
980 t6 = t5;
981 for (integer k = 0; k < l1; k++)
982 {
983 ch [t3] = cc [t6 - 1] + cc [t6 - 1];
984 ch [t4] = cc [t6] + cc [t6];
985 t3 += ido;
986 t4 += ido;
987 t6 += t10;
988 }
989 t5 += t7;
990 }
991
992 if (ido == 1)
993 goto L116;
994 if (nbd < l1)
995 goto L112;
996
997 t1 = 0;
998 t2 = ipp2 * t0;
999 t7 = 0;
1000 for (integer j = 1; j < ipph; j++)
1001 {
1002 t1 += t0;
1003 t2 -= t0;
1004 t3 = t1;
1005 t4 = t2;
1006
1007 t7 += (ido << 1);
1008 t8 = t7;
1009 for (integer k = 0; k < l1; k++)
1010 {
1011 t5 = t3;
1012 t6 = t4;
1013 t9 = t8;
1014 t11 = t8;
1015 for (integer i = 2; i < ido; i += 2)
1016 {
1017 t5 += 2;
1018 t6 += 2;
1019 t9 += 2;
1020 t11 -= 2;
1021 ch [t5 - 1] = cc [t9 - 1] + cc [t11 - 1];
1022 ch [t6 - 1] = cc [t9 - 1] - cc [t11 - 1];
1023 ch [t5] = cc [t9] - cc [t11];
1024 ch [t6] = cc [t9] + cc [t11];
1025 }
1026 t3 += ido;
1027 t4 += ido;
1028 t8 += t10;
1029 }
1030 }
1031 goto L116;
1032
1033 L112:
1034 t1 = 0;
1035 t2 = ipp2 * t0;
1036 t7 = 0;
1037 for (integer j = 1; j < ipph; j++)
1038 {
1039 t1 += t0;
1040 t2 -= t0;
1041 t3 = t1;
1042 t4 = t2;
1043 t7 += (ido << 1);
1044 t8 = t7;
1045 t9 = t7;
1046 for (integer i = 2; i < ido; i += 2)
1047 {
1048 t3 += 2;
1049 t4 += 2;
1050 t8 += 2;
1051 t9 -= 2;
1052 t5 = t3;
1053 t6 = t4;
1054 t11 = t8;
1055 t12 = t9;
1056 for (integer k = 0; k < l1; k++)
1057 {
1058 ch [t5 - 1] = cc [t11 - 1] + cc [t12 - 1];
1059 ch [t6 - 1] = cc [t11 - 1] - cc [t12 - 1];
1060 ch [t5] = cc [t11] - cc [t12];
1061 ch [t6] = cc [t11] + cc [t12];
1062 t5 += ido;
1063 t6 += ido;
1064 t11 += t10;
1065 t12 += t10;
1066 }
1067 }
1068 }
1069
1070 L116:
1071 double ar1 = 1.;
1072 double ai1 = 0.;
1073 t1 = 0;
1074 t9 = (t2 = ipp2 * idl1);
1075 t3 = (ip - 1) * idl1;
1076 for (integer l = 1; l < ipph; l++)
1077 {
1078 t1 += idl1;
1079 t2 -= idl1;
1080
1081 const double ar1h = dcp * ar1 - dsp * ai1;
1082 ai1 = dcp * ai1 + dsp * ar1;
1083 ar1 = ar1h;
1084 t4 = t1;
1085 t5 = t2;
1086 t6 = 0;
1087 t7 = idl1;
1088 t8 = t3;
1089 for (integer ik = 0; ik < idl1; ik++)
1090 {
1091 c2 [t4++] = ch2 [t6++] + ar1 * ch2 [t7++];
1092 c2 [t5++] = ai1 * ch2 [t8++];
1093 }
1094 const double dc2 = ar1;
1095 const double ds2 = ai1;
1096 double ar2 = ar1;
1097 double ai2 = ai1;
1098
1099 t6 = idl1;
1100 t7 = t9 - idl1;
1101 for (integer j = 2; j < ipph; j++)
1102 {
1103 t6 += idl1;
1104 t7 -= idl1;
1105 const double ar2h = dc2 * ar2 - ds2 * ai2;
1106 ai2 = dc2 * ai2 + ds2 * ar2;
1107 ar2 = ar2h;
1108 t4 = t1;
1109 t5 = t2;
1110 t11 = t6;
1111 t12 = t7;
1112 for (integer ik = 0; ik < idl1; ik++)
1113 {
1114 c2 [t4++] += ar2 * ch2 [t11++];
1115 c2 [t5++] += ai2 * ch2 [t12++];
1116 }
1117 }
1118 }
1119
1120 t1 = 0;
1121 for (integer j = 1; j < ipph; j++)
1122 {
1123 t1 += idl1;
1124 t2 = t1;
1125 for (integer ik = 0; ik < idl1; ik++)
1126 ch2 [ik] += ch2 [t2++];
1127 }
1128
1129 t1 = 0;
1130 t2 = ipp2 * t0;
1131 for (integer j = 1; j < ipph; j++)
1132 {
1133 t1 += t0;
1134 t2 -= t0;
1135 t3 = t1;
1136 t4 = t2;
1137 for (integer k = 0; k < l1; k++)
1138 {
1139 ch [t3] = c1 [t3] - c1 [t4];
1140 ch [t4] = c1 [t3] + c1 [t4];
1141 t3 += ido;
1142 t4 += ido;
1143 }
1144 }
1145
1146 if (ido == 1)
1147 goto L132;
1148 if (nbd < l1)
1149 goto L128;
1150
1151 t1 = 0;
1152 t2 = ipp2 * t0;
1153 for (integer j = 1; j < ipph; j++)
1154 {
1155 t1 += t0;
1156 t2 -= t0;
1157 t3 = t1;
1158 t4 = t2;
1159 for (integer k = 0; k < l1; k++)
1160 {
1161 t5 = t3;
1162 t6 = t4;
1163 for (integer i = 2; i < ido; i += 2)
1164 {
1165 t5 += 2;
1166 t6 += 2;
1167 ch [t5 - 1] = c1 [t5 - 1] - c1 [t6];
1168 ch [t6 - 1] = c1 [t5 - 1] + c1 [t6];
1169 ch [t5] = c1 [t5] + c1 [t6 - 1];
1170 ch [t6] = c1 [t5] - c1 [t6 - 1];
1171 }
1172 t3 += ido;
1173 t4 += ido;
1174 }
1175 }
1176 goto L132;
1177
1178 L128:
1179 t1 = 0;
1180 t2 = ipp2 * t0;
1181 for (integer j = 1; j < ipph; j++)
1182 {
1183 t1 += t0;
1184 t2 -= t0;
1185 t3 = t1;
1186 t4 = t2;
1187 for (integer i = 2; i < ido; i += 2)
1188 {
1189 t3 += 2;
1190 t4 += 2;
1191 t5 = t3;
1192 t6 = t4;
1193 for (integer k = 0; k < l1; k++)
1194 {
1195 ch [t5 - 1] = c1 [t5 - 1] - c1 [t6];
1196 ch [t6 - 1] = c1 [t5 - 1] + c1 [t6];
1197 ch [t5] = c1 [t5] + c1 [t6 - 1];
1198 ch [t6] = c1 [t5] - c1 [t6 - 1];
1199 t5 += ido;
1200 t6 += ido;
1201 }
1202 }
1203 }
1204
1205 L132:
1206 if (ido == 1)
1207 return;
1208
1209 for (integer ik = 0; ik < idl1; ik++)
1210 c2 [ik] = ch2 [ik];
1211
1212 t1 = 0;
1213 for (integer j = 1; j < ip; j++)
1214 {
1215 t2 = (t1 += t0);
1216 for (integer k = 0; k < l1; k++)
1217 {
1218 c1 [t2] = ch [t2];
1219 t2 += ido;
1220 }
1221 }
1222
1223 if (nbd > l1)
1224 goto L139;
1225
1226 is = -ido - 1;
1227 t1 = 0;
1228 for (integer j = 1; j < ip; j++)
1229 {
1230 is += ido;
1231 t1 += t0;
1232 integer idij = is;
1233 t2 = t1;
1234 for (integer i = 2; i < ido; i += 2)
1235 {
1236 t2 += 2;
1237 idij += 2;
1238 t3 = t2;
1239 for (integer k = 0; k < l1; k++)
1240 {
1241 c1 [t3 - 1] = wa [idij - 1] * ch [t3 - 1] - wa [idij] * ch [t3];
1242 c1 [t3] = wa [idij - 1] * ch [t3] + wa [idij] * ch [t3 - 1];
1243 t3 += ido;
1244 }
1245 }
1246 }
1247 return;
1248
1249 L139:
1250 is = -ido - 1;
1251 t1 = 0;
1252 for (integer j = 1; j < ip; j++)
1253 {
1254 is += ido;
1255 t1 += t0;
1256 t2 = t1;
1257 for (integer k = 0; k < l1; k++)
1258 {
1259 integer idij = is;
1260 t3 = t2;
1261 for (integer i = 2; i < ido; i += 2)
1262 {
1263 idij += 2;
1264 t3 += 2;
1265 c1 [t3 - 1] = wa [idij - 1] * ch [t3 - 1] - wa [idij] * ch [t3];
1266 c1 [t3] = wa [idij - 1] * ch [t3] + wa [idij] * ch [t3 - 1];
1267 }
1268 t2 += ido;
1269 }
1270 }
1271 }
1272
drftb1(integer n,FFT_DATA_TYPE * c,FFT_DATA_TYPE * ch,FFT_DATA_TYPE * wa,integer * ifac)1273 static void drftb1 (integer n, FFT_DATA_TYPE * c, FFT_DATA_TYPE * ch, FFT_DATA_TYPE * wa, integer *ifac)
1274 {
1275 const integer nf = ifac [1];
1276 integer na = 0;
1277 integer l1 = 1;
1278 integer iw = 1;
1279
1280 for (integer k1 = 0; k1 < nf; k1++)
1281 {
1282 const integer ip = ifac [k1 + 2];
1283 const integer l2 = ip * l1;
1284 const integer ido = n / l2;
1285 const integer idl1 = ido * l1;
1286 integer ix2, ix3;
1287 if (ip != 4)
1288 goto L103;
1289 ix2 = iw + ido;
1290 ix3 = ix2 + ido;
1291
1292 if (na != 0)
1293 dradb4 (ido, l1, ch, c, wa + iw - 1, wa + ix2 - 1, wa + ix3 - 1);
1294 else
1295 dradb4 (ido, l1, c, ch, wa + iw - 1, wa + ix2 - 1, wa + ix3 - 1);
1296 na = 1 - na;
1297 goto L115;
1298
1299 L103:
1300 if (ip != 2)
1301 goto L106;
1302
1303 if (na != 0)
1304 dradb2 (ido, l1, ch, c, wa + iw - 1);
1305 else
1306 dradb2 (ido, l1, c, ch, wa + iw - 1);
1307 na = 1 - na;
1308 goto L115;
1309
1310 L106:
1311 if (ip != 3)
1312 goto L109;
1313
1314 ix2 = iw + ido;
1315 if (na != 0)
1316 dradb3 (ido, l1, ch, c, wa + iw - 1, wa + ix2 - 1);
1317 else
1318 dradb3 (ido, l1, c, ch, wa + iw - 1, wa + ix2 - 1);
1319 na = 1 - na;
1320 goto L115;
1321
1322 L109:
1323 /* The radix five case can be translated later..... */
1324 /* if(ip!=5)goto L112;
1325
1326 ix2=iw+ido; ix3=ix2+ido; ix4=ix3+ido; if(na!=0)
1327 dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); else
1328 dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); na=1-na;
1329 goto L115;
1330
1331 L112: */
1332 if (na != 0)
1333 dradbg (ido, ip, l1, idl1, ch, ch, ch, c, c, wa + iw - 1);
1334 else
1335 dradbg (ido, ip, l1, idl1, c, c, c, ch, ch, wa + iw - 1);
1336 if (ido == 1)
1337 na = 1 - na;
1338
1339 L115:
1340 l1 = l2;
1341 iw += (ip - 1) * ido;
1342 }
1343
1344 if (na == 0)
1345 return;
1346
1347 for (integer i = 0; i < n; i++)
1348 c [i] = ch [i];
1349 }
1350
1351 /* End of file NUMfft_core.h */
1352