comparison fft.cpp @ 9:a1539d4e3b08

* Tidy code, start doc

+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#ifdef VST
+#include "devuvuzelator-vst.h"
+#else
+#ifdef LADSPA
+#include "devuvuzelator-ladspa.h"
+#else
+#error Need to define either VST or LADSPA
+#endif
+#endif
+
+// FFT implementation by Don Cross, public domain.
+// This version scales the forward transform.
+
+void Devuvuzelator::fft(unsigned int n, bool inverse,
+                        double *ri, double *ii, double *ro, double *io)
+{
+    if (!ri || !ro || !io) return;
+
+    unsigned int bits;
+    unsigned int i, j, k, m;
+    unsigned int blockSize, blockEnd;
+
+    double tr, ti;
+
+    if (n < 2) return;
+    if (n & (n-1)) return;
+
+    double angle = 2.0 * M_PI;
+    if (inverse) angle = -angle;
+
+    for (i = 0; ; ++i) {
+	if (n & (1 << i)) {
+	    bits = i;
+	    break;
+	}
+    }
+
+    static unsigned int tableSize = 0;
+    static int *table = 0;
+
+    if (tableSize != n) {
+
+	delete[] table;
+
+	table = new int[n];
+
+	for (i = 0; i < n; ++i) {
+	
+	    m = i;
+
+	    for (j = k = 0; j < bits; ++j) {
+		k = (k << 1) | (m & 1);
+		m >>= 1;
+	    }
+
+	    table[i] = k;
+	}
+
+	tableSize = n;
+    }
+
+    if (ii) {
+	for (i = 0; i < n; ++i) {
+	    ro[table[i]] = ri[i];
+	    io[table[i]] = ii[i];
+	}
+    } else {
+	for (i = 0; i < n; ++i) {
+	    ro[table[i]] = ri[i];
+	    io[table[i]] = 0.0;
+	}
+    }
+
+    blockEnd = 1;
+
+    for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
+
+	double delta = angle / (double)blockSize;
+	double sm2 = -sin(-2 * delta);
+	double sm1 = -sin(-delta);
+	double cm2 = cos(-2 * delta);
+	double cm1 = cos(-delta);
+	double w = 2 * cm1;
+	double ar[3], ai[3];
+
+	for (i = 0; i < n; i += blockSize) {
+
+	    ar[2] = cm2;
+	    ar[1] = cm1;
+
+	    ai[2] = sm2;
+	    ai[1] = sm1;
+
+	    for (j = i, m = 0; m < blockEnd; j++, m++) {
+
+		ar[0] = w * ar[1] - ar[2];
+		ar[2] = ar[1];
+		ar[1] = ar[0];
+
+		ai[0] = w * ai[1] - ai[2];
+		ai[2] = ai[1];
+		ai[1] = ai[0];
+
+		k = j + blockEnd;
+		tr = ar[0] * ro[k] - ai[0] * io[k];
+		ti = ar[0] * io[k] + ai[0] * ro[k];
+
+		ro[k] = ro[j] - tr;
+		io[k] = io[j] - ti;
+
+		ro[j] += tr;
+		io[j] += ti;
+	    }
+	}
+
+	blockEnd = blockSize;
+    }
+
+    if (!inverse) {
+
+	double denom = (double)n;
+
+	for (i = 0; i < n; i++) {
+	    ro[i] /= denom;
+	    io[i] /= denom;
+	}
+    }
+}
+