39 "category",
"Analysis",
40 "description",
"FFT_ooura",
43 static FactoryRegistrator<ProcessingFactory, FFT_ooura>
reg =
metadata;
46 bool FFT_ooura::ConcreteConfigure(
const ProcessingConfig& c) {
51 AddConfigErrorMessage(
"Configure failed: FFT Ooura algorithm only works for input buffers that are a power of two!");
55 if (
mSize != oldSize) {
65 void FFT_ooura::ReleaseMemory() {
66 if (ip) {
delete[] ip; ip = 0; }
67 if (w) {
delete[] w; w = 0; }
70 void FFT_ooura::SetupMemory() {
71 int ipSize = (
int)(2+(1<<(
int)(log(
mSize/2+0.5)/log(2.0))/2));
73 for (
int i=0; i<ipSize; i++) ip[i] = 0;
77 for (
int i=0; i<wSize; i++) w[i] = 0;
110 "FFT_ooura: Do(): Not in execution mode");
112 "FFT_ooura: Do(): Not a power of two");
114 out.SetSpectralRange(in.GetSampleRate()/2);
118 inbuffer = in.GetBuffer().GetPtr();
122 for (
int i=0; i<
mSize; i++)
128 inbuffer = in.GetBuffer().GetPtr();
132 for (
int i=0; i<
mSize; i++)
140 inbuffer = in.GetBuffer().GetPtr();
144 for (
int i=0; i<
mSize; i++)
151 CLAM_ASSERT(
false,
"FFT_ooura: Do(): Inconsistent state");
158 void FFT_ooura::ToComplex(
Spectrum &out) {
162 outbuffer[0].SetImag(0);
164 outbuffer[
mSize/2].SetImag(0);
166 for (
int i=1; i<
mSize/2; i++) {
174 void FFT_ooura::ToOther(Spectrum &out)
176 if(out.HasComplexArray()) {
221 a[1] = 0.5 * (a[0] - a[1]);
250 for (j = 2; j < nwh; j += 2) {
275 for (j = 1; j < nch; j++) {
277 c[nc - j] = 0.5 *
CLAM_sin(delta * j);
284 int j, j1, k, k1, l, m, m2;
285 TData xr, xi, yr, yi;
290 while ((m << 3) < l) {
292 for (j = 0; j < m; j++) {
293 ip[m + j] = ip[j] + l;
299 for (k = 0; k < m; k++) {
300 for (j = 0; j < k; j++) {
342 j1 = 2 * k + m2 + ip[k];
354 for (k = 1; k < m; k++) {
355 for (j = 0; j < k; j++) {
383 int j, j1, j2, j3, l;
384 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
390 while ((l << 2) < n) {
396 for (j = 0; j < l; j += 2) {
401 x0i = a[j + 1] + a[j1 + 1];
403 x1i = a[j + 1] - a[j1 + 1];
405 x2i = a[j2 + 1] + a[j3 + 1];
407 x3i = a[j2 + 1] - a[j3 + 1];
409 a[j + 1] = x0i + x2i;
411 a[j2 + 1] = x0i - x2i;
413 a[j1 + 1] = x1i + x3r;
415 a[j3 + 1] = x1i - x3r;
418 for (j = 0; j < l; j += 2) {
421 x0i = a[j + 1] - a[j1 + 1];
423 a[j + 1] += a[j1 + 1];
432 int j, j1, j2, j3, l;
433 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
439 while ((l << 2) < n) {
445 for (j = 0; j < l; j += 2) {
450 x0i = -a[j + 1] - a[j1 + 1];
452 x1i = -a[j + 1] + a[j1 + 1];
454 x2i = a[j2 + 1] + a[j3 + 1];
456 x3i = a[j2 + 1] - a[j3 + 1];
458 a[j + 1] = x0i - x2i;
460 a[j2 + 1] = x0i + x2i;
462 a[j1 + 1] = x1i - x3r;
464 a[j3 + 1] = x1i + x3r;
467 for (j = 0; j < l; j += 2) {
470 x0i = -a[j + 1] + a[j1 + 1];
472 a[j + 1] = -a[j + 1] - a[j1 + 1];
482 TData wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
483 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
516 a[10] = wk1r * (x0r - x0i);
517 a[11] = wk1r * (x0r + x0i);
520 a[14] = wk1r * (x0i - x0r);
521 a[15] = wk1r * (x0i + x0r);
523 for (j = 16; j < n; j += 16) {
530 wk3r = wk1r - 2 * wk2i * wk1i;
531 wk3i = 2 * wk2i * wk1r - wk1i;
532 x0r = a[j] + a[j + 2];
533 x0i = a[j + 1] + a[j + 3];
534 x1r = a[j] - a[j + 2];
535 x1i = a[j + 1] - a[j + 3];
536 x2r = a[j + 4] + a[j + 6];
537 x2i = a[j + 5] + a[j + 7];
538 x3r = a[j + 4] - a[j + 6];
539 x3i = a[j + 5] - a[j + 7];
541 a[j + 1] = x0i + x2i;
544 a[j + 4] = wk2r * x0r - wk2i * x0i;
545 a[j + 5] = wk2r * x0i + wk2i * x0r;
548 a[j + 2] = wk1r * x0r - wk1i * x0i;
549 a[j + 3] = wk1r * x0i + wk1i * x0r;
552 a[j + 6] = wk3r * x0r - wk3i * x0i;
553 a[j + 7] = wk3r * x0i + wk3i * x0r;
556 wk3r = wk1r - 2 * wk2r * wk1i;
557 wk3i = 2 * wk2r * wk1r - wk1i;
558 x0r = a[j + 8] + a[j + 10];
559 x0i = a[j + 9] + a[j + 11];
560 x1r = a[j + 8] - a[j + 10];
561 x1i = a[j + 9] - a[j + 11];
562 x2r = a[j + 12] + a[j + 14];
563 x2i = a[j + 13] + a[j + 15];
564 x3r = a[j + 12] - a[j + 14];
565 x3i = a[j + 13] - a[j + 15];
566 a[j + 8] = x0r + x2r;
567 a[j + 9] = x0i + x2i;
570 a[j + 12] = -wk2i * x0r - wk2r * x0i;
571 a[j + 13] = -wk2i * x0i + wk2r * x0r;
574 a[j + 10] = wk1r * x0r - wk1i * x0i;
575 a[j + 11] = wk1r * x0i + wk1i * x0r;
578 a[j + 14] = wk3r * x0r - wk3i * x0i;
579 a[j + 15] = wk3r * x0i + wk3i * x0r;
586 int j, j1, j2, j3, k, k1, k2, m, m2;
587 TData wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
588 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
591 for (j = 0; j < l; j += 2) {
596 x0i = a[j + 1] + a[j1 + 1];
598 x1i = a[j + 1] - a[j1 + 1];
600 x2i = a[j2 + 1] + a[j3 + 1];
602 x3i = a[j2 + 1] - a[j3 + 1];
604 a[j + 1] = x0i + x2i;
606 a[j2 + 1] = x0i - x2i;
608 a[j1 + 1] = x1i + x3r;
610 a[j3 + 1] = x1i - x3r;
613 for (j = m; j < l + m; j += 2) {
618 x0i = a[j + 1] + a[j1 + 1];
620 x1i = a[j + 1] - a[j1 + 1];
622 x2i = a[j2 + 1] + a[j3 + 1];
624 x3i = a[j2 + 1] - a[j3 + 1];
626 a[j + 1] = x0i + x2i;
628 a[j2 + 1] = x0r - x2r;
631 a[j1] = wk1r * (x0r - x0i);
632 a[j1 + 1] = wk1r * (x0r + x0i);
635 a[j3] = wk1r * (x0i - x0r);
636 a[j3 + 1] = wk1r * (x0i + x0r);
640 for (k = m2; k < n; k += m2) {
647 wk3r = wk1r - 2 * wk2i * wk1i;
648 wk3i = 2 * wk2i * wk1r - wk1i;
649 for (j = k; j < l + k; j += 2) {
654 x0i = a[j + 1] + a[j1 + 1];
656 x1i = a[j + 1] - a[j1 + 1];
658 x2i = a[j2 + 1] + a[j3 + 1];
660 x3i = a[j2 + 1] - a[j3 + 1];
662 a[j + 1] = x0i + x2i;
665 a[j2] = wk2r * x0r - wk2i * x0i;
666 a[j2 + 1] = wk2r * x0i + wk2i * x0r;
669 a[j1] = wk1r * x0r - wk1i * x0i;
670 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
673 a[j3] = wk3r * x0r - wk3i * x0i;
674 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
678 wk3r = wk1r - 2 * wk2r * wk1i;
679 wk3i = 2 * wk2r * wk1r - wk1i;
680 for (j = k + m; j < l + (k + m); j += 2) {
685 x0i = a[j + 1] + a[j1 + 1];
687 x1i = a[j + 1] - a[j1 + 1];
689 x2i = a[j2 + 1] + a[j3 + 1];
691 x3i = a[j2 + 1] - a[j3 + 1];
693 a[j + 1] = x0i + x2i;
696 a[j2] = -wk2i * x0r - wk2r * x0i;
697 a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
700 a[j1] = wk1r * x0r - wk1i * x0i;
701 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
704 a[j3] = wk3r * x0r - wk3i * x0i;
705 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
714 TData wkr, wki, xr, xi, yr, yi;
719 for (j = 2; j < m; j += 2) {
722 wkr = 0.5 - c[nc - kk];
725 xi = a[j + 1] + a[k + 1];
726 yr = wkr * xr - wki * xi;
727 yi = wkr * xi + wki * xr;
739 TData wkr, wki, xr, xi, yr, yi;
745 for (j = 2; j < m; j += 2) {
748 wkr = 0.5 - c[nc - kk];
751 xi = a[j + 1] + a[k + 1];
752 yr = wkr * xr + wki * xi;
753 yi = wkr * xi - wki * xr;
755 a[j + 1] = yi - a[j + 1];
757 a[k + 1] = yi - a[k + 1];
759 a[m + 1] = -a[m + 1];