# HG changeset patch
# User Chris Cannam
# Date 1443711058 -3600
# Node ID d7c216b6a84ffa4c703f157754faba236ff68b5b
Pull in some FFT implementations for test
diff -r 000000000000 -r d7c216b6a84f .hgsub
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgsub Thu Oct 01 15:50:58 2015 +0100
@@ -0,0 +1,4 @@
+fft-test/fft.js = [git]https://github.com/JensNockert/fft.js
+fft-test/timbre.js = [git]https://github.com/mohayonao/timbre.js
+fft-test/jsfft = [git]https://github.com/dntj/jsfft
+fft-test/dsp.js = [git]https://github.com/corbanbrook/dsp.js
diff -r 000000000000 -r d7c216b6a84f .hgsubstate
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgsubstate Thu Oct 01 15:50:58 2015 +0100
@@ -0,0 +1,4 @@
+a7b2e97b1385a43083e50ed6dc81d697f0e57e28 fft-test/dsp.js
+dd96e6bec5464fe6b30c91ddbfdef88ea9d70bd1 fft-test/fft.js
+b82257673c97f5f4bfc12fe445dcf1cf2169b46e fft-test/jsfft
+7d94de567ea1d758599b1e29a11afff304d44c55 fft-test/timbre.js
diff -r 000000000000 -r d7c216b6a84f fft/nayuki/fft-test.html
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/fft/nayuki/fft-test.html Thu Oct 01 15:50:58 2015 +0100
@@ -0,0 +1,180 @@
+
+
+
+
+
+FFT and convolution test (JavaScript)
+
+
+
+
+
+
+
+
diff -r 000000000000 -r d7c216b6a84f fft/nayuki/fft.js
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/fft/nayuki/fft.js Thu Oct 01 15:50:58 2015 +0100
@@ -0,0 +1,215 @@
+/*
+ * Free FFT and convolution (JavaScript)
+ *
+ * Copyright (c) 2014 Project Nayuki
+ * http://www.nayuki.io/page/free-small-fft-in-multiple-languages
+ *
+ * (MIT License)
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of
+ * this software and associated documentation files (the "Software"), to deal in
+ * the Software without restriction, including without limitation the rights to
+ * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+ * the Software, and to permit persons to whom the Software is furnished to do so,
+ * subject to the following conditions:
+ * - The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ * - The Software is provided "as is", without warranty of any kind, express or
+ * implied, including but not limited to the warranties of merchantability,
+ * fitness for a particular purpose and noninfringement. In no event shall the
+ * authors or copyright holders be liable for any claim, damages or other
+ * liability, whether in an action of contract, tort or otherwise, arising from,
+ * out of or in connection with the Software or the use or other dealings in the
+ * Software.
+ */
+
+"use strict";
+
+
+/*
+ * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
+ * The vector can have any length. This is a wrapper function.
+ */
+function transform(real, imag) {
+ if (real.length != imag.length)
+ throw "Mismatched lengths";
+
+ var n = real.length;
+ if (n == 0)
+ return;
+ else if ((n & (n - 1)) == 0) // Is power of 2
+ transformRadix2(real, imag);
+ else // More complicated algorithm for arbitrary sizes
+ transformBluestein(real, imag);
+}
+
+
+/*
+ * Computes the inverse discrete Fourier transform (IDFT) of the given complex vector, storing the result back into the vector.
+ * The vector can have any length. This is a wrapper function. This transform does not perform scaling, so the inverse is not a true inverse.
+ */
+function inverseTransform(real, imag) {
+ transform(imag, real);
+}
+
+
+/*
+ * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
+ * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm.
+ */
+function transformRadix2(real, imag) {
+ // Initialization
+ if (real.length != imag.length)
+ throw "Mismatched lengths";
+ var n = real.length;
+ if (n == 1) // Trivial transform
+ return;
+ var levels = -1;
+ for (var i = 0; i < 32; i++) {
+ if (1 << i == n)
+ levels = i; // Equal to log2(n)
+ }
+ if (levels == -1)
+ throw "Length is not a power of 2";
+ var cosTable = new Array(n / 2);
+ var sinTable = new Array(n / 2);
+ for (var i = 0; i < n / 2; i++) {
+ cosTable[i] = Math.cos(2 * Math.PI * i / n);
+ sinTable[i] = Math.sin(2 * Math.PI * i / n);
+ }
+
+ // Bit-reversed addressing permutation
+ for (var i = 0; i < n; i++) {
+ var j = reverseBits(i, levels);
+ if (j > i) {
+ var temp = real[i];
+ real[i] = real[j];
+ real[j] = temp;
+ temp = imag[i];
+ imag[i] = imag[j];
+ imag[j] = temp;
+ }
+ }
+
+ // Cooley-Tukey decimation-in-time radix-2 FFT
+ for (var size = 2; size <= n; size *= 2) {
+ var halfsize = size / 2;
+ var tablestep = n / size;
+ for (var i = 0; i < n; i += size) {
+ for (var j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
+ var tpre = real[j+halfsize] * cosTable[k] + imag[j+halfsize] * sinTable[k];
+ var tpim = -real[j+halfsize] * sinTable[k] + imag[j+halfsize] * cosTable[k];
+ real[j + halfsize] = real[j] - tpre;
+ imag[j + halfsize] = imag[j] - tpim;
+ real[j] += tpre;
+ imag[j] += tpim;
+ }
+ }
+ }
+
+ // Returns the integer whose value is the reverse of the lowest 'bits' bits of the integer 'x'.
+ function reverseBits(x, bits) {
+ var y = 0;
+ for (var i = 0; i < bits; i++) {
+ y = (y << 1) | (x & 1);
+ x >>>= 1;
+ }
+ return y;
+ }
+}
+
+
+/*
+ * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
+ * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function.
+ * Uses Bluestein's chirp z-transform algorithm.
+ */
+function transformBluestein(real, imag) {
+ // Find a power-of-2 convolution length m such that m >= n * 2 + 1
+ if (real.length != imag.length)
+ throw "Mismatched lengths";
+ var n = real.length;
+ var m = 1;
+ while (m < n * 2 + 1)
+ m *= 2;
+
+ // Trignometric tables
+ var cosTable = new Array(n);
+ var sinTable = new Array(n);
+ for (var i = 0; i < n; i++) {
+ var j = i * i % (n * 2); // This is more accurate than j = i * i
+ cosTable[i] = Math.cos(Math.PI * j / n);
+ sinTable[i] = Math.sin(Math.PI * j / n);
+ }
+
+ // Temporary vectors and preprocessing
+ var areal = new Array(m);
+ var aimag = new Array(m);
+ for (var i = 0; i < n; i++) {
+ areal[i] = real[i] * cosTable[i] + imag[i] * sinTable[i];
+ aimag[i] = -real[i] * sinTable[i] + imag[i] * cosTable[i];
+ }
+ for (var i = n; i < m; i++)
+ areal[i] = aimag[i] = 0;
+ var breal = new Array(m);
+ var bimag = new Array(m);
+ breal[0] = cosTable[0];
+ bimag[0] = sinTable[0];
+ for (var i = 1; i < n; i++) {
+ breal[i] = breal[m - i] = cosTable[i];
+ bimag[i] = bimag[m - i] = sinTable[i];
+ }
+ for (var i = n; i <= m - n; i++)
+ breal[i] = bimag[i] = 0;
+
+ // Convolution
+ var creal = new Array(m);
+ var cimag = new Array(m);
+ convolveComplex(areal, aimag, breal, bimag, creal, cimag);
+
+ // Postprocessing
+ for (var i = 0; i < n; i++) {
+ real[i] = creal[i] * cosTable[i] + cimag[i] * sinTable[i];
+ imag[i] = -creal[i] * sinTable[i] + cimag[i] * cosTable[i];
+ }
+}
+
+
+/*
+ * Computes the circular convolution of the given real vectors. Each vector's length must be the same.
+ */
+function convolveReal(x, y, out) {
+ if (x.length != y.length || x.length != out.length)
+ throw "Mismatched lengths";
+ var zeros = new Array(x.length);
+ for (var i = 0; i < zeros.length; i++)
+ zeros[i] = 0;
+ convolveComplex(x, zeros, y, zeros.slice(0), out, zeros.slice(0));
+}
+
+
+/*
+ * Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
+ */
+function convolveComplex(xreal, ximag, yreal, yimag, outreal, outimag) {
+ if (xreal.length != ximag.length || xreal.length != yreal.length || yreal.length != yimag.length || xreal.length != outreal.length || outreal.length != outimag.length)
+ throw "Mismatched lengths";
+
+ var n = xreal.length;
+ xreal = xreal.slice(0);
+ ximag = ximag.slice(0);
+ yreal = yreal.slice(0);
+ yimag = yimag.slice(0);
+
+ transform(xreal, ximag);
+ transform(yreal, yimag);
+ for (var i = 0; i < n; i++) {
+ var temp = xreal[i] * yreal[i] - ximag[i] * yimag[i];
+ ximag[i] = ximag[i] * yreal[i] + xreal[i] * yimag[i];
+ xreal[i] = temp;
+ }
+ inverseTransform(xreal, ximag);
+ for (var i = 0; i < n; i++) { // Scaling (because this FFT implementation omits it)
+ outreal[i] = xreal[i] / n;
+ outimag[i] = ximag[i] / n;
+ }
+}