Chris@25: /* Copyright (c) 2012, Jens Nockert , Jussi Kalliokoski Chris@25: * All rights reserved. Chris@25: * Chris@25: * Redistribution and use in source and binary forms, with or without Chris@25: * modification, are permitted provided that the following conditions are met: Chris@25: * Chris@25: * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Chris@25: * 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. Chris@25: * Chris@25: * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND Chris@25: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED Chris@25: * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE Chris@25: * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR Chris@25: * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES Chris@25: * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; Chris@25: * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND Chris@25: * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT Chris@25: * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS Chris@25: * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Chris@25: */ Chris@25: Chris@25: if (!FFT) { Chris@25: var FFT = {} Chris@25: } Chris@25: Chris@25: void function (namespace) { Chris@25: "use strict" Chris@25: Chris@25: function butterfly2(output, outputOffset, outputStride, fStride, state, m) { Chris@25: var t = state.twiddle Chris@25: Chris@25: for (var i = 0; i < m; i++) { Chris@25: var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] Chris@25: var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] Chris@25: Chris@25: var t1_r = t[2 * ((0) + (fStride) * (i))], t1_i = t[2 * ((0) + (fStride) * (i)) + 1] Chris@25: Chris@25: var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r Chris@25: Chris@25: var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i Chris@25: var r1_r = s0_r - v1_r, r1_i = s0_i - v1_i Chris@25: Chris@25: output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i Chris@25: output[2 * ((outputOffset) + (outputStride) * (i + m))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] = r1_i Chris@25: } Chris@25: } Chris@25: Chris@25: function butterfly3(output, outputOffset, outputStride, fStride, state, m) { Chris@25: var t = state.twiddle Chris@25: var m1 = m, m2 = 2 * m Chris@25: var fStride1 = fStride, fStride2 = 2 * fStride Chris@25: Chris@25: var e = t[2 * ((0) + (fStride) * (m)) + 1] Chris@25: Chris@25: for (var i = 0; i < m; i++) { Chris@25: var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] Chris@25: Chris@25: var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] Chris@25: var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1] Chris@25: var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r Chris@25: Chris@25: var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] Chris@25: var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1] Chris@25: var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r Chris@25: Chris@25: var i0_r = v1_r + v2_r, i0_i = v1_i + v2_i Chris@25: Chris@25: var r0_r = s0_r + i0_r, r0_i = s0_i + i0_i Chris@25: output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i Chris@25: Chris@25: var i1_r = s0_r - i0_r * 0.5 Chris@25: var i1_i = s0_i - i0_i * 0.5 Chris@25: Chris@25: var i2_r = (v1_r - v2_r) * e Chris@25: var i2_i = (v1_i - v2_i) * e Chris@25: Chris@25: var r1_r = i1_r - i2_i Chris@25: var r1_i = i1_i + i2_r Chris@25: output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i Chris@25: Chris@25: var r2_r = i1_r + i2_i Chris@25: var r2_i = i1_i - i2_r Chris@25: output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i Chris@25: } Chris@25: } Chris@25: Chris@25: function butterfly4(output, outputOffset, outputStride, fStride, state, m) { Chris@25: var t = state.twiddle Chris@25: var m1 = m, m2 = 2 * m, m3 = 3 * m Chris@25: var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride Chris@25: Chris@25: for (var i = 0; i < m; i++) { Chris@25: var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1] Chris@25: Chris@25: var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] Chris@25: var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1] Chris@25: var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r Chris@25: Chris@25: var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] Chris@25: var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1] Chris@25: var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r Chris@25: Chris@25: var s3_r = output[2 * ((outputOffset) + (outputStride) * (i + m3))], s3_i = output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] Chris@25: var t3_r = t[2 * ((0) + (fStride3) * (i))], t3_i = t[2 * ((0) + (fStride3) * (i)) + 1] Chris@25: var v3_r = s3_r * t3_r - s3_i * t3_i, v3_i = s3_r * t3_i + s3_i * t3_r Chris@25: Chris@25: var i0_r = s0_r + v2_r, i0_i = s0_i + v2_i Chris@25: var i1_r = s0_r - v2_r, i1_i = s0_i - v2_i Chris@25: var i2_r = v1_r + v3_r, i2_i = v1_i + v3_i Chris@25: var i3_r = v1_r - v3_r, i3_i = v1_i - v3_i Chris@25: Chris@25: var r0_r = i0_r + i2_r, r0_i = i0_i + i2_i Chris@25: Chris@25: if (state.inverse) { Chris@25: var r1_r = i1_r - i3_i Chris@25: var r1_i = i1_i + i3_r Chris@25: } else { Chris@25: var r1_r = i1_r + i3_i Chris@25: var r1_i = i1_i - i3_r Chris@25: } Chris@25: Chris@25: var r2_r = i0_r - i2_r, r2_i = i0_i - i2_i Chris@25: Chris@25: if (state.inverse) { Chris@25: var r3_r = i1_r + i3_i Chris@25: var r3_i = i1_i - i3_r Chris@25: } else { Chris@25: var r3_r = i1_r - i3_i Chris@25: var r3_i = i1_i + i3_r Chris@25: } Chris@25: Chris@25: output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i Chris@25: output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i Chris@25: output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i Chris@25: output[2 * ((outputOffset) + (outputStride) * (i + m3))] = r3_r, output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] = r3_i Chris@25: } Chris@25: } Chris@25: Chris@25: function butterfly(output, outputOffset, outputStride, fStride, state, m, p) { Chris@25: var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p) Chris@25: Chris@25: for (var u = 0; u < m; u++) { Chris@25: for (var q1 = 0, k = u; q1 < p; q1++, k += m) { Chris@25: var x0_r = output[2 * ((outputOffset) + (outputStride) * (k))], x0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1] Chris@25: scratch[2 * (q1)] = x0_r, scratch[2 * (q1) + 1] = x0_i Chris@25: } Chris@25: Chris@25: for (var q1 = 0, k = u; q1 < p; q1++, k += m) { Chris@25: var tOffset = 0 Chris@25: Chris@25: var x0_r = scratch[2 * (0)], x0_i = scratch[2 * (0) + 1] Chris@25: output[2 * ((outputOffset) + (outputStride) * (k))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = x0_i Chris@25: Chris@25: for (var q = 1; q < p; q++) { Chris@25: tOffset = (tOffset + fStride * k) % n Chris@25: Chris@25: var s0_r = output[2 * ((outputOffset) + (outputStride) * (k))], s0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1] Chris@25: Chris@25: var s1_r = scratch[2 * (q)], s1_i = scratch[2 * (q) + 1] Chris@25: var t1_r = t[2 * (tOffset)], t1_i = t[2 * (tOffset) + 1] Chris@25: var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r Chris@25: Chris@25: var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i Chris@25: output[2 * ((outputOffset) + (outputStride) * (k))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = r0_i Chris@25: } Chris@25: } Chris@25: } Chris@25: } Chris@25: Chris@25: function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) { Chris@25: var p = factors.shift() Chris@25: var m = factors.shift() Chris@25: Chris@25: if (m == 1) { Chris@25: for (var i = 0; i < p * m; i++) { Chris@25: var x0_r = f[2 * ((fOffset) + (fStride * inputStride) * (i))], x0_i = f[2 * ((fOffset) + (fStride * inputStride) * (i)) + 1] Chris@25: output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i Chris@25: } Chris@25: } else { Chris@25: for (var i = 0; i < p; i++) { Chris@25: work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state) Chris@25: } Chris@25: } Chris@25: Chris@25: switch (p) { Chris@25: case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break Chris@25: case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break Chris@25: case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break Chris@25: default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break Chris@25: } Chris@25: } Chris@25: Chris@25: var complex = function (n, inverse) { Chris@25: if (arguments.length < 2) { Chris@25: throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'") Chris@25: } Chris@25: Chris@25: var n = ~~n, inverse = !!inverse Chris@25: Chris@25: if (n < 1) { Chris@25: throw new RangeError("n is outside range, should be positive integer, was `" + n + "'") Chris@25: } Chris@25: Chris@25: var state = { Chris@25: n: n, Chris@25: inverse: inverse, Chris@25: Chris@25: factors: [], Chris@25: twiddle: new Float64Array(2 * n), Chris@25: scratch: new Float64Array(2 * n) Chris@25: } Chris@25: Chris@25: var t = state.twiddle, theta = 2 * Math.PI / n Chris@25: Chris@25: for (var i = 0; i < n; i++) { Chris@25: if (inverse) { Chris@25: var phase = theta * i Chris@25: } else { Chris@25: var phase = -theta * i Chris@25: } Chris@25: Chris@25: t[2 * (i)] = Math.cos(phase) Chris@25: t[2 * (i) + 1] = Math.sin(phase) Chris@25: } Chris@25: Chris@25: var p = 4, v = Math.floor(Math.sqrt(n)) Chris@25: Chris@25: while (n > 1) { Chris@25: while (n % p) { Chris@25: switch (p) { Chris@25: case 4: p = 2; break Chris@25: case 2: p = 3; break Chris@25: default: p += 2; break Chris@25: } Chris@25: Chris@25: if (p > v) { Chris@25: p = n Chris@25: } Chris@25: } Chris@25: Chris@25: n /= p Chris@25: Chris@25: state.factors.push(p) Chris@25: state.factors.push(n) Chris@25: } Chris@25: Chris@25: this.state = state Chris@25: } Chris@25: Chris@25: complex.prototype.simple = function (output, input, t) { Chris@25: this.process(output, 0, 1, input, 0, 1, t) Chris@25: } Chris@25: Chris@25: complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) { Chris@25: var outputStride = ~~outputStride, inputStride = ~~inputStride Chris@25: Chris@25: var type = t == 'real' ? t : 'complex' Chris@25: Chris@25: if (outputStride < 1) { Chris@25: throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'") Chris@25: } Chris@25: Chris@25: if (inputStride < 1) { Chris@25: throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'") Chris@25: } Chris@25: Chris@25: if (type == 'real') { Chris@25: for (var i = 0; i < this.state.n; i++) { Chris@25: var x0_r = input[inputOffset + inputStride * i] Chris@25: var x0_i = 0.0 Chris@25: Chris@25: this.state.scratch[2 * (i)] = x0_r, this.state.scratch[2 * (i) + 1] = x0_i Chris@25: } Chris@25: Chris@25: work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state) Chris@25: } else { Chris@25: if (input == output) { Chris@25: work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) Chris@25: Chris@25: for (var i = 0; i < this.state.n; i++) { Chris@25: var x0_r = this.state.scratch[2 * (i)], x0_i = this.state.scratch[2 * (i) + 1] Chris@25: Chris@25: output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i Chris@25: } Chris@25: } else { Chris@25: work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state) Chris@25: } Chris@25: } Chris@25: } Chris@25: Chris@25: namespace.complex = complex Chris@25: }(FFT) Chris@25: Chris@25: module.exports = FFT