Mercurial > hg > js-dsp-test

--- /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
--- /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
--- /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)
+  -
+  - 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.
+  -->
+<!DOCTYPE html>
+<html>
+<head>
+<meta charset="UTF-8">
+<title>FFT and convolution test (JavaScript)</title>
+</head>
+<body>
+<script src="fft.js" type="application/javascript"></script>
+<pre>
+<script>
+/* Main and test functions */
+
+function main() {
+    // Test power-of-2 size FFTs
+    for (var i = 0; i <= 12; i++)
+        testFft(1 << i);
+
+    // Test small size FFTs
+    for (var i = 0; i < 30; i++)
+        testFft(i);
+
+    // Test diverse size FFTs
+    var prev = 0;
+    for (var i = 0; i <= 100; i++) {
+        var n = Math.round(Math.pow(1500, i / 100.0));
+        if (n > prev) {
+            testFft(n);
+            prev = n;
+        }
+    }
+
+    // Test power-of-2 size convolutions
+    for (var i = 0; i <= 12; i++)
+        testConvolution(1 << i);
+
+    // Test diverse size convolutions
+    prev = 0;
+    for (var i = 0; i <= 100; i++) {
+        var n = Math.round(Math.pow(1500, i / 100.0));
+        if (n > prev) {
+            testConvolution(n);
+            prev = n;
+        }
+    }
+
+    document.write("\nMax log err = " + maxLogError.toFixed(1));
+    document.write("\nTest " + (maxLogError < -10 ? "passed" : "failed"));
+}
+
+
+function testFft(size) {
+    var inputreal = randomReals(size);
+    var inputimag = randomReals(size);
+
+    var refoutreal = new Array(size);
+    var refoutimag = new Array(size);
+    naiveDft(inputreal, inputimag, refoutreal, refoutimag, false);
+
+    var actualoutreal = inputreal.slice(0);
+    var actualoutimag = inputimag.slice(0);
+    transform(actualoutreal, actualoutimag);
+
+    document.write("fftsize=" + size + "  logerr=" + log10RmsErr(refoutreal, refoutimag, actualoutreal, actualoutimag).toFixed(1) + "\n");
+}
+
+
+function testConvolution(size) {
+    var input0real = randomReals(size);
+    var input0imag = randomReals(size);
+
+    var input1real = randomReals(size);
+    var input1imag = randomReals(size);
+
+    var refoutreal = new Array(size);
+    var refoutimag = new Array(size);
+    naiveConvolve(input0real, input0imag, input1real, input1imag, refoutreal, refoutimag);
+
+    var actualoutreal = new Array(size);
+    var actualoutimag = new Array(size);
+    convolveComplex(input0real, input0imag, input1real, input1imag, actualoutreal, actualoutimag);
+
+    document.write("convsize=" + size + "  logerr=" + log10RmsErr(refoutreal, refoutimag, actualoutreal, actualoutimag).toFixed(1) + "\n");
+}
+
+
+/* Naive reference computation functions */
+
+function naiveDft(inreal, inimag, outreal, outimag, inverse) {
+    if (inreal.length != inimag.length || inreal.length != outreal.length || outreal.length != outimag.length)
+        throw "Mismatched lengths";
+
+    var n = inreal.length;
+    var coef = (inverse ? 2 : -2) * Math.PI;
+    for (var k = 0; k < n; k++) {  // For each output element
+        var sumreal = 0;
+        var sumimag = 0;
+        for (var t = 0; t < n; t++) {  // For each input element
+            var angle = coef * (t * k % n) / n;  // This is more accurate than t * k
+            sumreal += inreal[t]*Math.cos(angle) - inimag[t]*Math.sin(angle);
+            sumimag += inreal[t]*Math.sin(angle) + inimag[t]*Math.cos(angle);
+        }
+        outreal[k] = sumreal;
+        outimag[k] = sumimag;
+    }
+}
+
+
+function naiveConvolve(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;
+    for (var i = 0; i < n; i++) {
+        var sumreal = 0;
+        var sumimag = 0;
+        for (var j = 0; j < n; j++) {
+            var k = (i - j + n) % n;
+            sumreal += xreal[k] * yreal[j] - ximag[k] * yimag[j];
+            sumimag += xreal[k] * yimag[j] + ximag[k] * yreal[j];
+        }
+        outreal[i] = sumreal;
+        outimag[i] = sumimag;
+    }
+}
+
+
+/* Utility functions */
+
+var maxLogError = Number.NEGATIVE_INFINITY;
+
+function log10RmsErr(xreal, ximag, yreal, yimag) {
+    if (xreal.length != ximag.length || xreal.length != yreal.length || yreal.length != yimag.length)
+        throw "Mismatched lengths";
+
+    var err = 0;
+    for (var i = 0; i < xreal.length; i++)
+        err += (xreal[i] - yreal[i]) * (xreal[i] - yreal[i]) + (ximag[i] - yimag[i]) * (ximag[i] - yimag[i]);
+    err = Math.sqrt(err / Math.max(xreal.length, 1));  // Now this is a root mean square (RMS) error
+    err = err > 0 ? Math.log(err) / Math.log(10) : -99;
+    maxLogError = Math.max(err, maxLogError);
+    return err;
+}
+
+
+function randomReals(size) {
+    var result = new Array(size);
+    for (var i = 0; i < result.length; i++)
+        result[i] = Math.random() * 2 - 1;
+    return result;
+}
+
+
+main();
+</script>
+</pre>
+</body>
+</html>
--- /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;
+    }
+}