r8brain-free-src
High-quality pro audio sample rate converter library
Loading...
Searching...
No Matches
fft4g.h
Go to the documentation of this file.
1//$ nobt
2//$ nocpp
3
17#ifndef R8B_FFT4G_INCLUDED
18#define R8B_FFT4G_INCLUDED
19
20/*
21Fast Fourier/Cosine/Sine Transform
22 dimension :one
23 data length :power of 2
24 decimation :frequency
25 radix :4, 2
26 data :inplace
27 table :use
28functions
29 cdft: Complex Discrete Fourier Transform
30 rdft: Real Discrete Fourier Transform
31 ddct: Discrete Cosine Transform
32 ddst: Discrete Sine Transform
33 dfct: Cosine Transform of RDFT (Real Symmetric DFT)
34 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
35function prototypes
36 void cdft(int, int, FPType *, int *, FPType *);
37 void rdft(int, int, FPType *, int *, FPType *);
38 void ddct(int, int, FPType *, int *, FPType *);
39 void ddst(int, int, FPType *, int *, FPType *);
40 void dfct(int, FPType *, FPType *, int *, FPType *);
41 void dfst(int, FPType *, FPType *, int *, FPType *);
42
43
44-------- Complex DFT (Discrete Fourier Transform) --------
45 [definition]
46 <case1>
47 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
48 <case2>
49 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
50 (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
51 [usage]
52 <case1>
53 ip[0] = 0; // first time only
54 cdft(2*n, 1, a, ip, w);
55 <case2>
56 ip[0] = 0; // first time only
57 cdft(2*n, -1, a, ip, w);
58 [parameters]
59 2*n :data length (int)
60 n >= 1, n = power of 2
61 a[0...2*n-1] :input/output data (FPType *)
62 input data
63 a[2*j] = Re(x[j]),
64 a[2*j+1] = Im(x[j]), 0<=j<n
65 output data
66 a[2*k] = Re(X[k]),
67 a[2*k+1] = Im(X[k]), 0<=k<n
68 ip[0...*] :work area for bit reversal (int *)
69 length of ip >= 2+sqrt(n)
70 strictly,
71 length of ip >=
72 2+(1<<(int)(log(n+0.5)/log(2))/2).
73 ip[0],ip[1] are pointers of the cos/sin table.
74 w[0...n/2-1] :cos/sin table (FPType *)
75 w[],ip[] are initialized if ip[0] == 0.
76 [remark]
77 Inverse of
78 cdft(2*n, -1, a, ip, w);
79 is
80 cdft(2*n, 1, a, ip, w);
81 for (j = 0; j <= 2 * n - 1; j++) {
82 a[j] *= 1.0 / n;
83 }
84 .
85
86
87-------- Real DFT / Inverse of Real DFT --------
88 [definition]
89 <case1> RDFT
90 R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
91 I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
92 <case2> IRDFT (excluding scale)
93 a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
94 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
95 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
96 [usage]
97 <case1>
98 ip[0] = 0; // first time only
99 rdft(n, 1, a, ip, w);
100 <case2>
101 ip[0] = 0; // first time only
102 rdft(n, -1, a, ip, w);
103 [parameters]
104 n :data length (int)
105 n >= 2, n = power of 2
106 a[0...n-1] :input/output data (FPType *)
107 <case1>
108 output data
109 a[2*k] = R[k], 0<=k<n/2
110 a[2*k+1] = I[k], 0<k<n/2
111 a[1] = R[n/2]
112 <case2>
113 input data
114 a[2*j] = R[j], 0<=j<n/2
115 a[2*j+1] = I[j], 0<j<n/2
116 a[1] = R[n/2]
117 ip[0...*] :work area for bit reversal (int *)
118 length of ip >= 2+sqrt(n/2)
119 strictly,
120 length of ip >=
121 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
122 ip[0],ip[1] are pointers of the cos/sin table.
123 w[0...n/2-1] :cos/sin table (FPType *)
124 w[],ip[] are initialized if ip[0] == 0.
125 [remark]
126 Inverse of
127 rdft(n, 1, a, ip, w);
128 is
129 rdft(n, -1, a, ip, w);
130 for (j = 0; j <= n - 1; j++) {
131 a[j] *= 2.0 / n;
132 }
133 .
134
135
136-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
137 [definition]
138 <case1> IDCT (excluding scale)
139 C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
140 <case2> DCT
141 C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
142 [usage]
143 <case1>
144 ip[0] = 0; // first time only
145 ddct(n, 1, a, ip, w);
146 <case2>
147 ip[0] = 0; // first time only
148 ddct(n, -1, a, ip, w);
149 [parameters]
150 n :data length (int)
151 n >= 2, n = power of 2
152 a[0...n-1] :input/output data (FPType *)
153 output data
154 a[k] = C[k], 0<=k<n
155 ip[0...*] :work area for bit reversal (int *)
156 length of ip >= 2+sqrt(n/2)
157 strictly,
158 length of ip >=
159 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
160 ip[0],ip[1] are pointers of the cos/sin table.
161 w[0...n*5/4-1] :cos/sin table (FPType *)
162 w[],ip[] are initialized if ip[0] == 0.
163 [remark]
164 Inverse of
165 ddct(n, -1, a, ip, w);
166 is
167 a[0] *= 0.5;
168 ddct(n, 1, a, ip, w);
169 for (j = 0; j <= n - 1; j++) {
170 a[j] *= 2.0 / n;
171 }
172 .
173
174
175-------- DST (Discrete Sine Transform) / Inverse of DST --------
176 [definition]
177 <case1> IDST (excluding scale)
178 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
179 <case2> DST
180 S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
181 [usage]
182 <case1>
183 ip[0] = 0; // first time only
184 ddst(n, 1, a, ip, w);
185 <case2>
186 ip[0] = 0; // first time only
187 ddst(n, -1, a, ip, w);
188 [parameters]
189 n :data length (int)
190 n >= 2, n = power of 2
191 a[0...n-1] :input/output data (FPType *)
192 <case1>
193 input data
194 a[j] = A[j], 0<j<n
195 a[0] = A[n]
196 output data
197 a[k] = S[k], 0<=k<n
198 <case2>
199 output data
200 a[k] = S[k], 0<k<n
201 a[0] = S[n]
202 ip[0...*] :work area for bit reversal (int *)
203 length of ip >= 2+sqrt(n/2)
204 strictly,
205 length of ip >=
206 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
207 ip[0],ip[1] are pointers of the cos/sin table.
208 w[0...n*5/4-1] :cos/sin table (FPType *)
209 w[],ip[] are initialized if ip[0] == 0.
210 [remark]
211 Inverse of
212 ddst(n, -1, a, ip, w);
213 is
214 a[0] *= 0.5;
215 ddst(n, 1, a, ip, w);
216 for (j = 0; j <= n - 1; j++) {
217 a[j] *= 2.0 / n;
218 }
219 .
220
221
222-------- Cosine Transform of RDFT (Real Symmetric DFT) --------
223 [definition]
224 C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
225 [usage]
226 ip[0] = 0; // first time only
227 dfct(n, a, t, ip, w);
228 [parameters]
229 n :data length - 1 (int)
230 n >= 2, n = power of 2
231 a[0...n] :input/output data (FPType *)
232 output data
233 a[k] = C[k], 0<=k<=n
234 t[0...n/2] :work area (FPType *)
235 ip[0...*] :work area for bit reversal (int *)
236 length of ip >= 2+sqrt(n/4)
237 strictly,
238 length of ip >=
239 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
240 ip[0],ip[1] are pointers of the cos/sin table.
241 w[0...n*5/8-1] :cos/sin table (FPType *)
242 w[],ip[] are initialized if ip[0] == 0.
243 [remark]
244 Inverse of
245 a[0] *= 0.5;
246 a[n] *= 0.5;
247 dfct(n, a, t, ip, w);
248 is
249 a[0] *= 0.5;
250 a[n] *= 0.5;
251 dfct(n, a, t, ip, w);
252 for (j = 0; j <= n; j++) {
253 a[j] *= 2.0 / n;
254 }
255 .
256
257
258-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
259 [definition]
260 S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
261 [usage]
262 ip[0] = 0; // first time only
263 dfst(n, a, t, ip, w);
264 [parameters]
265 n :data length + 1 (int)
266 n >= 2, n = power of 2
267 a[0...n-1] :input/output data (FPType *)
268 output data
269 a[k] = S[k], 0<k<n
270 (a[0] is used for work area)
271 t[0...n/2-1] :work area (FPType *)
272 ip[0...*] :work area for bit reversal (int *)
273 length of ip >= 2+sqrt(n/4)
274 strictly,
275 length of ip >=
276 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
277 ip[0],ip[1] are pointers of the cos/sin table.
278 w[0...n*5/8-1] :cos/sin table (FPType *)
279 w[],ip[] are initialized if ip[0] == 0.
280 [remark]
281 Inverse of
282 dfst(n, a, t, ip, w);
283 is
284 dfst(n, a, t, ip, w);
285 for (j = 1; j <= n - 1; j++) {
286 a[j] *= 2.0 / n;
287 }
288 .
289
290
291Appendix :
292 The cos/sin table is recalculated when the larger table required.
293 w[] and ip[] are compatible with all routines.
294*/
295
296namespace r8b {
297
305class ooura_fft
306{
307public:
308typedef double FPType;
309
310static void cdft(int n, int isgn, FPType *a, int *ip, FPType *w)
311{
312 if (n > (ip[0] << 2)) {
313 makewt(n >> 2, ip, w);
314 }
315 if (n > 4) {
316 if (isgn >= 0) {
317 bitrv2(n, ip + 2, a);
318 cftfsub(n, a, w);
319 } else {
320 bitrv2conj(n, ip + 2, a);
321 cftbsub(n, a, w);
322 }
323 } else if (n == 4) {
324 cftfsub(n, a, w);
325 }
326}
327
328static void rdft(int n, int isgn, FPType *a, int *ip, FPType *w)
329{
330 int nw, nc;
331 double xi;
332
333 nw = ip[0];
334 if (n > (nw << 2)) {
335 nw = n >> 2;
336 makewt(nw, ip, w);
337 }
338 nc = ip[1];
339 if (n > (nc << 2)) {
340 nc = n >> 2;
341 makect(nc, ip, w + nw);
342 }
343 if (isgn >= 0) {
344 if (n > 4) {
345 bitrv2(n, ip + 2, a);
346 cftfsub(n, a, w);
347 rftfsub(n, a, nc, w + nw);
348 } else if (n == 4) {
349 cftfsub(n, a, w);
350 }
351 xi = a[0] - a[1];
352 a[0] += a[1];
353 a[1] = xi;
354 } else {
355 a[1] = 0.5 * (a[0] - a[1]);
356 a[0] -= a[1];
357 if (n > 4) {
358 rftbsub(n, a, nc, w + nw);
359 bitrv2(n, ip + 2, a);
360 cftbsub(n, a, w);
361 } else if (n == 4) {
362 cftfsub(n, a, w);
363 }
364 }
365}
366
367static void ddct(int n, int isgn, FPType *a, int *ip, FPType *w)
368{
369 int j, nw, nc;
370 double xr;
371
372 nw = ip[0];
373 if (n > (nw << 2)) {
374 nw = n >> 2;
375 makewt(nw, ip, w);
376 }
377 nc = ip[1];
378 if (n > nc) {
379 nc = n;
380 makect(nc, ip, w + nw);
381 }
382 if (isgn < 0) {
383 xr = a[n - 1];
384 for (j = n - 2; j >= 2; j -= 2) {
385 a[j + 1] = a[j] - a[j - 1];
386 a[j] += a[j - 1];
387 }
388 a[1] = a[0] - xr;
389 a[0] += xr;
390 if (n > 4) {
391 rftbsub(n, a, nc, w + nw);
392 bitrv2(n, ip + 2, a);
393 cftbsub(n, a, w);
394 } else if (n == 4) {
395 cftfsub(n, a, w);
396 }
397 }
398 dctsub(n, a, nc, w + nw);
399 if (isgn >= 0) {
400 if (n > 4) {
401 bitrv2(n, ip + 2, a);
402 cftfsub(n, a, w);
403 rftfsub(n, a, nc, w + nw);
404 } else if (n == 4) {
405 cftfsub(n, a, w);
406 }
407 xr = a[0] - a[1];
408 a[0] += a[1];
409 for (j = 2; j < n; j += 2) {
410 a[j - 1] = a[j] - a[j + 1];
411 a[j] += a[j + 1];
412 }
413 a[n - 1] = xr;
414 }
415}
416
417static void ddst(int n, int isgn, FPType *a, int *ip, FPType *w)
418{
419 int j, nw, nc;
420 double xr;
421
422 nw = ip[0];
423 if (n > (nw << 2)) {
424 nw = n >> 2;
425 makewt(nw, ip, w);
426 }
427 nc = ip[1];
428 if (n > nc) {
429 nc = n;
430 makect(nc, ip, w + nw);
431 }
432 if (isgn < 0) {
433 xr = a[n - 1];
434 for (j = n - 2; j >= 2; j -= 2) {
435 a[j + 1] = -a[j] - a[j - 1];
436 a[j] -= a[j - 1];
437 }
438 a[1] = a[0] + xr;
439 a[0] -= xr;
440 if (n > 4) {
441 rftbsub(n, a, nc, w + nw);
442 bitrv2(n, ip + 2, a);
443 cftbsub(n, a, w);
444 } else if (n == 4) {
445 cftfsub(n, a, w);
446 }
447 }
448 dstsub(n, a, nc, w + nw);
449 if (isgn >= 0) {
450 if (n > 4) {
451 bitrv2(n, ip + 2, a);
452 cftfsub(n, a, w);
453 rftfsub(n, a, nc, w + nw);
454 } else if (n == 4) {
455 cftfsub(n, a, w);
456 }
457 xr = a[0] - a[1];
458 a[0] += a[1];
459 for (j = 2; j < n; j += 2) {
460 a[j - 1] = -a[j] - a[j + 1];
461 a[j] -= a[j + 1];
462 }
463 a[n - 1] = -xr;
464 }
465}
466
467static void dfct(int n, FPType *a, FPType *t, int *ip, FPType *w)
468{
469 int j, k, l, m, mh, nw, nc;
470 double xr, xi, yr, yi;
471
472 nw = ip[0];
473 if (n > (nw << 3)) {
474 nw = n >> 3;
475 makewt(nw, ip, w);
476 }
477 nc = ip[1];
478 if (n > (nc << 1)) {
479 nc = n >> 1;
480 makect(nc, ip, w + nw);
481 }
482 m = n >> 1;
483 yi = a[m];
484 xi = a[0] + a[n];
485 a[0] -= a[n];
486 t[0] = xi - yi;
487 t[m] = xi + yi;
488 if (n > 2) {
489 mh = m >> 1;
490 for (j = 1; j < mh; j++) {
491 k = m - j;
492 xr = a[j] - a[n - j];
493 xi = a[j] + a[n - j];
494 yr = a[k] - a[n - k];
495 yi = a[k] + a[n - k];
496 a[j] = xr;
497 a[k] = yr;
498 t[j] = xi - yi;
499 t[k] = xi + yi;
500 }
501 t[mh] = a[mh] + a[n - mh];
502 a[mh] -= a[n - mh];
503 dctsub(m, a, nc, w + nw);
504 if (m > 4) {
505 bitrv2(m, ip + 2, a);
506 cftfsub(m, a, w);
507 rftfsub(m, a, nc, w + nw);
508 } else if (m == 4) {
509 cftfsub(m, a, w);
510 }
511 a[n - 1] = a[0] - a[1];
512 a[1] = a[0] + a[1];
513 for (j = m - 2; j >= 2; j -= 2) {
514 a[2 * j + 1] = a[j] + a[j + 1];
515 a[2 * j - 1] = a[j] - a[j + 1];
516 }
517 l = 2;
518 m = mh;
519 while (m >= 2) {
520 dctsub(m, t, nc, w + nw);
521 if (m > 4) {
522 bitrv2(m, ip + 2, t);
523 cftfsub(m, t, w);
524 rftfsub(m, t, nc, w + nw);
525 } else if (m == 4) {
526 cftfsub(m, t, w);
527 }
528 a[n - l] = t[0] - t[1];
529 a[l] = t[0] + t[1];
530 k = 0;
531 for (j = 2; j < m; j += 2) {
532 k += l << 2;
533 a[k - l] = t[j] - t[j + 1];
534 a[k + l] = t[j] + t[j + 1];
535 }
536 l <<= 1;
537 mh = m >> 1;
538 for (j = 0; j < mh; j++) {
539 k = m - j;
540 t[j] = t[m + k] - t[m + j];
541 t[k] = t[m + k] + t[m + j];
542 }
543 t[mh] = t[m + mh];
544 m = mh;
545 }
546 a[l] = t[0];
547 a[n] = t[2] - t[1];
548 a[0] = t[2] + t[1];
549 } else {
550 a[1] = a[0];
551 a[2] = t[0];
552 a[0] = t[1];
553 }
554}
555
556static void dfst(int n, FPType *a, FPType *t, int *ip, FPType *w)
557{
558 int j, k, l, m, mh, nw, nc;
559 double xr, xi, yr, yi;
560
561 nw = ip[0];
562 if (n > (nw << 3)) {
563 nw = n >> 3;
564 makewt(nw, ip, w);
565 }
566 nc = ip[1];
567 if (n > (nc << 1)) {
568 nc = n >> 1;
569 makect(nc, ip, w + nw);
570 }
571 if (n > 2) {
572 m = n >> 1;
573 mh = m >> 1;
574 for (j = 1; j < mh; j++) {
575 k = m - j;
576 xr = a[j] + a[n - j];
577 xi = a[j] - a[n - j];
578 yr = a[k] + a[n - k];
579 yi = a[k] - a[n - k];
580 a[j] = xr;
581 a[k] = yr;
582 t[j] = xi + yi;
583 t[k] = xi - yi;
584 }
585 t[0] = a[mh] - a[n - mh];
586 a[mh] += a[n - mh];
587 a[0] = a[m];
588 dstsub(m, a, nc, w + nw);
589 if (m > 4) {
590 bitrv2(m, ip + 2, a);
591 cftfsub(m, a, w);
592 rftfsub(m, a, nc, w + nw);
593 } else if (m == 4) {
594 cftfsub(m, a, w);
595 }
596 a[n - 1] = a[1] - a[0];
597 a[1] = a[0] + a[1];
598 for (j = m - 2; j >= 2; j -= 2) {
599 a[2 * j + 1] = a[j] - a[j + 1];
600 a[2 * j - 1] = -a[j] - a[j + 1];
601 }
602 l = 2;
603 m = mh;
604 while (m >= 2) {
605 dstsub(m, t, nc, w + nw);
606 if (m > 4) {
607 bitrv2(m, ip + 2, t);
608 cftfsub(m, t, w);
609 rftfsub(m, t, nc, w + nw);
610 } else if (m == 4) {
611 cftfsub(m, t, w);
612 }
613 a[n - l] = t[1] - t[0];
614 a[l] = t[0] + t[1];
615 k = 0;
616 for (j = 2; j < m; j += 2) {
617 k += l << 2;
618 a[k - l] = -t[j] - t[j + 1];
619 a[k + l] = t[j] - t[j + 1];
620 }
621 l <<= 1;
622 mh = m >> 1;
623 for (j = 1; j < mh; j++) {
624 k = m - j;
625 t[j] = t[m + k] + t[m + j];
626 t[k] = t[m + k] - t[m + j];
627 }
628 t[0] = t[m + mh];
629 m = mh;
630 }
631 a[l] = t[0];
632 }
633 a[0] = 0;
634}
635
636/* -------- initializing routines -------- */
637
638static void makewt(int nw, int *ip, FPType *w)
639{
640 int j, nwh;
641 double delta, x, y;
642
643 ip[0] = nw;
644 ip[1] = 1;
645 if (nw > 2) {
646 nwh = nw >> 1;
647 delta = atan(1.0) / nwh;
648 w[0] = 1;
649 w[1] = 0;
650 w[nwh] = cos(delta * nwh);
651 w[nwh + 1] = w[nwh];
652 if (nwh > 2) {
653 for (j = 2; j < nwh; j += 2) {
654 x = cos(delta * j);
655 y = sin(delta * j);
656 w[j] = x;
657 w[j + 1] = y;
658 w[nw - j] = y;
659 w[nw - j + 1] = x;
660 }
661 bitrv2(nw, ip + 2, w);
662 }
663 }
664}
665
666static void makect(int nc, int *ip, FPType *c)
667{
668 int j, nch;
669 double delta;
670
671 ip[1] = nc;
672 if (nc > 1) {
673 nch = nc >> 1;
674 delta = atan(1.0) / nch;
675 c[0] = cos(delta * nch);
676 c[nch] = 0.5 * c[0];
677 for (j = 1; j < nch; j++) {
678 c[j] = 0.5 * cos(delta * j);
679 c[nc - j] = 0.5 * sin(delta * j);
680 }
681 }
682}
683
684/* -------- child routines -------- */
685
686static void bitrv2(int n, int *ip, FPType *a)
687{
688 int j, j1, k, k1, l, m, m2;
689 double xr, xi, yr, yi;
690
691 ip[0] = 0;
692 l = n;
693 m = 1;
694 while ((m << 3) < l) {
695 l >>= 1;
696 for (j = 0; j < m; j++) {
697 ip[m + j] = ip[j] + l;
698 }
699 m <<= 1;
700 }
701 m2 = 2 * m;
702 if ((m << 3) == l) {
703 for (k = 0; k < m; k++) {
704 for (j = 0; j < k; j++) {
705 j1 = 2 * j + ip[k];
706 k1 = 2 * k + ip[j];
707 xr = a[j1];
708 xi = a[j1 + 1];
709 yr = a[k1];
710 yi = a[k1 + 1];
711 a[j1] = yr;
712 a[j1 + 1] = yi;
713 a[k1] = xr;
714 a[k1 + 1] = xi;
715 j1 += m2;
716 k1 += 2 * m2;
717 xr = a[j1];
718 xi = a[j1 + 1];
719 yr = a[k1];
720 yi = a[k1 + 1];
721 a[j1] = yr;
722 a[j1 + 1] = yi;
723 a[k1] = xr;
724 a[k1 + 1] = xi;
725 j1 += m2;
726 k1 -= m2;
727 xr = a[j1];
728 xi = a[j1 + 1];
729 yr = a[k1];
730 yi = a[k1 + 1];
731 a[j1] = yr;
732 a[j1 + 1] = yi;
733 a[k1] = xr;
734 a[k1 + 1] = xi;
735 j1 += m2;
736 k1 += 2 * m2;
737 xr = a[j1];
738 xi = a[j1 + 1];
739 yr = a[k1];
740 yi = a[k1 + 1];
741 a[j1] = yr;
742 a[j1 + 1] = yi;
743 a[k1] = xr;
744 a[k1 + 1] = xi;
745 }
746 j1 = 2 * k + m2 + ip[k];
747 k1 = j1 + m2;
748 xr = a[j1];
749 xi = a[j1 + 1];
750 yr = a[k1];
751 yi = a[k1 + 1];
752 a[j1] = yr;
753 a[j1 + 1] = yi;
754 a[k1] = xr;
755 a[k1 + 1] = xi;
756 }
757 } else {
758 for (k = 1; k < m; k++) {
759 for (j = 0; j < k; j++) {
760 j1 = 2 * j + ip[k];
761 k1 = 2 * k + ip[j];
762 xr = a[j1];
763 xi = a[j1 + 1];
764 yr = a[k1];
765 yi = a[k1 + 1];
766 a[j1] = yr;
767 a[j1 + 1] = yi;
768 a[k1] = xr;
769 a[k1 + 1] = xi;
770 j1 += m2;
771 k1 += m2;
772 xr = a[j1];
773 xi = a[j1 + 1];
774 yr = a[k1];
775 yi = a[k1 + 1];
776 a[j1] = yr;
777 a[j1 + 1] = yi;
778 a[k1] = xr;
779 a[k1 + 1] = xi;
780 }
781 }
782 }
783}
784
785static void bitrv2conj(int n, int *ip, FPType *a)
786{
787 int j, j1, k, k1, l, m, m2;
788 double xr, xi, yr, yi;
789
790 ip[0] = 0;
791 l = n;
792 m = 1;
793 while ((m << 3) < l) {
794 l >>= 1;
795 for (j = 0; j < m; j++) {
796 ip[m + j] = ip[j] + l;
797 }
798 m <<= 1;
799 }
800 m2 = 2 * m;
801 if ((m << 3) == l) {
802 for (k = 0; k < m; k++) {
803 for (j = 0; j < k; j++) {
804 j1 = 2 * j + ip[k];
805 k1 = 2 * k + ip[j];
806 xr = a[j1];
807 xi = -a[j1 + 1];
808 yr = a[k1];
809 yi = -a[k1 + 1];
810 a[j1] = yr;
811 a[j1 + 1] = yi;
812 a[k1] = xr;
813 a[k1 + 1] = xi;
814 j1 += m2;
815 k1 += 2 * m2;
816 xr = a[j1];
817 xi = -a[j1 + 1];
818 yr = a[k1];
819 yi = -a[k1 + 1];
820 a[j1] = yr;
821 a[j1 + 1] = yi;
822 a[k1] = xr;
823 a[k1 + 1] = xi;
824 j1 += m2;
825 k1 -= m2;
826 xr = a[j1];
827 xi = -a[j1 + 1];
828 yr = a[k1];
829 yi = -a[k1 + 1];
830 a[j1] = yr;
831 a[j1 + 1] = yi;
832 a[k1] = xr;
833 a[k1 + 1] = xi;
834 j1 += m2;
835 k1 += 2 * m2;
836 xr = a[j1];
837 xi = -a[j1 + 1];
838 yr = a[k1];
839 yi = -a[k1 + 1];
840 a[j1] = yr;
841 a[j1 + 1] = yi;
842 a[k1] = xr;
843 a[k1 + 1] = xi;
844 }
845 k1 = 2 * k + ip[k];
846 a[k1 + 1] = -a[k1 + 1];
847 j1 = k1 + m2;
848 k1 = j1 + m2;
849 xr = a[j1];
850 xi = -a[j1 + 1];
851 yr = a[k1];
852 yi = -a[k1 + 1];
853 a[j1] = yr;
854 a[j1 + 1] = yi;
855 a[k1] = xr;
856 a[k1 + 1] = xi;
857 k1 += m2;
858 a[k1 + 1] = -a[k1 + 1];
859 }
860 } else {
861 a[1] = -a[1];
862 a[m2 + 1] = -a[m2 + 1];
863 for (k = 1; k < m; k++) {
864 for (j = 0; j < k; j++) {
865 j1 = 2 * j + ip[k];
866 k1 = 2 * k + ip[j];
867 xr = a[j1];
868 xi = -a[j1 + 1];
869 yr = a[k1];
870 yi = -a[k1 + 1];
871 a[j1] = yr;
872 a[j1 + 1] = yi;
873 a[k1] = xr;
874 a[k1 + 1] = xi;
875 j1 += m2;
876 k1 += m2;
877 xr = a[j1];
878 xi = -a[j1 + 1];
879 yr = a[k1];
880 yi = -a[k1 + 1];
881 a[j1] = yr;
882 a[j1 + 1] = yi;
883 a[k1] = xr;
884 a[k1 + 1] = xi;
885 }
886 k1 = 2 * k + ip[k];
887 a[k1 + 1] = -a[k1 + 1];
888 a[k1 + m2 + 1] = -a[k1 + m2 + 1];
889 }
890 }
891}
892
893static void cftfsub(int n, FPType *a, const FPType *w)
894{
895 int j, j1, j2, j3, l;
896 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
897
898 l = 2;
899 if (n > 8) {
900 cft1st(n, a, w);
901 l = 8;
902 while ((l << 2) < n) {
903 cftmdl(n, l, a, w);
904 l <<= 2;
905 }
906 }
907 if ((l << 2) == n) {
908 for (j = 0; j < l; j += 2) {
909 j1 = j + l;
910 j2 = j1 + l;
911 j3 = j2 + l;
912 x0r = a[j] + a[j1];
913 x0i = a[j + 1] + a[j1 + 1];
914 x1r = a[j] - a[j1];
915 x1i = a[j + 1] - a[j1 + 1];
916 x2r = a[j2] + a[j3];
917 x2i = a[j2 + 1] + a[j3 + 1];
918 x3r = a[j2] - a[j3];
919 x3i = a[j2 + 1] - a[j3 + 1];
920 a[j] = x0r + x2r;
921 a[j + 1] = x0i + x2i;
922 a[j2] = x0r - x2r;
923 a[j2 + 1] = x0i - x2i;
924 a[j1] = x1r - x3i;
925 a[j1 + 1] = x1i + x3r;
926 a[j3] = x1r + x3i;
927 a[j3 + 1] = x1i - x3r;
928 }
929 } else {
930 for (j = 0; j < l; j += 2) {
931 j1 = j + l;
932 x0r = a[j] - a[j1];
933 x0i = a[j + 1] - a[j1 + 1];
934 a[j] += a[j1];
935 a[j + 1] += a[j1 + 1];
936 a[j1] = x0r;
937 a[j1 + 1] = x0i;
938 }
939 }
940}
941
942static void cftbsub(int n, FPType *a, const FPType *w)
943{
944 int j, j1, j2, j3, l;
945 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
946
947 l = 2;
948 if (n > 8) {
949 cft1st(n, a, w);
950 l = 8;
951 while ((l << 2) < n) {
952 cftmdl(n, l, a, w);
953 l <<= 2;
954 }
955 }
956 if ((l << 2) == n) {
957 for (j = 0; j < l; j += 2) {
958 j1 = j + l;
959 j2 = j1 + l;
960 j3 = j2 + l;
961 x0r = a[j] + a[j1];
962 x0i = -a[j + 1] - a[j1 + 1];
963 x1r = a[j] - a[j1];
964 x1i = -a[j + 1] + a[j1 + 1];
965 x2r = a[j2] + a[j3];
966 x2i = a[j2 + 1] + a[j3 + 1];
967 x3r = a[j2] - a[j3];
968 x3i = a[j2 + 1] - a[j3 + 1];
969 a[j] = x0r + x2r;
970 a[j + 1] = x0i - x2i;
971 a[j2] = x0r - x2r;
972 a[j2 + 1] = x0i + x2i;
973 a[j1] = x1r - x3i;
974 a[j1 + 1] = x1i - x3r;
975 a[j3] = x1r + x3i;
976 a[j3 + 1] = x1i + x3r;
977 }
978 } else {
979 for (j = 0; j < l; j += 2) {
980 j1 = j + l;
981 x0r = a[j] - a[j1];
982 x0i = -a[j + 1] + a[j1 + 1];
983 a[j] += a[j1];
984 a[j + 1] = -a[j + 1] - a[j1 + 1];
985 a[j1] = x0r;
986 a[j1 + 1] = x0i;
987 }
988 }
989}
990
991static void cft1st(int n, FPType *a, const FPType *w)
992{
993 int j, k1, k2;
994 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
995 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
996
997 x0r = a[0] + a[2];
998 x0i = a[1] + a[3];
999 x1r = a[0] - a[2];
1000 x1i = a[1] - a[3];
1001 x2r = a[4] + a[6];
1002 x2i = a[5] + a[7];
1003 x3r = a[4] - a[6];
1004 x3i = a[5] - a[7];
1005 a[0] = x0r + x2r;
1006 a[1] = x0i + x2i;
1007 a[4] = x0r - x2r;
1008 a[5] = x0i - x2i;
1009 a[2] = x1r - x3i;
1010 a[3] = x1i + x3r;
1011 a[6] = x1r + x3i;
1012 a[7] = x1i - x3r;
1013 wk1r = w[2];
1014 x0r = a[8] + a[10];
1015 x0i = a[9] + a[11];
1016 x1r = a[8] - a[10];
1017 x1i = a[9] - a[11];
1018 x2r = a[12] + a[14];
1019 x2i = a[13] + a[15];
1020 x3r = a[12] - a[14];
1021 x3i = a[13] - a[15];
1022 a[8] = x0r + x2r;
1023 a[9] = x0i + x2i;
1024 a[12] = x2i - x0i;
1025 a[13] = x0r - x2r;
1026 x0r = x1r - x3i;
1027 x0i = x1i + x3r;
1028 a[10] = wk1r * (x0r - x0i);
1029 a[11] = wk1r * (x0r + x0i);
1030 x0r = x3i + x1r;
1031 x0i = x3r - x1i;
1032 a[14] = wk1r * (x0i - x0r);
1033 a[15] = wk1r * (x0i + x0r);
1034 k1 = 0;
1035 for (j = 16; j < n; j += 16) {
1036 k1 += 2;
1037 k2 = 2 * k1;
1038 wk2r = w[k1];
1039 wk2i = w[k1 + 1];
1040 wk1r = w[k2];
1041 wk1i = w[k2 + 1];
1042 wk3r = wk1r - 2 * wk2i * wk1i;
1043 wk3i = 2 * wk2i * wk1r - wk1i;
1044 x0r = a[j] + a[j + 2];
1045 x0i = a[j + 1] + a[j + 3];
1046 x1r = a[j] - a[j + 2];
1047 x1i = a[j + 1] - a[j + 3];
1048 x2r = a[j + 4] + a[j + 6];
1049 x2i = a[j + 5] + a[j + 7];
1050 x3r = a[j + 4] - a[j + 6];
1051 x3i = a[j + 5] - a[j + 7];
1052 a[j] = x0r + x2r;
1053 a[j + 1] = x0i + x2i;
1054 x0r -= x2r;
1055 x0i -= x2i;
1056 a[j + 4] = wk2r * x0r - wk2i * x0i;
1057 a[j + 5] = wk2r * x0i + wk2i * x0r;
1058 x0r = x1r - x3i;
1059 x0i = x1i + x3r;
1060 a[j + 2] = wk1r * x0r - wk1i * x0i;
1061 a[j + 3] = wk1r * x0i + wk1i * x0r;
1062 x0r = x1r + x3i;
1063 x0i = x1i - x3r;
1064 a[j + 6] = wk3r * x0r - wk3i * x0i;
1065 a[j + 7] = wk3r * x0i + wk3i * x0r;
1066 wk1r = w[k2 + 2];
1067 wk1i = w[k2 + 3];
1068 wk3r = wk1r - 2 * wk2r * wk1i;
1069 wk3i = 2 * wk2r * wk1r - wk1i;
1070 x0r = a[j + 8] + a[j + 10];
1071 x0i = a[j + 9] + a[j + 11];
1072 x1r = a[j + 8] - a[j + 10];
1073 x1i = a[j + 9] - a[j + 11];
1074 x2r = a[j + 12] + a[j + 14];
1075 x2i = a[j + 13] + a[j + 15];
1076 x3r = a[j + 12] - a[j + 14];
1077 x3i = a[j + 13] - a[j + 15];
1078 a[j + 8] = x0r + x2r;
1079 a[j + 9] = x0i + x2i;
1080 x0r -= x2r;
1081 x0i -= x2i;
1082 a[j + 12] = -wk2i * x0r - wk2r * x0i;
1083 a[j + 13] = -wk2i * x0i + wk2r * x0r;
1084 x0r = x1r - x3i;
1085 x0i = x1i + x3r;
1086 a[j + 10] = wk1r * x0r - wk1i * x0i;
1087 a[j + 11] = wk1r * x0i + wk1i * x0r;
1088 x0r = x1r + x3i;
1089 x0i = x1i - x3r;
1090 a[j + 14] = wk3r * x0r - wk3i * x0i;
1091 a[j + 15] = wk3r * x0i + wk3i * x0r;
1092 }
1093}
1094
1095static void cftmdl(int n, int l, FPType *a, const FPType *w)
1096{
1097 int j, j1, j2, j3, k, k1, k2, m, m2;
1098 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
1099 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1100
1101 m = l << 2;
1102 for (j = 0; j < l; j += 2) {
1103 j1 = j + l;
1104 j2 = j1 + l;
1105 j3 = j2 + l;
1106 x0r = a[j] + a[j1];
1107 x0i = a[j + 1] + a[j1 + 1];
1108 x1r = a[j] - a[j1];
1109 x1i = a[j + 1] - a[j1 + 1];
1110 x2r = a[j2] + a[j3];
1111 x2i = a[j2 + 1] + a[j3 + 1];
1112 x3r = a[j2] - a[j3];
1113 x3i = a[j2 + 1] - a[j3 + 1];
1114 a[j] = x0r + x2r;
1115 a[j + 1] = x0i + x2i;
1116 a[j2] = x0r - x2r;
1117 a[j2 + 1] = x0i - x2i;
1118 a[j1] = x1r - x3i;
1119 a[j1 + 1] = x1i + x3r;
1120 a[j3] = x1r + x3i;
1121 a[j3 + 1] = x1i - x3r;
1122 }
1123 wk1r = w[2];
1124 for (j = m; j < l + m; j += 2) {
1125 j1 = j + l;
1126 j2 = j1 + l;
1127 j3 = j2 + l;
1128 x0r = a[j] + a[j1];
1129 x0i = a[j + 1] + a[j1 + 1];
1130 x1r = a[j] - a[j1];
1131 x1i = a[j + 1] - a[j1 + 1];
1132 x2r = a[j2] + a[j3];
1133 x2i = a[j2 + 1] + a[j3 + 1];
1134 x3r = a[j2] - a[j3];
1135 x3i = a[j2 + 1] - a[j3 + 1];
1136 a[j] = x0r + x2r;
1137 a[j + 1] = x0i + x2i;
1138 a[j2] = x2i - x0i;
1139 a[j2 + 1] = x0r - x2r;
1140 x0r = x1r - x3i;
1141 x0i = x1i + x3r;
1142 a[j1] = wk1r * (x0r - x0i);
1143 a[j1 + 1] = wk1r * (x0r + x0i);
1144 x0r = x3i + x1r;
1145 x0i = x3r - x1i;
1146 a[j3] = wk1r * (x0i - x0r);
1147 a[j3 + 1] = wk1r * (x0i + x0r);
1148 }
1149 k1 = 0;
1150 m2 = 2 * m;
1151 for (k = m2; k < n; k += m2) {
1152 k1 += 2;
1153 k2 = 2 * k1;
1154 wk2r = w[k1];
1155 wk2i = w[k1 + 1];
1156 wk1r = w[k2];
1157 wk1i = w[k2 + 1];
1158 wk3r = wk1r - 2 * wk2i * wk1i;
1159 wk3i = 2 * wk2i * wk1r - wk1i;
1160 for (j = k; j < l + k; j += 2) {
1161 j1 = j + l;
1162 j2 = j1 + l;
1163 j3 = j2 + l;
1164 x0r = a[j] + a[j1];
1165 x0i = a[j + 1] + a[j1 + 1];
1166 x1r = a[j] - a[j1];
1167 x1i = a[j + 1] - a[j1 + 1];
1168 x2r = a[j2] + a[j3];
1169 x2i = a[j2 + 1] + a[j3 + 1];
1170 x3r = a[j2] - a[j3];
1171 x3i = a[j2 + 1] - a[j3 + 1];
1172 a[j] = x0r + x2r;
1173 a[j + 1] = x0i + x2i;
1174 x0r -= x2r;
1175 x0i -= x2i;
1176 a[j2] = wk2r * x0r - wk2i * x0i;
1177 a[j2 + 1] = wk2r * x0i + wk2i * x0r;
1178 x0r = x1r - x3i;
1179 x0i = x1i + x3r;
1180 a[j1] = wk1r * x0r - wk1i * x0i;
1181 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1182 x0r = x1r + x3i;
1183 x0i = x1i - x3r;
1184 a[j3] = wk3r * x0r - wk3i * x0i;
1185 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1186 }
1187 wk1r = w[k2 + 2];
1188 wk1i = w[k2 + 3];
1189 wk3r = wk1r - 2 * wk2r * wk1i;
1190 wk3i = 2 * wk2r * wk1r - wk1i;
1191 for (j = k + m; j < l + (k + m); j += 2) {
1192 j1 = j + l;
1193 j2 = j1 + l;
1194 j3 = j2 + l;
1195 x0r = a[j] + a[j1];
1196 x0i = a[j + 1] + a[j1 + 1];
1197 x1r = a[j] - a[j1];
1198 x1i = a[j + 1] - a[j1 + 1];
1199 x2r = a[j2] + a[j3];
1200 x2i = a[j2 + 1] + a[j3 + 1];
1201 x3r = a[j2] - a[j3];
1202 x3i = a[j2 + 1] - a[j3 + 1];
1203 a[j] = x0r + x2r;
1204 a[j + 1] = x0i + x2i;
1205 x0r -= x2r;
1206 x0i -= x2i;
1207 a[j2] = -wk2i * x0r - wk2r * x0i;
1208 a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
1209 x0r = x1r - x3i;
1210 x0i = x1i + x3r;
1211 a[j1] = wk1r * x0r - wk1i * x0i;
1212 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
1213 x0r = x1r + x3i;
1214 x0i = x1i - x3r;
1215 a[j3] = wk3r * x0r - wk3i * x0i;
1216 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
1217 }
1218 }
1219}
1220
1221static void rftfsub(int n, FPType *a, int nc, const FPType *c)
1222{
1223 int j, k, kk, ks, m;
1224 double wkr, wki, xr, xi, yr, yi;
1225
1226 m = n >> 1;
1227 ks = 2 * nc / m;
1228 kk = 0;
1229 for (j = 2; j < m; j += 2) {
1230 k = n - j;
1231 kk += ks;
1232 wkr = 0.5 - c[nc - kk];
1233 wki = c[kk];
1234 xr = a[j] - a[k];
1235 xi = a[j + 1] + a[k + 1];
1236 yr = wkr * xr - wki * xi;
1237 yi = wkr * xi + wki * xr;
1238 a[j] -= yr;
1239 a[j + 1] -= yi;
1240 a[k] += yr;
1241 a[k + 1] -= yi;
1242 }
1243}
1244
1245static void rftbsub(int n, FPType *a, int nc, const FPType *c)
1246{
1247 int j, k, kk, ks, m;
1248 double wkr, wki, xr, xi, yr, yi;
1249
1250 a[1] = -a[1];
1251 m = n >> 1;
1252 ks = 2 * nc / m;
1253 kk = 0;
1254 for (j = 2; j < m; j += 2) {
1255 k = n - j;
1256 kk += ks;
1257 wkr = 0.5 - c[nc - kk];
1258 wki = c[kk];
1259 xr = a[j] - a[k];
1260 xi = a[j + 1] + a[k + 1];
1261 yr = wkr * xr + wki * xi;
1262 yi = wkr * xi - wki * xr;
1263 a[j] -= yr;
1264 a[j + 1] = yi - a[j + 1];
1265 a[k] += yr;
1266 a[k + 1] = yi - a[k + 1];
1267 }
1268 a[m + 1] = -a[m + 1];
1269}
1270
1271static void dctsub(int n, FPType *a, int nc, const FPType *c)
1272{
1273 int j, k, kk, ks, m;
1274 double wkr, wki, xr;
1275
1276 m = n >> 1;
1277 ks = nc / n;
1278 kk = 0;
1279 for (j = 1; j < m; j++) {
1280 k = n - j;
1281 kk += ks;
1282 wkr = c[kk] - c[nc - kk];
1283 wki = c[kk] + c[nc - kk];
1284 xr = wki * a[j] - wkr * a[k];
1285 a[j] = wkr * a[j] + wki * a[k];
1286 a[k] = xr;
1287 }
1288 a[m] *= c[0];
1289}
1290
1291static void dstsub(int n, FPType *a, int nc, const FPType *c)
1292{
1293 int j, k, kk, ks, m;
1294 double wkr, wki, xr;
1295
1296 m = n >> 1;
1297 ks = nc / n;
1298 kk = 0;
1299 for (j = 1; j < m; j++) {
1300 k = n - j;
1301 kk += ks;
1302 wkr = c[kk] - c[nc - kk];
1303 wki = c[kk] + c[nc - kk];
1304 xr = wki * a[k] - wkr * a[j];
1305 a[k] = wkr * a[k] + wki * a[j];
1306 a[j] = xr;
1307 }
1308 a[m] *= c[0];
1309}
1310};
1311
1312} // namespace r8b
1313
1314#endif // R8B_FFT4G_INCLUDED
r8b
The "r8brain-free-src" library namespace.
Definition: CDSPBlockConvolver.h:21
r8b::ooura_fft
A wrapper class around Takuya OOURA's FFT functions.
Definition: fft4g.h:306
Generated by doxygen 1.9.5