# HG changeset patch # User Chris Cannam # Date 1444221998 -3600 # Node ID 66f9fd5ac61156c7d8493e65f3ff35f844b55f0d # Parent e705de983b67d3899a87143068e6265d3be46d0e Bring in some more of the third-party code diff -r e705de983b67 -r 66f9fd5ac611 .hgsub --- a/.hgsub Wed Oct 07 08:57:11 2015 +0100 +++ b/.hgsub Wed Oct 07 13:46:38 2015 +0100 @@ -1,4 +0,0 @@ -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 e705de983b67 -r 66f9fd5ac611 .hgsubstate --- a/.hgsubstate Wed Oct 07 08:57:11 2015 +0100 +++ b/.hgsubstate Wed Oct 07 13:46:38 2015 +0100 @@ -1,4 +0,0 @@ -a7b2e97b1385a43083e50ed6dc81d697f0e57e28 fft-test/dsp.js -dd96e6bec5464fe6b30c91ddbfdef88ea9d70bd1 fft-test/fft.js -b82257673c97f5f4bfc12fe445dcf1cf2169b46e fft-test/jsfft -7d94de567ea1d758599b1e29a11afff304d44c55 fft-test/timbre.js diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/LICENSE --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/LICENSE Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,20 @@ +/* Copyright (c) 2012, Jens Nockert , Jussi Kalliokoski + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/Makefile --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/Makefile Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,18 @@ +TARGETS := lib/complex.js lib/node.erb.js lib/node.js + +all: $(TARGETS) +node: lib/node.js + +lib/complex.js: src/complex.erb.js + erb $^ > $@ + +lib/node.erb.js: src/node.erb.js + erb $^ > $@ + +lib/node.js: lib/node.erb.js + erb $^ > $@ + +clean: + rm -rf $(TARGETS) + +.PHONY: all clean diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/README.md --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/README.md Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,41 @@ +fft.js +================================================================================ + +FFT in JavaScript, it works, I think. + +No promises, but I tested it against Wolfram Alpha once, and it was reasonably accurate. + +There are optimized kernels for prime factors, 2, 3, 4, so if you want high performance, use lengths that are a factor of those. + +Notice that the DFT is not normalized, so `ifft(fft(x)) / n ~= x` + + +Usage +--------------------------------------------------------------------------------- + +```javascript + +/* Create a new FFT object */ + +var fft = new FFT.complex(n, inverse) + +/* Output and input should be float arrays (of the right length), type is either 'complex' (default) or 'real' */ +fft.process(output, outputOffset, outputStride, input, inputOffset, inputStride, type) + +/* Or the simplified interface, which just sets the offsets to 0, and the strides to 1 */ +fft.simple(output, input, type) + +``` + + +Installing via npm +--------------------------------------------------------------------------------- + +You can also install via npm, the name is `fft` in the registy. + + +Credits +--------------------------------------------------------------------------------- + +I was too lazy to calculate the butterflies myself, so they are inspired by [kissfft](http://sourceforge.net/projects/kissfft/), which is a small library for doing discrete fourier transforms. + diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/lib/complex.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/lib/complex.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,292 @@ +/* Copyright (c) 2012, Jens Nockert , Jussi Kalliokoski + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +if (!FFT) { + var FFT = {} +} + +void function (namespace) { + "use strict" + + function butterfly2(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + + for (var i = 0; i < m; i++) { + var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] + var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] + + var t1_r = t[2 * ((0) + (fStride) * (i))], t1_i = t[2 * ((0) + (fStride) * (i)) + 1] + + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i + var r1_r = s0_r - v1_r, r1_i = s0_i - v1_i + + output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i + output[2 * ((outputOffset) + (outputStride) * (i + m))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] = r1_i + } + } + + function butterfly3(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m + var fStride1 = fStride, fStride2 = 2 * fStride + + var e = t[2 * ((0) + (fStride) * (m)) + 1] + + for (var i = 0; i < m; i++) { + var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] + + var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] + var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] + var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1] + var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r + + var i0_r = v1_r + v2_r, i0_i = v1_i + v2_i + + var r0_r = s0_r + i0_r, r0_i = s0_i + i0_i + output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i + + var i1_r = s0_r - i0_r * 0.5 + var i1_i = s0_i - i0_i * 0.5 + + var i2_r = (v1_r - v2_r) * e + var i2_i = (v1_i - v2_i) * e + + var r1_r = i1_r - i2_i + var r1_i = i1_i + i2_r + output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i + + var r2_r = i1_r + i2_i + var r2_i = i1_i - i2_r + output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i + } + } + + function butterfly4(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m, m3 = 3 * m + var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride + + for (var i = 0; i < m; i++) { + var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] + + var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] + var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] + var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1] + var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r + + var s3_r = output[2 * ((outputOffset) + (outputStride) * (i + m3))], s3_i = output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] + var t3_r = t[2 * ((0) + (fStride3) * (i))], t3_i = t[2 * ((0) + (fStride3) * (i)) + 1] + var v3_r = s3_r * t3_r - s3_i * t3_i, v3_i = s3_r * t3_i + s3_i * t3_r + + var i0_r = s0_r + v2_r, i0_i = s0_i + v2_i + var i1_r = s0_r - v2_r, i1_i = s0_i - v2_i + var i2_r = v1_r + v3_r, i2_i = v1_i + v3_i + var i3_r = v1_r - v3_r, i3_i = v1_i - v3_i + + var r0_r = i0_r + i2_r, r0_i = i0_i + i2_i + + if (state.inverse) { + var r1_r = i1_r - i3_i + var r1_i = i1_i + i3_r + } else { + var r1_r = i1_r + i3_i + var r1_i = i1_i - i3_r + } + + var r2_r = i0_r - i2_r, r2_i = i0_i - i2_i + + if (state.inverse) { + var r3_r = i1_r + i3_i + var r3_i = i1_i - i3_r + } else { + var r3_r = i1_r - i3_i + var r3_i = i1_i + i3_r + } + + output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i + output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i + output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i + output[2 * ((outputOffset) + (outputStride) * (i + m3))] = r3_r, output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] = r3_i + } + } + + function butterfly(output, outputOffset, outputStride, fStride, state, m, p) { + var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p) + + for (var u = 0; u < m; u++) { + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + var x0_r = output[2 * ((outputOffset) + (outputStride) * (k))], x0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1] + scratch[2 * (q1)] = x0_r, scratch[2 * (q1) + 1] = x0_i + } + + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + var tOffset = 0 + + var x0_r = scratch[2 * (0)], x0_i = scratch[2 * (0) + 1] + output[2 * ((outputOffset) + (outputStride) * (k))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = x0_i + + for (var q = 1; q < p; q++) { + tOffset = (tOffset + fStride * k) % n + + var s0_r = output[2 * ((outputOffset) + (outputStride) * (k))], s0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1] + + var s1_r = scratch[2 * (q)], s1_i = scratch[2 * (q) + 1] + var t1_r = t[2 * (tOffset)], t1_i = t[2 * (tOffset) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i + output[2 * ((outputOffset) + (outputStride) * (k))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = r0_i + } + } + } + } + + function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { + var p = factors.shift() + var m = factors.shift() + + if (m == 1) { + for (var i = 0; i < p * m; i++) { + var x0_r = f[2 * ((fOffset) + (fStride * inputStride) * (i))], x0_i = f[2 * ((fOffset) + (fStride * inputStride) * (i)) + 1] + output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i + } + } else { + for (var i = 0; i < p; i++) { + work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) + } + } + + switch (p) { + case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break + case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break + case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break + default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break + } + } + + var complex = function (n, inverse) { + if (arguments.length < 2) { + throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'") + } + + var n = ~~n, inverse = !!inverse + + if (n < 1) { + throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") + } + + var state = { + n: n, + inverse: inverse, + + factors: [], + twiddle: new Float64Array(2 * n), + scratch: new Float64Array(2 * n) + } + + var t = state.twiddle, theta = 2 * Math.PI / n + + for (var i = 0; i < n; i++) { + if (inverse) { + var phase = theta * i + } else { + var phase = -theta * i + } + + t[2 * (i)] = Math.cos(phase) + t[2 * (i) + 1] = Math.sin(phase) + } + + var p = 4, v = Math.floor(Math.sqrt(n)) + + while (n > 1) { + while (n % p) { + switch (p) { + case 4: p = 2; break + case 2: p = 3; break + default: p += 2; break + } + + if (p > v) { + p = n + } + } + + n /= p + + state.factors.push(p) + state.factors.push(n) + } + + this.state = state + } + + complex.prototype.simple = function (output, input, t) { + this.process(output, 0, 1, input, 0, 1, t) + } + + complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) { + var outputStride = ~~outputStride, inputStride = ~~inputStride + + var type = t == 'real' ? t : 'complex' + + if (outputStride < 1) { + throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") + } + + if (inputStride < 1) { + throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") + } + + if (type == 'real') { + for (var i = 0; i < this.state.n; i++) { + var x0_r = input[inputOffset + inputStride * i] + var x0_i = 0.0 + + this.state.scratch[2 * (i)] = x0_r, this.state.scratch[2 * (i) + 1] = x0_i + } + + work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state) + } else { + if (input == output) { + work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + + for (var i = 0; i < this.state.n; i++) { + var x0_r = this.state.scratch[2 * (i)], x0_i = this.state.scratch[2 * (i) + 1] + + output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i + } + } else { + work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + } + } + } + + namespace.complex = complex +}(FFT) diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/lib/node.erb.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/lib/node.erb.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,275 @@ +<%= $:.unshift('.'); require "#{File.dirname(__FILE__)}/../src/complex.rb"; File.read "#{File.dirname(__FILE__)}/../LICENSE" %> + +if (!FFT) { + var FFT = {} +} + +void function (namespace) { + "use strict" + + function butterfly2(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + + for (var i = 0; i < m; i++) { + <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %> + <%= load('s1', 'output', 'outputOffset', 'i + m', 'outputStride') %> + + <%= load('t1', 't', 0, 'i', 'fStride') %> + + <%= cmul('v1', 's1', 't1') %> + + <%= cadd('r0', 's0', 'v1') %> + <%= csub('r1', 's0', 'v1') %> + + <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %> + <%= store('r1', 'output', 'outputOffset', 'i + m', 'outputStride') %> + } + } + + function butterfly3(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m + var fStride1 = fStride, fStride2 = 2 * fStride + + var e = <%= imag('t', 0, 'm', 'fStride') %> + + for (var i = 0; i < m; i++) { + <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %> + + <%= load('s1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + <%= load('t1', 't', 0, 'i', 'fStride1') %> + <%= cmul('v1', 's1', 't1') %> + + <%= load('s2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + <%= load('t2', 't', 0, 'i', 'fStride2') %> + <%= cmul('v2', 's2', 't2') %> + + <%= cadd('i0', 'v1', 'v2') %> + + <%= cadd('r0', 's0', 'i0') %> + <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %> + + var i1_r = s0_r - i0_r * 0.5 + var i1_i = s0_i - i0_i * 0.5 + + var i2_r = (v1_r - v2_r) * e + var i2_i = (v1_i - v2_i) * e + + var r1_r = i1_r - i2_i + var r1_i = i1_i + i2_r + <%= store('r1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + + var r2_r = i1_r + i2_i + var r2_i = i1_i - i2_r + <%= store('r2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + } + } + + function butterfly4(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m, m3 = 3 * m + var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride + + for (var i = 0; i < m; i++) { + <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %> + + <%= load('s1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + <%= load('t1', 't', 0, 'i', 'fStride1') %> + <%= cmul('v1', 's1', 't1') %> + + <%= load('s2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + <%= load('t2', 't', 0, 'i', 'fStride2') %> + <%= cmul('v2', 's2', 't2') %> + + <%= load('s3', 'output', 'outputOffset', 'i + m3', 'outputStride') %> + <%= load('t3', 't', 0, 'i', 'fStride3') %> + <%= cmul('v3', 's3', 't3') %> + + <%= cadd('i0', 's0', 'v2') %> + <%= csub('i1', 's0', 'v2') %> + <%= cadd('i2', 'v1', 'v3') %> + <%= csub('i3', 'v1', 'v3') %> + + <%= cadd('r0', 'i0', 'i2') %> + + if (state.inverse) { + var r1_r = i1_r - i3_i + var r1_i = i1_i + i3_r + } else { + var r1_r = i1_r + i3_i + var r1_i = i1_i - i3_r + } + + <%= csub('r2', 'i0', 'i2') %> + + if (state.inverse) { + var r3_r = i1_r + i3_i + var r3_i = i1_i - i3_r + } else { + var r3_r = i1_r - i3_i + var r3_i = i1_i + i3_r + } + + <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %> + <%= store('r1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + <%= store('r2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + <%= store('r3', 'output', 'outputOffset', 'i + m3', 'outputStride') %> + } + } + + function butterfly(output, outputOffset, outputStride, fStride, state, m, p) { + var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p) + + for (var u = 0; u < m; u++) { + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + <%= load('x0', 'output', 'outputOffset', 'k', 'outputStride') %> + <%= store('x0', 'scratch', 'q1') %> + } + + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + var tOffset = 0 + + <%= load('x0', 'scratch', 0) %> + <%= store('x0', 'output', 'outputOffset', 'k', 'outputStride') %> + + for (var q = 1; q < p; q++) { + tOffset = (tOffset + fStride * k) % n + + <%= load('s0', 'output', 'outputOffset', 'k', 'outputStride') %> + + <%= load('s1', 'scratch', 'q') %> + <%= load('t1', 't', 'tOffset') %> + <%= cmul('v1', 's1', 't1') %> + + <%= cadd('r0', 's0', 'v1') %> + <%= store('r0', 'output', 'outputOffset', 'k', 'outputStride') %> + } + } + } + } + + function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { + var p = factors.shift() + var m = factors.shift() + + if (m == 1) { + for (var i = 0; i < p * m; i++) { + <%= load('x0', 'f', 'fOffset', 'i', 'fStride * inputStride') %> + <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %> + } + } else { + for (var i = 0; i < p; i++) { + work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) + } + } + + switch (p) { + case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break + case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break + case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break + default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break + } + } + + var complex = function (n, inverse) { + if (arguments.length < 2) { + throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'") + } + + var n = ~~n, inverse = !!inverse + + if (n < 1) { + throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") + } + + var state = { + n: n, + inverse: inverse, + + factors: [], + twiddle: new Float64Array(2 * n), + scratch: new Float64Array(2 * n) + } + + var t = state.twiddle, theta = 2 * Math.PI / n + + for (var i = 0; i < n; i++) { + if (inverse) { + var phase = theta * i + } else { + var phase = -theta * i + } + + <%= real('t', 'i') %> = Math.cos(phase) + <%= imag('t', 'i') %> = Math.sin(phase) + } + + var p = 4, v = Math.floor(Math.sqrt(n)) + + while (n > 1) { + while (n % p) { + switch (p) { + case 4: p = 2; break + case 2: p = 3; break + default: p += 2; break + } + + if (p > v) { + p = n + } + } + + n /= p + + state.factors.push(p) + state.factors.push(n) + } + + this.state = state + } + + complex.prototype.simple = function (output, input, t) { + this.process(output, 0, 1, input, 0, 1, t) + } + + complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) { + var outputStride = ~~outputStride, inputStride = ~~inputStride + + var type = t == 'real' ? t : 'complex' + + if (outputStride < 1) { + throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") + } + + if (inputStride < 1) { + throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") + } + + if (type == 'real') { + for (var i = 0; i < this.state.n; i++) { + var x0_r = input[inputOffset + inputStride * i] + var x0_i = 0.0 + + <%= store('x0', 'this.state.scratch', 'i') %> + } + + work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state) + } else { + if (input == output) { + work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + + for (var i = 0; i < this.state.n; i++) { + <%= load('x0', 'this.state.scratch', 'i') %> + + <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %> + } + } else { + work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + } + } + } + + namespace.complex = complex +}(FFT) + +module.exports = FFT \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/lib/node.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/lib/node.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,294 @@ +/* Copyright (c) 2012, Jens Nockert , Jussi Kalliokoski + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +if (!FFT) { + var FFT = {} +} + +void function (namespace) { + "use strict" + + function butterfly2(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + + for (var i = 0; i < m; i++) { + var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] + var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] + + var t1_r = t[2 * ((0) + (fStride) * (i))], t1_i = t[2 * ((0) + (fStride) * (i)) + 1] + + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i + var r1_r = s0_r - v1_r, r1_i = s0_i - v1_i + + output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i + output[2 * ((outputOffset) + (outputStride) * (i + m))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] = r1_i + } + } + + function butterfly3(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m + var fStride1 = fStride, fStride2 = 2 * fStride + + var e = t[2 * ((0) + (fStride) * (m)) + 1] + + for (var i = 0; i < m; i++) { + var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] + + var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] + var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] + var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1] + var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r + + var i0_r = v1_r + v2_r, i0_i = v1_i + v2_i + + var r0_r = s0_r + i0_r, r0_i = s0_i + i0_i + output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i + + var i1_r = s0_r - i0_r * 0.5 + var i1_i = s0_i - i0_i * 0.5 + + var i2_r = (v1_r - v2_r) * e + var i2_i = (v1_i - v2_i) * e + + var r1_r = i1_r - i2_i + var r1_i = i1_i + i2_r + output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i + + var r2_r = i1_r + i2_i + var r2_i = i1_i - i2_r + output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i + } + } + + function butterfly4(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m, m3 = 3 * m + var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride + + for (var i = 0; i < m; i++) { + var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] + + var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] + var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] + var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1] + var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r + + var s3_r = output[2 * ((outputOffset) + (outputStride) * (i + m3))], s3_i = output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] + var t3_r = t[2 * ((0) + (fStride3) * (i))], t3_i = t[2 * ((0) + (fStride3) * (i)) + 1] + var v3_r = s3_r * t3_r - s3_i * t3_i, v3_i = s3_r * t3_i + s3_i * t3_r + + var i0_r = s0_r + v2_r, i0_i = s0_i + v2_i + var i1_r = s0_r - v2_r, i1_i = s0_i - v2_i + var i2_r = v1_r + v3_r, i2_i = v1_i + v3_i + var i3_r = v1_r - v3_r, i3_i = v1_i - v3_i + + var r0_r = i0_r + i2_r, r0_i = i0_i + i2_i + + if (state.inverse) { + var r1_r = i1_r - i3_i + var r1_i = i1_i + i3_r + } else { + var r1_r = i1_r + i3_i + var r1_i = i1_i - i3_r + } + + var r2_r = i0_r - i2_r, r2_i = i0_i - i2_i + + if (state.inverse) { + var r3_r = i1_r + i3_i + var r3_i = i1_i - i3_r + } else { + var r3_r = i1_r - i3_i + var r3_i = i1_i + i3_r + } + + output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i + output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i + output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i + output[2 * ((outputOffset) + (outputStride) * (i + m3))] = r3_r, output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] = r3_i + } + } + + function butterfly(output, outputOffset, outputStride, fStride, state, m, p) { + var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p) + + for (var u = 0; u < m; u++) { + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + var x0_r = output[2 * ((outputOffset) + (outputStride) * (k))], x0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1] + scratch[2 * (q1)] = x0_r, scratch[2 * (q1) + 1] = x0_i + } + + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + var tOffset = 0 + + var x0_r = scratch[2 * (0)], x0_i = scratch[2 * (0) + 1] + output[2 * ((outputOffset) + (outputStride) * (k))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = x0_i + + for (var q = 1; q < p; q++) { + tOffset = (tOffset + fStride * k) % n + + var s0_r = output[2 * ((outputOffset) + (outputStride) * (k))], s0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1] + + var s1_r = scratch[2 * (q)], s1_i = scratch[2 * (q) + 1] + var t1_r = t[2 * (tOffset)], t1_i = t[2 * (tOffset) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i + output[2 * ((outputOffset) + (outputStride) * (k))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = r0_i + } + } + } + } + + function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { + var p = factors.shift() + var m = factors.shift() + + if (m == 1) { + for (var i = 0; i < p * m; i++) { + var x0_r = f[2 * ((fOffset) + (fStride * inputStride) * (i))], x0_i = f[2 * ((fOffset) + (fStride * inputStride) * (i)) + 1] + output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i + } + } else { + for (var i = 0; i < p; i++) { + work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) + } + } + + switch (p) { + case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break + case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break + case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break + default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break + } + } + + var complex = function (n, inverse) { + if (arguments.length < 2) { + throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'") + } + + var n = ~~n, inverse = !!inverse + + if (n < 1) { + throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") + } + + var state = { + n: n, + inverse: inverse, + + factors: [], + twiddle: new Float64Array(2 * n), + scratch: new Float64Array(2 * n) + } + + var t = state.twiddle, theta = 2 * Math.PI / n + + for (var i = 0; i < n; i++) { + if (inverse) { + var phase = theta * i + } else { + var phase = -theta * i + } + + t[2 * (i)] = Math.cos(phase) + t[2 * (i) + 1] = Math.sin(phase) + } + + var p = 4, v = Math.floor(Math.sqrt(n)) + + while (n > 1) { + while (n % p) { + switch (p) { + case 4: p = 2; break + case 2: p = 3; break + default: p += 2; break + } + + if (p > v) { + p = n + } + } + + n /= p + + state.factors.push(p) + state.factors.push(n) + } + + this.state = state + } + + complex.prototype.simple = function (output, input, t) { + this.process(output, 0, 1, input, 0, 1, t) + } + + complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) { + var outputStride = ~~outputStride, inputStride = ~~inputStride + + var type = t == 'real' ? t : 'complex' + + if (outputStride < 1) { + throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") + } + + if (inputStride < 1) { + throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") + } + + if (type == 'real') { + for (var i = 0; i < this.state.n; i++) { + var x0_r = input[inputOffset + inputStride * i] + var x0_i = 0.0 + + this.state.scratch[2 * (i)] = x0_r, this.state.scratch[2 * (i) + 1] = x0_i + } + + work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state) + } else { + if (input == output) { + work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + + for (var i = 0; i < this.state.n; i++) { + var x0_r = this.state.scratch[2 * (i)], x0_i = this.state.scratch[2 * (i) + 1] + + output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i + } + } else { + work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + } + } + } + + namespace.complex = complex +}(FFT) + +module.exports = FFT \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/lib/real.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/lib/real.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,300 @@ +/* Copyright (c) 2012, Jens Nockert , Jussi Kalliokoski + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +if (!FFT) { + var FFT = {} +} + +void function (namespace) { + "use strict" + + function forwardButterfly2(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride) { + var m = n / 2, q = n / product, old = product / 2 + + for (var i = 0; i < q; i++) { + var a0 = old * i + var a1 = a0 + m + + var s0 = input[inputOffset + inputStride * a0] + var s1 = input[inputOffset + inputStride * a1] + + var r0 = s0 + s1 + var r1 = s0 - s1 + + var a0 = product * i + var a1 = a0 + product - 1 + + output[outputOffset + outputStride * a0] = r0 + output[outputOffset + outputStride * a1] = r1 + } + + if (old == 1) { return } + + for (var i = 0; i < old / 2; i++) { + var t1_r = twiddle[2 * ((-1) + (i))], t1_i = twiddle[2 * ((-1) + (i)) + 1] + + for (var j = 0; j < q; j++) { + var a0 = j * old + 2 * i - 1 + var a1 = a0 + m + + var s0_r = input[2 * ((inputOffset) + (inputStride) * (a0))], s0_i = input[2 * ((inputOffset) + (inputStride) * (a0)) + 1] + + var s1_r = input[2 * ((inputOffset) + (inputStride) * (a1))], s1_i = input[2 * ((inputOffset) + (inputStride) * (a1)) + 1] + var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r + + var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i + var r1_r = s0_r - v1_r, r1_i = s0_i - v1_i; r1_i = -r1_i + + var a0 = j * product + 2 * i - 1 + var a1 = (j - 1) * product - 2 * i - 1 + + output[2 * ((outputOffset) + (outputStride) * (a0))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (a0)) + 1] = r0_i + output[2 * ((outputOffset) + (outputStride) * (a1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (a1)) + 1] = r1_i + } + } + + if (old % 2 == 1) { return } + + for (var i = 0; i < q; i++) { + var a0 = (i + 1) * old - 1 + var a1 = a0 + m + + var r0_r = input[2 * ((inputOffset) + (inputStride) * (a0))] + var r1_i = -input[2 * ((inputOffset) + (inputStride) * (a1))] + + var a0 = i * product + old - 1 + + output[2 * ((outputOffset) + (outputStride) * (a0))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (a0)) + 1] = r0_i + } + } + + function backwardButterfly2(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride) { + var m = n / 2, q = n / product, old = product / 2 + + for (var i = 0; i < q; i++) { + var a0 = (2 * i) * q + var a1 = (2 * i + 2) * q - 1 + + var s0 = input[inputOffset + inputStride * a0] + var s1 = input[inputOffset + inputStride * a1] + + var r0 = s0 + s1 + var r1 = s0 - s1 + + var a0 = q * i + var a1 = q * i + m + + output[outputOffset + outputStride * a0] = r0 + output[outputOffset + outputStride * a1] = r1 + } + + if (q == 1) { return } + + for (var i = 0; i < q / 2; i++) { + var t1_r = twiddle[2 * ((-1) + (i))], t1_i = twiddle[2 * ((-1) + (i)) + 1] + + for (var j = 0; j < old; j++) { + var a0 = 2 * j * q + 2 * i - 1 + var a1 = 2 * (j + 1) * q - 2 * i - 1 + + var s0_r = input[2 * ((inputOffset) + (inputStride) * (a0))], s0_i = input[2 * ((inputOffset) + (inputStride) * (a0)) + 1] + var s1_r = input[2 * ((inputOffset) + (inputStride) * (a1))], s1_i = input[2 * ((inputOffset) + (inputStride) * (a1)) + 1] + + var r0_r = s0_r + s1_r + var r0_i = s0_i - s1_i + + var v1_r = s0_r - s1_r + var v1_i = s0_i + s1_i + + var r1_r = v1_r * t1_r - v1_i * t1_i, r1_i = v1_r * t1_i + v1_i * t1_r + + var a0 = j * q + 2 * i - 1 + var a1 = a0 + m + + output[2 * ((outputOffset) + (outputStride) * (a0))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (a0)) + 1] = r0_i + output[2 * ((outputOffset) + (outputStride) * (a1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (a1)) + 1] = r1_i + } + } + + if (q % 2 == 1) { return } + + for (var i = 0; i < q; i++) { + var a0 = 2 * (i + 1) * q - 1 + + var r0_r = input[2 * ((inputOffset) + (inputStride) * (a0))], r0_i = input[2 * ((inputOffset) + (inputStride) * (a0)) + 1] + + input[2 * ((inputOffset) + (inputStride) * (a0))] = 2 * r0_r + input[2 * ((inputOffset) + (inputStride) * (a1)) + 1] = -2 * r0_i + } + } + + function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { + var p = factors.shift() + var m = factors.shift() + + if (m == 1) { + for (var i = 0; i < p * m; i++) { + var x0_r = f[2 * ((fOffset) + (fStride * inputStride) * (i))], x0_i = f[2 * ((fOffset) + (fStride * inputStride) * (i)) + 1] + output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i + } + } else { + for (var i = 0; i < p; i++) { + work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) + } + } + + switch (p) { + case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break + case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break + case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break + default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break + } + } + + var real = function (n, inverse) { + var n = ~~n, inverse = !!inverse + + if (n < 1) { + throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") + } + + var state = { + n: n, + inverse: inverse, + + factors: [], + twiddle: [], + scratch: new Float64Array(n) + } + + var t = new Float64Array(n) + + var p = 4, v = Math.floor(Math.sqrt(n)) + + while (n > 1) { + while (n % p) { + switch (p) { + case 4: p = 2; break + case 2: p = 3; break + default: p += 2; break + } + + if (p > v) { + p = n + } + } + + n /= p + + state.factors.push(p) + } + + var theta = 2 * Math.PI / n, product = 1, twiddle = new Float64Array(n) + + for (var i = 0, t = 0; i < state.factors.length; i++) { + var phase = theta * i, factor = state.factors[i] + + var old = product, product = product * factor, q = n / product + + state.twiddle.push(new Float64Array(twiddle, t)) + + if (inverse) { + var counter = q, multiplier = old + } else { + var counter = old, multiplier = q + } + + for (var j = 1; j < factor; j++) { + var m = 0 + + for (var k = 1; k < counter / 2; k++, t++) { + m = (m + j * multiplier) % n + + var phase = theta * m + + t[2 * (i)] = Math.cos(phase) + t[2 * (i) + 1] = Math.sin(phase) + } + } + } + + this.state = state + } + + real.prototype.process = function(output, outputStride, input, inputStride) { + var outputStride = ~~outputStride, inputStride = ~~inputStride + + if (outputStride < 1) { + throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") + } + + if (inputStride < 1) { + throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") + } + + var product = 1, state = 0, inverse = this.state.inverse + + var n = this.state.n, factors = this.state.factors + var twiddle = this.state.twiddle, scratch = this.state.scratch + + for (var i = 0; i < factors.length; i++) { + var factor = factors[i], old = product, product = product * factor + + var q = n / product, fStride = Math.ceil(old / 2) - 1 + + if (state == 0) { + var inBuffer = input, inStride = inputStride + + if (this.state.factors.length % 2 == 0) { + var outBuffer = scratch, outStride = 1, state = 1 + } else { + var outBuffer = output, outStride = outputStride, state = 2 + } + } else if (state == 1) { + var inBuffer = scratch, inStride = 1, outBuffer = output, outStride = outputStride, state = 2 + } else if (state == 2) { + var inBuffer = output, inStride = outputStride, outBuffer = scratch, outStride = 1, state = 1 + } else { + throw new RangeError("state somehow is not in the range (0 .. 2)") + } + + if (inverse) { + switch (factor) { + case 2: backwardButterfly2(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 3: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 4: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 5: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + default: backwardButterfly(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + } + } else { + switch (factor) { + case 2: forwardButterfly2(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 3: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 4: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 5: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + default: forwardButterfly(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + } + } + } + } + + namespace.real = real +}(FFT) diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/package.json --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/package.json Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,20 @@ +{ + "name": "fft", + "authors": ["Jens Nockert", "Jussi Kalliokoski"], + "version": "0.2.1", + "description": "A Fast Fourier Transform library for JS.", + "repository": { + "type": "git", + "url": "http://github.com/JensNockert/fft.js" + }, + "engines": { + "node": ">=0.6.0-0" + }, + "licenses": [{ + "type": "FreeBSD" + }], + "bugs": { + "url": "http://github.com/JensNockert/fft.js/issues" + }, + "main": "./lib/node.js" +} diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/src/complex.erb.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/src/complex.erb.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,273 @@ +<%= $:.unshift('.'); require "#{File.dirname(__FILE__)}/../src/complex.rb"; File.read "#{File.dirname(__FILE__)}/../LICENSE" %> + +if (!FFT) { + var FFT = {} +} + +void function (namespace) { + "use strict" + + function butterfly2(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + + for (var i = 0; i < m; i++) { + <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %> + <%= load('s1', 'output', 'outputOffset', 'i + m', 'outputStride') %> + + <%= load('t1', 't', 0, 'i', 'fStride') %> + + <%= cmul('v1', 's1', 't1') %> + + <%= cadd('r0', 's0', 'v1') %> + <%= csub('r1', 's0', 'v1') %> + + <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %> + <%= store('r1', 'output', 'outputOffset', 'i + m', 'outputStride') %> + } + } + + function butterfly3(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m + var fStride1 = fStride, fStride2 = 2 * fStride + + var e = <%= imag('t', 0, 'm', 'fStride') %> + + for (var i = 0; i < m; i++) { + <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %> + + <%= load('s1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + <%= load('t1', 't', 0, 'i', 'fStride1') %> + <%= cmul('v1', 's1', 't1') %> + + <%= load('s2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + <%= load('t2', 't', 0, 'i', 'fStride2') %> + <%= cmul('v2', 's2', 't2') %> + + <%= cadd('i0', 'v1', 'v2') %> + + <%= cadd('r0', 's0', 'i0') %> + <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %> + + var i1_r = s0_r - i0_r * 0.5 + var i1_i = s0_i - i0_i * 0.5 + + var i2_r = (v1_r - v2_r) * e + var i2_i = (v1_i - v2_i) * e + + var r1_r = i1_r - i2_i + var r1_i = i1_i + i2_r + <%= store('r1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + + var r2_r = i1_r + i2_i + var r2_i = i1_i - i2_r + <%= store('r2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + } + } + + function butterfly4(output, outputOffset, outputStride, fStride, state, m) { + var t = state.twiddle + var m1 = m, m2 = 2 * m, m3 = 3 * m + var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride + + for (var i = 0; i < m; i++) { + <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %> + + <%= load('s1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + <%= load('t1', 't', 0, 'i', 'fStride1') %> + <%= cmul('v1', 's1', 't1') %> + + <%= load('s2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + <%= load('t2', 't', 0, 'i', 'fStride2') %> + <%= cmul('v2', 's2', 't2') %> + + <%= load('s3', 'output', 'outputOffset', 'i + m3', 'outputStride') %> + <%= load('t3', 't', 0, 'i', 'fStride3') %> + <%= cmul('v3', 's3', 't3') %> + + <%= cadd('i0', 's0', 'v2') %> + <%= csub('i1', 's0', 'v2') %> + <%= cadd('i2', 'v1', 'v3') %> + <%= csub('i3', 'v1', 'v3') %> + + <%= cadd('r0', 'i0', 'i2') %> + + if (state.inverse) { + var r1_r = i1_r - i3_i + var r1_i = i1_i + i3_r + } else { + var r1_r = i1_r + i3_i + var r1_i = i1_i - i3_r + } + + <%= csub('r2', 'i0', 'i2') %> + + if (state.inverse) { + var r3_r = i1_r + i3_i + var r3_i = i1_i - i3_r + } else { + var r3_r = i1_r - i3_i + var r3_i = i1_i + i3_r + } + + <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %> + <%= store('r1', 'output', 'outputOffset', 'i + m1', 'outputStride') %> + <%= store('r2', 'output', 'outputOffset', 'i + m2', 'outputStride') %> + <%= store('r3', 'output', 'outputOffset', 'i + m3', 'outputStride') %> + } + } + + function butterfly(output, outputOffset, outputStride, fStride, state, m, p) { + var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p) + + for (var u = 0; u < m; u++) { + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + <%= load('x0', 'output', 'outputOffset', 'k', 'outputStride') %> + <%= store('x0', 'scratch', 'q1') %> + } + + for (var q1 = 0, k = u; q1 < p; q1++, k += m) { + var tOffset = 0 + + <%= load('x0', 'scratch', 0) %> + <%= store('x0', 'output', 'outputOffset', 'k', 'outputStride') %> + + for (var q = 1; q < p; q++) { + tOffset = (tOffset + fStride * k) % n + + <%= load('s0', 'output', 'outputOffset', 'k', 'outputStride') %> + + <%= load('s1', 'scratch', 'q') %> + <%= load('t1', 't', 'tOffset') %> + <%= cmul('v1', 's1', 't1') %> + + <%= cadd('r0', 's0', 'v1') %> + <%= store('r0', 'output', 'outputOffset', 'k', 'outputStride') %> + } + } + } + } + + function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { + var p = factors.shift() + var m = factors.shift() + + if (m == 1) { + for (var i = 0; i < p * m; i++) { + <%= load('x0', 'f', 'fOffset', 'i', 'fStride * inputStride') %> + <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %> + } + } else { + for (var i = 0; i < p; i++) { + work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) + } + } + + switch (p) { + case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break + case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break + case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break + default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break + } + } + + var complex = function (n, inverse) { + if (arguments.length < 2) { + throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'") + } + + var n = ~~n, inverse = !!inverse + + if (n < 1) { + throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") + } + + var state = { + n: n, + inverse: inverse, + + factors: [], + twiddle: new Float64Array(2 * n), + scratch: new Float64Array(2 * n) + } + + var t = state.twiddle, theta = 2 * Math.PI / n + + for (var i = 0; i < n; i++) { + if (inverse) { + var phase = theta * i + } else { + var phase = -theta * i + } + + <%= real('t', 'i') %> = Math.cos(phase) + <%= imag('t', 'i') %> = Math.sin(phase) + } + + var p = 4, v = Math.floor(Math.sqrt(n)) + + while (n > 1) { + while (n % p) { + switch (p) { + case 4: p = 2; break + case 2: p = 3; break + default: p += 2; break + } + + if (p > v) { + p = n + } + } + + n /= p + + state.factors.push(p) + state.factors.push(n) + } + + this.state = state + } + + complex.prototype.simple = function (output, input, t) { + this.process(output, 0, 1, input, 0, 1, t) + } + + complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) { + var outputStride = ~~outputStride, inputStride = ~~inputStride + + var type = t == 'real' ? t : 'complex' + + if (outputStride < 1) { + throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") + } + + if (inputStride < 1) { + throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") + } + + if (type == 'real') { + for (var i = 0; i < this.state.n; i++) { + var x0_r = input[inputOffset + inputStride * i] + var x0_i = 0.0 + + <%= store('x0', 'this.state.scratch', 'i') %> + } + + work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state) + } else { + if (input == output) { + work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + + for (var i = 0; i < this.state.n; i++) { + <%= load('x0', 'this.state.scratch', 'i') %> + + <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %> + } + } else { + work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) + } + } + } + + namespace.complex = complex +}(FFT) diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/src/complex.rb --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/src/complex.rb Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,43 @@ +def real(x, offset = 0, i = nil, stride = nil) + if stride + "#{x}[2 * ((#{offset}) + (#{stride}) * (#{i}))]" + elsif i + "#{x}[2 * ((#{offset}) + (#{i}))]" + elsif offset + "#{x}[2 * (#{offset})]" + else + "#{x}[0]" + end +end + +def imag(x, offset = 0, i = nil, stride = nil) + if stride + "#{x}[2 * ((#{offset}) + (#{stride}) * (#{i})) + 1]" + elsif i + "#{x}[2 * ((#{offset}) + (#{i})) + 1]" + elsif offset + "#{x}[2 * (#{offset}) + 1]" + else + "#{x}[1]" + end +end + +def load(value, x, offset = 0, i = nil, stride = nil) + "var #{value}_r = #{real(x, offset, i, stride)}, #{value}_i = #{imag(x, offset, i, stride)}" +end + +def store(value, x, offset = 0, i = nil, stride = nil) + "#{real(x, offset, i, stride)} = #{value}_r, #{imag(x, offset, i, stride)} = #{value}_i" +end + +def cadd(result, a, b) + "var #{result}_r = #{a}_r + #{b}_r, #{result}_i = #{a}_i + #{b}_i" +end + +def csub(result, a, b) + "var #{result}_r = #{a}_r - #{b}_r, #{result}_i = #{a}_i - #{b}_i" +end + +def cmul(result, a, b) + "var #{result}_r = #{a}_r * #{b}_r - #{a}_i * #{b}_i, #{result}_i = #{a}_r * #{b}_i + #{a}_i * #{b}_r" +end \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/src/node.erb.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/src/node.erb.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,2 @@ +<%= File.read "#{File.dirname(__FILE__)}/../src/complex.erb.js" %> +module.exports = FFT \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/src/real.erb.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/src/real.erb.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,541 @@ +<%= $:.unshift('.'); require "#{File.dirname(__FILE__)}/../src/complex.rb"; File.read "#{File.dirname(__FILE__)}/../LICENSE" %> + +if (!FFT) { + var FFT = {} +} + +void function (namespace) { + "use strict" + + function forwardButterfly2(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride) { + var m = n / 2, q = n / product, old = product / 2 + + for (var i = 0; i < q; i++) { + var a0 = old * i + var a1 = a0 + m + + var s0 = input[inputOffset + inputStride * a0] + var s1 = input[inputOffset + inputStride * a1] + + var r0 = s0 + s1 + var r1 = s0 - s1 + + var a0 = product * i + var a1 = a0 + product - 1 + + output[outputOffset + outputStride * a0] = r0 + output[outputOffset + outputStride * a1] = r1 + } + + if (old == 1) { return } + + for (var i = 0; i < old / 2; i++) { + <%= load('t1', 'twiddle', -1, 'i') %> + + for (var j = 0; j < q; j++) { + var a0 = j * old + 2 * i - 1 + var a1 = a0 + m + + <%= load('s0', 'input', 'inputOffset', 'a0', 'inputStride') %> + + <%= load('s1', 'input', 'inputOffset', 'a1', 'inputStride') %> + <%= cmul('v1', 's1', 't1') %> + + <%= cadd('r0', 's0', 'v1') %> + <%= csub('r1', 's0', 'v1') %>; r1_i = -r1_i + + var a0 = j * product + 2 * i - 1 + var a1 = (j - 1) * product - 2 * i - 1 + + <%= store('r0', 'output', 'outputOffset', 'a0', 'outputStride') %> + <%= store('r1', 'output', 'outputOffset', 'a1', 'outputStride') %> + } + } + + if (old % 2 == 1) { return } + + for (var i = 0; i < q; i++) { + var a0 = (i + 1) * old - 1 + var a1 = a0 + m + + var r0_r = <%= real('input', 'inputOffset', 'a0', 'inputStride') %> + var r1_i = -<%= real('input', 'inputOffset', 'a1', 'inputStride') %> + + var a0 = i * product + old - 1 + + <%= store('r0', 'output', 'outputOffset', 'a0', 'outputStride') %> + } + } + + function forwardButterfly(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride, factor) { + var m = n / 2, q = n / product, old = product / 2 + + var theta = 2.0 * Math.PI / factor + + var theta_r = Math.cos(theta) + var theta_i = -Math.sin(theta) + + for (var i = 0; i < q; i++) { + var d_r = 1.0 + var d_i = 0.0 + + for (var j = 0; j <= Math.floor(factor / 2); j++) { + var sum_r = 0.0 + var sum_i = 0.0 + + var w_r = 1.0 + var w_i = 0.0 + + if (j > 0) { + <%= cmul('t0', 'd', 'theta') %> + + d_r = t0_r + d_i = t0_i + } + + for (var k = 0; k < factor; k++) { + var z = <%= real('input', 'inputOffset', 'i * old + k * m', 'inputStride') %> + + if (k > 0) { + <%= cmul('t0', 'd', 'w') %> + + d_r = t0_r + d_i = t0_i + } + + /* TODO: Use Kahan summation..? */ + s_r += w_r * z + s_i += w_i * z + } + + if (j == 0) { + var a0 = product * i + + output[outputOffset + outputStride * a0] = sum_r + } else if (j < factor / 2) { + var a0 = product * i + + output[outputOffset + outputStride * a0] = sum_r + } else if (j == factor / 2) { + + } + } + for (e1 = 0; e1 <= factor - e1; e1++) + { + if (e1 == 0) + { + const size_t to0 = product * k1; + VECTOR(out,ostride,to0) = sum_real; + } + else if (e1 < factor - e1) + { + const size_t to0 = k1 * product + 2 * e1 * product_1 - 1; + VECTOR(out,ostride,to0) = sum_real; + VECTOR(out,ostride,to0 + 1) = sum_imag; + } + else if (e1 == factor - e1) + { + const size_t to0 = k1 * product + 2 * e1 * product_1 - 1; + VECTOR(out,ostride,to0) = sum_real; + } + + } + + } + + if (old == 1) { return } + + for (var i = 0; i < old / 2; i++) { + + } + + if (old % 2 == 1) { return } + + var t_arg = Math.PI / factor + + var t_r = cos(t_arg) + var t_i = -sin(t_arg) + + for (var i = 0; i < q; i++) { + + } + } + + if (product_1 == 1) + return; + + for (k = 1; k < (product_1 + 1) / 2; k++) + { + for (k1 = 0; k1 < q; k1++) + { + + ATOMIC dw_real = 1.0, dw_imag = 0.0; + + for (e1 = 0; e1 < factor; e1++) + { + ATOMIC sum_real = 0.0, sum_imag = 0.0; + + ATOMIC w_real = 1.0, w_imag = 0.0; + + if (e1 > 0) + { + const ATOMIC tmp_real = dw_real * cos_d_theta + dw_imag * sin_d_theta; + const ATOMIC tmp_imag = -dw_real * sin_d_theta + dw_imag * cos_d_theta; + dw_real = tmp_real; + dw_imag = tmp_imag; + } + + for (e2 = 0; e2 < factor; e2++) + { + + int tskip = (product_1 + 1) / 2 - 1; + const size_t from0 = k1 * product_1 + 2 * k + e2 * m - 1; + ATOMIC tw_real, tw_imag; + ATOMIC z_real, z_imag; + + if (e2 == 0) + { + tw_real = 1.0; + tw_imag = 0.0; + } + else + { + const size_t t_index = (k - 1) + (e2 - 1) * tskip; + tw_real = GSL_REAL(twiddle[t_index]); + tw_imag = -GSL_IMAG(twiddle[t_index]); + } + + { + const ATOMIC f0_real = VECTOR(in,istride,from0); + const ATOMIC f0_imag = VECTOR(in,istride,from0 + 1); + + z_real = tw_real * f0_real - tw_imag * f0_imag; + z_imag = tw_real * f0_imag + tw_imag * f0_real; + } + + if (e2 > 0) + { + const ATOMIC tmp_real = dw_real * w_real - dw_imag * w_imag; + const ATOMIC tmp_imag = dw_real * w_imag + dw_imag * w_real; + w_real = tmp_real; + w_imag = tmp_imag; + } + + sum_real += w_real * z_real - w_imag * z_imag; + sum_imag += w_real * z_imag + w_imag * z_real; + } + + if (e1 < factor - e1) + { + const size_t to0 = k1 * product - 1 + 2 * e1 * product_1 + 2 * k; + VECTOR(out,ostride,to0) = sum_real; + VECTOR(out,ostride,to0 + 1) = sum_imag; + } + else + { + const size_t to0 = k1 * product - 1 + 2 * (factor - e1) * product_1 - 2 * k; + VECTOR(out,ostride,to0) = sum_real; + VECTOR(out,ostride,to0 + 1) = -sum_imag; + } + + } + } + } + + + if (product_1 % 2 == 1) + return; + + { + double tw_arg = M_PI / ((double) factor); + ATOMIC cos_tw_arg = cos (tw_arg); + ATOMIC sin_tw_arg = -sin (tw_arg); + + for (k1 = 0; k1 < q; k1++) + { + ATOMIC dw_real = 1.0, dw_imag = 0.0; + + for (e1 = 0; e1 < factor; e1++) + { + ATOMIC z_real, z_imag; + + ATOMIC sum_real = 0.0; + ATOMIC sum_imag = 0.0; + + ATOMIC w_real = 1.0, w_imag = 0.0; + ATOMIC tw_real = 1.0, tw_imag = 0.0; + + if (e1 > 0) + { + ATOMIC t_real = dw_real * cos_d_theta + dw_imag * sin_d_theta; + ATOMIC t_imag = -dw_real * sin_d_theta + dw_imag * cos_d_theta; + dw_real = t_real; + dw_imag = t_imag; + } + + for (e2 = 0; e2 < factor; e2++) + { + + if (e2 > 0) + { + ATOMIC tmp_real = tw_real * cos_tw_arg - tw_imag * sin_tw_arg; + ATOMIC tmp_imag = tw_real * sin_tw_arg + tw_imag * cos_tw_arg; + tw_real = tmp_real; + tw_imag = tmp_imag; + } + + if (e2 > 0) + { + ATOMIC tmp_real = dw_real * w_real - dw_imag * w_imag; + ATOMIC tmp_imag = dw_real * w_imag + dw_imag * w_real; + w_real = tmp_real; + w_imag = tmp_imag; + } + + + { + const size_t from0 = k1 * product_1 + 2 * k + e2 * m - 1; + const ATOMIC f0_real = VECTOR(in,istride,from0); + z_real = tw_real * f0_real; + z_imag = tw_imag * f0_real; + } + + sum_real += w_real * z_real - w_imag * z_imag; + sum_imag += w_real * z_imag + w_imag * z_real; + } + + if (e1 + 1 < factor - e1) + { + const size_t to0 = k1 * product - 1 + 2 * e1 * product_1 + 2 * k; + VECTOR(out,ostride,to0) = sum_real; + VECTOR(out,ostride,to0 + 1) = sum_imag; + } + else if (e1 + 1 == factor - e1) + { + const size_t to0 = k1 * product - 1 + 2 * e1 * product_1 + 2 * k; + VECTOR(out,ostride,to0) = sum_real; + } + else + { + const size_t to0 = k1 * product - 1 + 2 * (factor - e1) * product_1 - 2 * k; + VECTOR(out,ostride,to0) = sum_real; + VECTOR(out,ostride,to0 + 1) = -sum_imag; + } + + } + } + } + return; + + function backwardButterfly2(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride) { + var m = n / 2, q = n / product, old = product / 2 + + for (var i = 0; i < q; i++) { + var a0 = (2 * i) * q + var a1 = (2 * i + 2) * q - 1 + + var s0 = input[inputOffset + inputStride * a0] + var s1 = input[inputOffset + inputStride * a1] + + var r0 = s0 + s1 + var r1 = s0 - s1 + + var a0 = q * i + var a1 = q * i + m + + output[outputOffset + outputStride * a0] = r0 + output[outputOffset + outputStride * a1] = r1 + } + + if (q == 1) { return } + + for (var i = 0; i < q / 2; i++) { + <%= load('t1', 'twiddle', -1, 'i') %> + + for (var j = 0; j < old; j++) { + var a0 = 2 * j * q + 2 * i - 1 + var a1 = 2 * (j + 1) * q - 2 * i - 1 + + <%= load('s0', 'input', 'inputOffset', 'a0', 'inputStride') %> + <%= load('s1', 'input', 'inputOffset', 'a1', 'inputStride') %> + + var r0_r = s0_r + s1_r + var r0_i = s0_i - s1_i + + var v1_r = s0_r - s1_r + var v1_i = s0_i + s1_i + + <%= cmul('r1', 'v1', 't1') %> + + var a0 = j * q + 2 * i - 1 + var a1 = a0 + m + + <%= store('r0', 'output', 'outputOffset', 'a0', 'outputStride') %> + <%= store('r1', 'output', 'outputOffset', 'a1', 'outputStride') %> + } + } + + if (q % 2 == 1) { return } + + for (var i = 0; i < q; i++) { + var a0 = 2 * (i + 1) * q - 1 + + <%= load('r0', 'input', 'inputOffset', 'a0', 'inputStride') %> + + <%= real('input', 'inputOffset', 'a0', 'inputStride') %> = 2 * r0_r + <%= imag('input', 'inputOffset', 'a1', 'inputStride') %> = -2 * r0_i + } + } + + function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { + var p = factors.shift() + var m = factors.shift() + + if (m == 1) { + for (var i = 0; i < p * m; i++) { + <%= load('x0', 'f', 'fOffset', 'i', 'fStride * inputStride') %> + <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %> + } + } else { + for (var i = 0; i < p; i++) { + work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) + } + } + + switch (p) { + case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break + case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break + case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break + default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break + } + } + + var real = function (n, inverse) { + var n = ~~n, inverse = !!inverse + + if (n < 1) { + throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") + } + + var state = { + n: n, + inverse: inverse, + + factors: [], + twiddle: [], + scratch: new Float64Array(n) + } + + var t = new Float64Array(n) + + var p = 4, v = Math.floor(Math.sqrt(n)) + + while (n > 1) { + while (n % p) { + switch (p) { + case 4: p = 2; break + case 2: p = 3; break + default: p += 2; break + } + + if (p > v) { + p = n + } + } + + n /= p + + state.factors.push(p) + } + + var theta = 2 * Math.PI / n, product = 1, twiddle = new Float64Array(n) + + for (var i = 0, t = 0; i < state.factors.length; i++) { + var phase = theta * i, factor = state.factors[i] + + var old = product, product = product * factor, q = n / product + + state.twiddle.push(new Float64Array(twiddle, t)) + + if (inverse) { + var counter = q, multiplier = old + } else { + var counter = old, multiplier = q + } + + for (var j = 1; j < factor; j++) { + var m = 0 + + for (var k = 1; k < counter / 2; k++, t++) { + m = (m + j * multiplier) % n + + var phase = theta * m + + <%= real('t', 'i') %> = Math.cos(phase) + <%= imag('t', 'i') %> = Math.sin(phase) + } + } + } + + this.state = state + } + + real.prototype.process = function(output, outputStride, input, inputStride) { + var outputStride = ~~outputStride, inputStride = ~~inputStride + + if (outputStride < 1) { + throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") + } + + if (inputStride < 1) { + throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") + } + + var product = 1, state = 0, inverse = this.state.inverse + + var n = this.state.n, factors = this.state.factors + var twiddle = this.state.twiddle, scratch = this.state.scratch + + for (var i = 0; i < factors.length; i++) { + var factor = factors[i], old = product, product = product * factor + + var q = n / product, fStride = Math.ceil(old / 2) - 1 + + if (state == 0) { + var inBuffer = input, inStride = inputStride + + if (this.state.factors.length % 2 == 0) { + var outBuffer = scratch, outStride = 1, state = 1 + } else { + var outBuffer = output, outStride = outputStride, state = 2 + } + } else if (state == 1) { + var inBuffer = scratch, inStride = 1, outBuffer = output, outStride = outputStride, state = 2 + } else if (state == 2) { + var inBuffer = output, inStride = outputStride, outBuffer = scratch, outStride = 1, state = 1 + } else { + throw new RangeError("state somehow is not in the range (0 .. 2)") + } + + if (inverse) { + switch (factor) { + case 2: backwardButterfly2(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 3: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 4: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 5: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + default: backwardButterfly(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride, factor); break + } + } else { + switch (factor) { + case 2: forwardButterfly2(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 3: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 4: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + case 5: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break + default: forwardButterfly(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride, factor); break + } + } + } + } + + namespace.real = real +}(FFT) diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/src/real.rb --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/src/real.rb Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,47 @@ +def load(offset = 0, i = nil, stride = nil) + "#{real(x, offset, i, stride)}" +end + +def creal(x, offset = 0, i = nil, stride = nil) + if stride + "#{x}[2 * ((#{offset}) + (#{stride}) * (#{i}))]" + elsif i + "#{x}[2 * ((#{offset}) + (#{i}))]" + elsif offset + "#{x}[2 * (#{offset})]" + else + "#{x}[0]" + end +end + +def cimag(x, offset = 0, i = nil, stride = nil) + if stride + "#{x}[2 * ((#{offset}) + (#{stride}) * (#{i})) + 1]" + elsif i + "#{x}[2 * ((#{offset}) + (#{i})) + 1]" + elsif offset + "#{x}[2 * (#{offset}) + 1]" + else + "#{x}[1]" + end +end + +def cload(value, x, offset = 0, i = nil, stride = nil) + "var #{value}_r = #{real(x, offset, i, stride)}, #{value}_i = #{imag(x, offset, i, stride)}" +end + +def cstore(value, x, offset = 0, i = nil, stride = nil) + "#{real(x, offset, i, stride)} = #{value}_r, #{imag(x, offset, i, stride)} = #{value}_i" +end + +def cadd(result, a, b) + "var #{result}_r = #{a}_r + #{b}_r, #{result}_i = #{a}_i + #{b}_i" +end + +def csub(result, a, b) + "var #{result}_r = #{a}_r - #{b}_r, #{result}_i = #{a}_i - #{b}_i" +end + +def cmul(result, a, b) + "var #{result}_r = #{a}_r * #{b}_r - #{a}_i * #{b}_i, #{result}_i = #{a}_r * #{b}_i + #{a}_i * #{b}_r" +end diff -r e705de983b67 -r 66f9fd5ac611 fft/fft.js/test.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/fft.js/test.html Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,171 @@ + + + + fft.js + + + + \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/index.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/index.html Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + +

Results

+ +

+ + + + + + + + + + + + + + + + + + + + +
ImplementationResultTime (first half)Time (second half)Rate (second half)
Nayuki
Nayuki (obj)
Nockert
Dntj
Cross
KissFFT
FFTW
+ +

Notes

+ +
    +
  • Nayuki: in-place single-precision complex-complex. Around 7kb.
    • +
  • Nayuki (obj): Nayuki with the sin/cos tables pre-calculated on object construction. Around 4kb.
    • +
  • Nockert: double-precision real-complex. Around 25kb.
    • +
  • Dntj: double-precision complex-complex. Forward + transform is scaled and I've scaled it back again here, which may + introduce rounding error. Around 10kb.
    • +
  • Cross: double-precision real-complex in C, compiled + with Emscripten. This is considered a slow implementation amongst + native code ones. Around 60kb.
    • +
  • KissFFT: single-precision real-complex in C, compiled + with Emscripten. A reasonably sophisticated implementation. Around + 70kb.
    • +
  • FFTW: single-precision real-complex in C, compiled with + Emscripten. GPL licensed. Around 3Mb.
    • +
+ + + diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/LICENSE --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/LICENSE Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,21 @@ +The MIT License + +Copyright (c) 2012 Nick Jones + +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. \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/README.md --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/README.md Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,45 @@ +# jsfft + +Small, efficient Javascript FFT implementation for node or the browser. + +## Usage + +JSFFT ships with a **complex_array** and a **fft** module. + +```javascript +var data = new complex_array.ComplexArray(512) +// Use the in-place mapper to populate the data. +data.map(function(value, i, n) { + value.real = (i > n/3 && i < 2*n/3) ? 1 : 0 +}) +``` + +Including the **fft** module attaches FFT methods to ComplexArray. FFT and +InvFFT perform in-place transforms on the underlying data: + +```javascript +var frequencies = data.FFT() +// Implement a low-pass filter using the in-place mapper. +frequencies.map(function(frequency, i, n) { + if (i > n/5 && i < 4*n/5) { + frequency.real = 0 + frequency.imag = 0 + } +}) +``` + +Alternatively, frequency-space filters can be implemented via the frequencyMap: + +```javascript +var filtered = data.frequencyMap(function(frequency, i, n) { + if (i > n/5 && i < 4*n/5) { + frequency.real = 0 + frequency.imag = 0 + } +}) +``` + +## Other Implementations + +[DSP](https://github.com/corbanbrook/dsp.js) is a full featured Digital Signal +Processing library in JS which includes a JS FFT implementation. diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/example/example.css --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/example/example.css Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,34 @@ +body { + font-family: sans-serif; +} + +canvas { + background-color: #eee; +} + +div.panel { + float: left; + width: 312px; + height: 312px; + border: 1px solid #ccc; + margin: 1%; +} + +div.panel div.example { + width: 256px; + height: 128px; + left: 28px; + position: relative; +} + +div.explanation { + left: 0; + right: 0; + height: 50%; + padding: 20px 10px 0 10px; + overflow: scroll; +} + +pre, code { + font-size: 85%; +} diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/example/index.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/example/index.html Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,126 @@ + + + + + + + + + + +
+
+
+ Sample data + 
+var data = new ComplexArray(128)
+data.map(function(value, i, n) {
+  value.real = (i > n/3 && i < 2*n/3) ? 1 : 0
+})
+
+
+
+
+
+
+ Transform (in place):
+ + data.FFT() + +
+
+
+
+
+ Simple low pass filter: + 
+data.map(function(freq, i, n) {
+  if (i > n/5 && i < 4*n/5) {
+    freq.real = 0
+    freq.imag = 0
+  }
+})
+
+
+
+
+
+
+ Transform back:
+ + data.InvFFT() + +
+
+
+
+
+ ... or all in one step + 
+data.frequencyMap(function(freq, i, n) {
+  if (i > n/5 && i < 4*n/5) {
+    freq.real = 0
+    freq.imag = 0
+  }
+})
+
+
+
+ + + diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/lib/complex_array.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/lib/complex_array.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,116 @@ +'use strict'; + +!function(exports, undefined) { + + var + // If the typed array is unspecified, use this. + DefaultArrayType = Float32Array, + // Simple math functions we need. + sqrt = Math.sqrt, + sqr = function(number) {return Math.pow(number, 2)}, + // Internal convenience copies of the exported functions + isComplexArray, + ComplexArray + + exports.isComplexArray = isComplexArray = function(obj) { + return obj !== undefined && + obj.hasOwnProperty !== undefined && + obj.hasOwnProperty('real') && + obj.hasOwnProperty('imag') + } + + exports.ComplexArray = ComplexArray = function(other, opt_array_type){ + if (isComplexArray(other)) { + // Copy constuctor. + this.ArrayType = other.ArrayType + this.real = new this.ArrayType(other.real) + this.imag = new this.ArrayType(other.imag) + } else { + this.ArrayType = opt_array_type || DefaultArrayType + // other can be either an array or a number. + this.real = new this.ArrayType(other) + this.imag = new this.ArrayType(this.real.length) + } + + this.length = this.real.length + } + + ComplexArray.prototype.toString = function() { + var components = [] + + this.forEach(function(c_value, i) { + components.push( + '(' + + c_value.real.toFixed(2) + ',' + + c_value.imag.toFixed(2) + + ')' + ) + }) + + return '[' + components.join(',') + ']' + } + + // In-place mapper. + ComplexArray.prototype.map = function(mapper) { + var + i, + n = this.length, + // For GC efficiency, pass a single c_value object to the mapper. + c_value = {} + + for (i = 0; i < n; i++) { + c_value.real = this.real[i] + c_value.imag = this.imag[i] + mapper(c_value, i, n) + this.real[i] = c_value.real + this.imag[i] = c_value.imag + } + + return this + } + + ComplexArray.prototype.forEach = function(iterator) { + var + i, + n = this.length, + // For consistency with .map. + c_value = {} + + for (i = 0; i < n; i++) { + c_value.real = this.real[i] + c_value.imag = this.imag[i] + iterator(c_value, i, n) + } + } + + ComplexArray.prototype.conjugate = function() { + return (new ComplexArray(this)).map(function(value) { + value.imag *= -1 + }) + } + + // Helper so we can make ArrayType objects returned have similar interfaces + // to ComplexArrays. + function iterable(obj) { + if (!obj.forEach) + obj.forEach = function(iterator) { + var i, n = this.length + + for (i = 0; i < n; i++) + iterator(this[i], i, n) + } + + return obj + } + + ComplexArray.prototype.magnitude = function() { + var mags = new this.ArrayType(this.length) + + this.forEach(function(value, i) { + mags[i] = sqrt(sqr(value.real) + sqr(value.imag)) + }) + + // ArrayType will not necessarily be iterable: make it so. + return iterable(mags) + } +}(typeof exports === 'undefined' && (this.complex_array = {}) || exports) diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/lib/fft.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/lib/fft.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,225 @@ +'use strict'; + +!function(exports, complex_array) { + + var + ComplexArray = complex_array.ComplexArray, + // Math constants and functions we need. + PI = Math.PI, + SQRT1_2 = Math.SQRT1_2, + sqrt = Math.sqrt, + cos = Math.cos, + sin = Math.sin + + ComplexArray.prototype.FFT = function() { + return FFT(this, false) + } + + exports.FFT = function(input) { + return ensureComplexArray(input).FFT() + } + + ComplexArray.prototype.InvFFT = function() { + return FFT(this, true) + } + + exports.InvFFT = function(input) { + return ensureComplexArray(input).InvFFT() + } + + // Applies a frequency-space filter to input, and returns the real-space + // filtered input. + // filterer accepts freq, i, n and modifies freq.real and freq.imag. + ComplexArray.prototype.frequencyMap = function(filterer) { + return this.FFT().map(filterer).InvFFT() + } + + exports.frequencyMap = function(input, filterer) { + return ensureComplexArray(input).frequencyMap(filterer) + } + + function ensureComplexArray(input) { + return complex_array.isComplexArray(input) && input || + new ComplexArray(input) + } + + function FFT(input, inverse) { + var n = input.length + + if (n & (n - 1)) { + return FFT_Recursive(input, inverse) + } else { + return FFT_2_Iterative(input, inverse) + } + } + + function FFT_Recursive(input, inverse) { + var + n = input.length, + // Counters. + i, j, + output, + // Complex multiplier and its delta. + f_r, f_i, del_f_r, del_f_i, + // Lowest divisor and remainder. + p, m, + normalisation, + recursive_result, + _swap, _real, _imag + + if (n === 1) { + return input + } + + output = new ComplexArray(n, input.ArrayType) + + // Use the lowest odd factor, so we are able to use FFT_2_Iterative in the + // recursive transforms optimally. + p = LowestOddFactor(n) + m = n / p + normalisation = 1 / sqrt(p) + recursive_result = new ComplexArray(m, input.ArrayType) + + // Loops go like O(n Σ p_i), where p_i are the prime factors of n. + // for a power of a prime, p, this reduces to O(n p log_p n) + for(j = 0; j < p; j++) { + for(i = 0; i < m; i++) { + recursive_result.real[i] = input.real[i * p + j] + recursive_result.imag[i] = input.imag[i * p + j] + } + // Don't go deeper unless necessary to save allocs. + if (m > 1) { + recursive_result = FFT(recursive_result, inverse) + } + + del_f_r = cos(2*PI*j/n) + del_f_i = (inverse ? -1 : 1) * sin(2*PI*j/n) + f_r = 1 + f_i = 0 + + for(i = 0; i < n; i++) { + _real = recursive_result.real[i % m] + _imag = recursive_result.imag[i % m] + + output.real[i] += f_r * _real - f_i * _imag + output.imag[i] += f_r * _imag + f_i * _real + + _swap = f_r * del_f_r - f_i * del_f_i + f_i = f_r * del_f_i + f_i * del_f_r + f_r = _swap + } + } + + // Copy back to input to match FFT_2_Iterative in-placeness + // TODO: faster way of making this in-place? + for(i = 0; i < n; i++) { + input.real[i] = normalisation * output.real[i] + input.imag[i] = normalisation * output.imag[i] + } + + return input + } + + function FFT_2_Iterative(input, inverse) { + var + n = input.length, + // Counters. + i, j, + output, output_r, output_i, + // Complex multiplier and its delta. + f_r, f_i, del_f_r, del_f_i, temp, + // Temporary loop variables. + l_index, r_index, + left_r, left_i, right_r, right_i, + // width of each sub-array for which we're iteratively calculating FFT. + width + + output = BitReverseComplexArray(input) + output_r = output.real + output_i = output.imag + // Loops go like O(n log n): + // width ~ log n; i,j ~ n + width = 1 + while (width < n) { + del_f_r = cos(PI/width) + del_f_i = (inverse ? -1 : 1) * sin(PI/width) + for (i = 0; i < n/(2*width); i++) { + f_r = 1 + f_i = 0 + for (j = 0; j < width; j++) { + l_index = 2*i*width + j + r_index = l_index + width + + left_r = output_r[l_index] + left_i = output_i[l_index] + right_r = f_r * output_r[r_index] - f_i * output_i[r_index] + right_i = f_i * output_r[r_index] + f_r * output_i[r_index] + + output_r[l_index] = SQRT1_2 * (left_r + right_r) + output_i[l_index] = SQRT1_2 * (left_i + right_i) + output_r[r_index] = SQRT1_2 * (left_r - right_r) + output_i[r_index] = SQRT1_2 * (left_i - right_i) + temp = f_r * del_f_r - f_i * del_f_i + f_i = f_r * del_f_i + f_i * del_f_r + f_r = temp + } + } + width <<= 1 + } + + return output + } + + function BitReverseIndex(index, n) { + var bitreversed_index = 0 + + while (n > 1) { + bitreversed_index <<= 1 + bitreversed_index += index & 1 + index >>= 1 + n >>= 1 + } + return bitreversed_index + } + + function BitReverseComplexArray(array) { + var n = array.length, + flips = {}, + swap, + i + + for(i = 0; i < n; i++) { + var r_i = BitReverseIndex(i, n) + + if (flips.hasOwnProperty(i) || flips.hasOwnProperty(r_i)) continue + + swap = array.real[r_i] + array.real[r_i] = array.real[i] + array.real[i] = swap + + swap = array.imag[r_i] + array.imag[r_i] = array.imag[i] + array.imag[i] = swap + + flips[i] = flips[r_i] = true + } + + return array + } + + function LowestOddFactor(n) { + var factor = 3, + sqrt_n = sqrt(n) + + while(factor <= sqrt_n) { + if (n % factor === 0) return factor + factor = factor + 2 + } + return n + } + +}( + typeof exports === 'undefined' && (this.fft = {}) || exports, + typeof require === 'undefined' && (this.complex_array) || + require('./complex_array') +) diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/lib/fft_image.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/lib/fft_image.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,87 @@ +'use strict'; + +!function(complex_array, fft) { + + var ComplexArray = complex_array.ComplexArray + + fft.FFTImageDataRGBA = function(data, nx, ny) { + var rgb = splitRGB(data) + + return mergeRGB( + FFT2D(new ComplexArray(rgb[0], Float32Array), nx, ny), + FFT2D(new ComplexArray(rgb[1], Float32Array), nx, ny), + FFT2D(new ComplexArray(rgb[2], Float32Array), nx, ny) + ) + } + + function splitRGB(data) { + var n = data.length / 4, + r = new Uint8ClampedArray(n), + g = new Uint8ClampedArray(n), + b = new Uint8ClampedArray(n), + i + + for(i = 0; i < n; i++) { + r[i] = data[4 * i ] + g[i] = data[4 * i + 1] + b[i] = data[4 * i + 2] + } + + return [r, g, b] + } + + function mergeRGB(r, g, b) { + var n = r.length, + output = new ComplexArray(n * 4), + i + + for(i = 0; i < n; i++) { + output.real[4 * i ] = r.real[i] + output.imag[4 * i ] = r.imag[i] + output.real[4 * i + 1] = g.real[i] + output.imag[4 * i + 1] = g.imag[i] + output.real[4 * i + 2] = b.real[i] + output.imag[4 * i + 2] = b.imag[i] + } + + return output + } + + function FFT2D(input, nx, ny, inverse) { + var i, j, + transform = inverse ? 'InvFFT' : 'FFT', + output = new ComplexArray(input.length, input.ArrayType), + row = new ComplexArray(nx, input.ArrayType), + col = new ComplexArray(ny, input.ArrayType) + + for(j = 0; j < ny; j++) { + row.map(function(v, i) { + v.real = input.real[i + j * nx] + v.imag = input.imag[i + j * nx] + }) + row[transform]().forEach(function(v, i) { + output.real[i + j * nx] = v.real + output.imag[i + j * nx] = v.imag + }) + } + + for(i = 0; i < nx; i++) { + col.map(function(v, j) { + v.real = output.real[i + j * nx] + v.imag = output.imag[i + j * nx] + }) + col[transform]().forEach(function(v, j) { + output.real[i + j * nx] = v.real + output.imag[i + j * nx] = v.imag + }) + } + + return output + } + +}( + typeof require === 'undefined' && (this.complex_array) || + require('./complex_array'), + typeof require === 'undefined' && (this.fft) || + require('./fft') +) diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/package.json --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/package.json Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,23 @@ +{ + "name": "jsfft", + "version": "0.0.1", + "author": "Nick Jones", + "description": "A lightweight FFT implementation for Javascript", + "repository": { + "type": "git", + "url": "https://github.com/dntj/jsfft.git" + }, + "engines": { + "node": "*" + }, + + "main": "./lib/fft", + + "scripts": { + "test": "mocha" + }, + + "devDependencies": { + "mocha": ">=1.1.0" + } +} diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/test/test-complex_array.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/test/test-complex_array.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,125 @@ +var assert = require('assert'), + complex_array = require('../lib/complex_array'), + ComplexArray = complex_array.ComplexArray + +function assertArrayEquals(first, second) { + var message = first + ' != ' + second + + first.forEach(function(item, i) { + assert.equal(item, second[i], message) + }) +} + +describe('isComplexArray', function() { + it('should correctly identify complex arrays', function() { + assert.ok(!complex_array.isComplexArray(1)) + assert.ok(!complex_array.isComplexArray([1,2,3])) + assert.ok(!complex_array.isComplexArray('string')) + + assert.ok(complex_array.isComplexArray(new ComplexArray(1))) + }) +}) + +describe('ComplexArray', function() { + describe('#__constructor__()', function() { + it('should construct from a number', function() { + var a = new ComplexArray(10) + assert.equal(10, a.real.length) + assert.equal(10, a.imag.length) + assert.equal(0, a.real[0]) + assert.equal(0, a.imag[0]) + }) + + it('should construct from a number with a type', function() { + var a = new ComplexArray(10, Int32Array) + assert.equal(Int32Array, a.ArrayType) + assert.equal(10, a.real.length) + assert.equal(10, a.imag.length) + assert.equal(0, a.real[0]) + assert.equal(0, a.imag[0]) + }) + + it('should contruct from a real array', function() { + var a = new ComplexArray([1, 2]) + assertArrayEquals([1, 2], a.real) + assertArrayEquals([0, 0], a.imag) + }) + + it('should contruct from a real array with a type', function() { + var a = new ComplexArray([1, 2], Int32Array) + assert.equal(Int32Array, a.ArrayType) + assertArrayEquals([1, 2], a.real) + assertArrayEquals([0, 0], a.imag) + }) + + it('should contruct from another complex array', function() { + var a = new ComplexArray(new ComplexArray([1, 2])) + assertArrayEquals([1, 2], a.real) + assertArrayEquals([0, 0], a.imag) + }) + }) + + describe('#map()', function() { + it('should alter all values', function() { + var a = new ComplexArray([1, 2]) + + a.map(function(value, i) { + value.real *= 10 + value.imag = i + }) + assertArrayEquals([10, 20], a.real) + assertArrayEquals([0, 1], a.imag) + }) + }) + + describe('#forEach()', function() { + it('should touch every value', function() { + var + a = new ComplexArray([1, 2]), + sum = 0 + + a.imag[0] = 4 + a.imag[1] = 8 + a.forEach(function(value, i) { + sum += value.real + sum += value.imag + }) + assert.equal(15, sum) + }) + }) + + describe('#conjugate()', function() { + it('should multiply a number', function() { + var + a = new ComplexArray([1, 2]), + b + + a.imag[0] = 1 + a.imag[1] = -2 + b = a.conjugate() + assertArrayEquals([1, 2], b.real) + assertArrayEquals([-1, 2], b.imag) + }) + }) + + describe('#magnitude()', function() { + it('should give the an array of magnitudes', function() { + var a = new ComplexArray([1, 3]) + + a.imag[0] = 0 + a.imag[1] = 4 + assertArrayEquals([1, 5], a.magnitude()) + }) + + it('should return an iterable ArrayType object', function() { + var + sum = 0, + a = new ComplexArray([1, 2]) + + a.magnitude().forEach(function(value, i) { + sum += value + }) + assert.equal(3, sum) + }) + }) +}) \ No newline at end of file diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/test/test-fft.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/test/test-fft.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,136 @@ +require('./test_helper') + +var ComplexArray = require('../lib/complex_array').ComplexArray + +describe('fft', function() { + describe('#FFT()', function() { + describe('on N=4 Arrays', function() { + it('should return a single frequency given a constant array', function() { + assertFFTMatches([1, 1, 1, 1], new ComplexArray([2, 0, 0, 0])) + }) + + it('should return flat with a delta function input', function() { + assertFFTMatches([1, 0, 0, 0], new ComplexArray([0.5, 0.5, 0.5, 0.5])) + }) + + it('should return a single high freq', function() { + assertFFTMatches([1, -1, 1, -1], new ComplexArray([0, 0, 2, 0])) + }) + + it('should return a single low freq', function() { + assertFFTMatches([1, 0, -1, -0], new ComplexArray([0, 1, 0, 1])) + }) + + it('should return a high freq and DC', function() { + assertFFTMatches([1, 0, 1, 0], new ComplexArray([1, 0, 1, 0])) + }) + }) + + describe('on N=6 Arrays', function() { + it('should return a single frequency given a constant array', function() { + assertFFTMatches( + [1, 1, 1, 1, 1, 1], + new ComplexArray([Math.sqrt(6), 0, 0, 0, 0, 0]) + ) + }) + + it('should return flat with a delta function input', function() { + var a = 1 / Math.sqrt(6) + assertFFTMatches( + [1, 0, 0, 0, 0, 0], + new ComplexArray([a, a, a, a, a, a]) + ) + }) + }) + + describe('on N=`prime` Arrays', function() { + it('should match the DFT', function() { + var a = new ComplexArray(13) + + a.map(function(value) { + value.real = Math.random() + value.imag = Math.random() + }) + + assertFFTMatchesDFT(a) + }) + }) + + describe('on N=512 Arrays', function() { + it('should match the DFT', function() { + var a = new ComplexArray(512) + + a.map(function(value) { + value.real = Math.random() + value.imag = Math.random() + }) + + assertFFTMatchesDFT(a) + }) + }) + + describe('on N=900 Arrays', function() { + it('should match the DFT', function() { + var a = new ComplexArray(900) + + a.map(function(value) { + value.real = Math.random() + value.imag = Math.random() + }) + + assertFFTMatchesDFT(a) + }) + }) + }) + + describe('#frequencyMap()', function() { + it('should not modify the original', function() { + var + original = new ComplexArray([1, 2, 3, 4]), + filtered = original.frequencyMap(function() {}) + + assertComplexArraysAlmostEqual(original, filtered) + }) + + it('should halve the original', function() { + var + original = new ComplexArray([1, 2, 3, 4]), + filtered = original.frequencyMap(function(value, i) { + value.real /= 2 + value.imag /= 2 + }) + + assertComplexArraysAlmostEqual( + new ComplexArray([0.5, 1, 1.5, 2]), filtered) + }) + + it('should return zeroed ComplexArray', function() { + var + original = new ComplexArray([1, 2, 3, 4]), + filtered = original.frequencyMap(function(value, i) { + value.real = value.imag = 0 + }) + + assertComplexArraysAlmostEqual(new ComplexArray([0, 0, 0, 0]), filtered) + }) + + it('should shift the original', function() { + var + original = new ComplexArray([1, 2, 3, 4]), + filtered = original.frequencyMap(function(value, i) { + var + // Multiply by a phase to shift the original. + real_multiplier = i % 2 ? 0 : (1 - i), + imag_multiplier = i % 2 ? (2 - i) : 0, + swap_real = value.real, + swap_imag = value.imag + + value.real = real_multiplier * swap_real - imag_multiplier * swap_imag + value.imag = real_multiplier * swap_imag + imag_multiplier * swap_real + }) + + assertComplexArraysAlmostEqual(new ComplexArray([4, 1, 2, 3]), filtered) + }) + }) +}) + diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/test/test-fft_image.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/test/test-fft_image.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,70 @@ +require('../lib/fft_image') +require('./test_helper') + +var assert = require('assert'), + fft_lib = require('../lib/fft'), + complex_array_lib = require('../lib/complex_array'), + ComplexArray = complex_array_lib.ComplexArray + +function randomImageData(n) { + var array = new Uint8ClampedArray(n), + i + + for(i = 0; i < n; i++) { + array[i] = Math.random() * 128 + } + return array +} + +function mergeRGBA(r, g, b, a) { + var n = r.length, + output = new Uint8ClampedArray(4 * n), + i + + for(i = 0; i < n; i++) { + output[4 * i ] = r[i] + output[4 * i + 1] = g[i] + output[4 * i + 2] = b[i] + output[4 * i + 3] = a[i] + } + return output +} + +function mergeComplexRGBA(r, g, b, a) { + var n = r.length, + output = new ComplexArray(4 * n), + i + + for(i = 0; i < n; i++) { + output.real[4 * i ] = r.real[i] + output.imag[4 * i ] = r.imag[i] + output.real[4 * i + 1] = g.real[i] + output.imag[4 * i + 1] = g.imag[i] + output.real[4 * i + 2] = b.real[i] + output.imag[4 * i + 2] = b.imag[i] + output.real[4 * i + 3] = a.real[i] + output.imag[4 * i + 3] = a.imag[i] + } + + return output +} + +describe('fft', function() { + describe('#FFTImageDataRGBA()', function() { + it('transforms independent channels', function() { + var n = 48, + r = randomImageData(n), + g = randomImageData(n), + b = randomImageData(n), + a = new Uint8ClampedArray(n), + data = mergeRGBA(r, g, b, a), + expected = mergeComplexRGBA(DFT(r), DFT(g), DFT(b), DFT(a)), + output = fft_lib.FFTImageDataRGBA(data, n, 1) + + assertComplexArraysAlmostEqual(expected, output) + }) + + xit('transforms in 2D') + }) +}) + diff -r e705de983b67 -r 66f9fd5ac611 fft/jsfft/test/test_helper.js --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/jsfft/test/test_helper.js Wed Oct 07 13:46:38 2015 +0100 @@ -0,0 +1,83 @@ +!function() { + + var EPSILON = 1e-4, + assert = require('assert'), + complex_array_lib = require('../lib/complex_array'), + fft_lib = require('../lib/fft'), + ComplexArray = complex_array_lib.ComplexArray, + isComplexArray = complex_array_lib.isComplexArray, + PI = Math.PI, + SQRT2 = Math.SQRT2, + SQRT1_2 = Math.SQRT1_2, + cos = Math.cos, + sin = Math.sin, + sqrt = Math.sqrt + + global.assertComplexArraysAlmostEqual = function(first, second) { + var message = second + ' != ' + first + + assert.equal(first.length, second.length, message) + + first.forEach(function(value, i) { + assertApproximatelyEqual(value.real, second.real[i], message) + assertApproximatelyEqual(value.imag, second.imag[i], message) + }) + } + + global.assertFFTMatches = function(original, expected) { + var transformed, copy + + if (!isComplexArray(expected)) { + throw TypeError('expected match should be a ComplexArray') + } + + copy = new ComplexArray(original) + transformed = fft_lib.FFT(original) + assertComplexArraysAlmostEqual(expected, transformed) + assertComplexArraysAlmostEqual(copy, fft_lib.InvFFT(transformed)) + } + + global.assertFFTMatchesDFT = function(input) { + input = new ComplexArray(input) + + assertComplexArraysAlmostEqual(DFT(input), fft_lib.FFT(input)) + } + + global.DFT = function(input) { + var n = input.length, + amplitude = 1 / sqrt(n), + output = new ComplexArray(input), + phase = {real: 0, imag: 0}, + delta = {real: 0, imag: 0}, + i, j, + _swap + + if (!isComplexArray(input)) { + input = new ComplexArray(input) + } + + for(i = 0; i < n; i++) { + output.real[i] = 0, output.imag[i] = 0 + phase.real = 1, phase.imag = 0 + delta.real = cos(2*PI*i/n), delta.imag = sin(2*PI*i/n) + + for(j = 0; j < n; j++) { + output.real[i] += phase.real * input.real[j] - phase.imag * input.imag[j] + output.imag[i] += phase.real * input.imag[j] + phase.imag * input.real[j] + _swap = phase.real + phase.real = phase.real * delta.real - phase.imag * delta.imag + phase.imag = _swap * delta.imag + phase.imag * delta.real + } + output.real[i] *= amplitude + output.imag[i] *= amplitude + } + + return output + } + + function assertApproximatelyEqual(first, second, message) { + var delta = Math.abs(first - second) + assert.ok(delta < EPSILON, message) + } + +}() diff -r e705de983b67 -r 66f9fd5ac611 fft/test.html --- a/fft/test.html Wed Oct 07 08:57:11 2015 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - -

Results

- -

- - - - - - - - - - - - - - - - - - - - -
ImplementationResultTime (first half)Time (second half)Rate (second half)
Nayuki
Nayuki (obj)
Nockert
Dntj
Cross
KissFFT
FFTW
- -

Notes

- -
    -
  • Nayuki: in-place single-precision complex-complex. Around 7kb.
    • -
  • Nayuki (obj): Nayuki with the sin/cos tables pre-calculated on object construction. Around 4kb.
    • -
  • Nockert: double-precision real-complex. Around 25kb.
    • -
  • Dntj: double-precision complex-complex. Forward - transform is scaled and I've scaled it back again here, which may - introduce rounding error. Around 10kb.
    • -
  • Cross: double-precision real-complex in C, compiled - with Emscripten. This is considered a slow implementation amongst - native code ones. Around 60kb.
    • -
  • KissFFT: single-precision real-complex in C, compiled - with Emscripten. A reasonably sophisticated implementation. Around - 70kb.
    • -
  • FFTW: single-precision real-complex in C, compiled with - Emscripten. GPL licensed. Around 3Mb.
    • -
- - -