fft.js/lib/node.js

   1 /* Copyright (c) 2012, Jens Nockert <jens@ofmlabs.org>, Jussi Kalliokoski <jussi@ofmlabs.org>
   2  * All rights reserved.
   3  *
   4  * Redistribution and use in source and binary forms, with or without
   5  * modification, are permitted provided that the following conditions are met:
   6  *
   7  * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
   8  * 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.
   9  *
  10  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  11  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  12  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  13  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
  14  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  15  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  16  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  17  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  18  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  19  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  20  */
  21
  22 if (!FFT) {
  23         var FFT = {}
  24 }
  25
  26 void function (namespace) {
  27         "use strict"
  28
  29         function butterfly2(output, outputOffset, outputStride, fStride, state, m) {
  30                 var t = state.twiddle
  31
  32                 for (var i = 0; i < m; i++) {
  33                         var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1]
  34                         var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1]
  35
  36                         var t1_r = t[2 * ((0) + (fStride) * (i))], t1_i = t[2 * ((0) + (fStride) * (i)) + 1]
  37
  38                         var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r
  39
  40                         var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i
  41                         var r1_r = s0_r - v1_r, r1_i = s0_i - v1_i
  42
  43                         output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i
  44                         output[2 * ((outputOffset) + (outputStride) * (i + m))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m)) + 1] = r1_i
  45                 }
  46         }
  47
  48         function butterfly3(output, outputOffset, outputStride, fStride, state, m) {
  49                 var t = state.twiddle
  50                 var m1 = m, m2 = 2 * m
  51                 var fStride1 = fStride, fStride2 = 2 * fStride
  52
  53                 var e = t[2 * ((0) + (fStride) * (m)) + 1]
  54
  55                 for (var i = 0; i < m; i++) {
  56                         var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1]
  57
  58                         var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1]
  59                         var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1]
  60                         var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r
  61
  62                         var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1]
  63                         var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1]
  64                         var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r
  65
  66                         var i0_r = v1_r + v2_r, i0_i = v1_i + v2_i
  67
  68                         var r0_r = s0_r + i0_r, r0_i = s0_i + i0_i
  69                         output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i
  70
  71                         var i1_r = s0_r - i0_r * 0.5
  72                         var i1_i = s0_i - i0_i * 0.5
  73
  74                         var i2_r = (v1_r - v2_r) * e
  75                         var i2_i = (v1_i - v2_i) * e
  76
  77                         var r1_r = i1_r - i2_i
  78                         var r1_i = i1_i + i2_r
  79                         output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i
  80
  81                         var r2_r = i1_r + i2_i
  82                         var r2_i = i1_i - i2_r
  83                         output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i
  84                 }
  85         }
  86
  87         function butterfly4(output, outputOffset, outputStride, fStride, state, m) {
  88                 var t = state.twiddle
  89                 var m1 = m, m2 = 2 * m, m3 = 3 * m
  90                 var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride
  91
  92                 for (var i = 0; i < m; i++) {
  93                         var s0_r = output[2 * ((outputOffset) + (outputStride) * (i))], s0_i = output[2 * ((outputOffset) + (outputStride) * (i)) + 1]
  94
  95                         var s1_r = output[2 * ((outputOffset) + (outputStride) * (i + m1))], s1_i = output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1]
  96                         var t1_r = t[2 * ((0) + (fStride1) * (i))], t1_i = t[2 * ((0) + (fStride1) * (i)) + 1]
  97                         var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r
  98
  99                         var s2_r = output[2 * ((outputOffset) + (outputStride) * (i + m2))], s2_i = output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1]
 100                         var t2_r = t[2 * ((0) + (fStride2) * (i))], t2_i = t[2 * ((0) + (fStride2) * (i)) + 1]
 101                         var v2_r = s2_r * t2_r - s2_i * t2_i, v2_i = s2_r * t2_i + s2_i * t2_r
 102
 103                         var s3_r = output[2 * ((outputOffset) + (outputStride) * (i + m3))], s3_i = output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1]
 104                         var t3_r = t[2 * ((0) + (fStride3) * (i))], t3_i = t[2 * ((0) + (fStride3) * (i)) + 1]
 105                         var v3_r = s3_r * t3_r - s3_i * t3_i, v3_i = s3_r * t3_i + s3_i * t3_r
 106
 107                         var i0_r = s0_r + v2_r, i0_i = s0_i + v2_i
 108                         var i1_r = s0_r - v2_r, i1_i = s0_i - v2_i
 109                         var i2_r = v1_r + v3_r, i2_i = v1_i + v3_i
 110                         var i3_r = v1_r - v3_r, i3_i = v1_i - v3_i
 111
 112                         var r0_r = i0_r + i2_r, r0_i = i0_i + i2_i
 113
 114                         if (state.inverse) {
 115                                 var r1_r = i1_r - i3_i
 116                                 var r1_i = i1_i + i3_r
 117                         } else {
 118                                 var r1_r = i1_r + i3_i
 119                                 var r1_i = i1_i - i3_r
 120                         }
 121
 122                         var r2_r = i0_r - i2_r, r2_i = i0_i - i2_i
 123
 124                         if (state.inverse) {
 125                                 var r3_r = i1_r + i3_i
 126                                 var r3_i = i1_i - i3_r
 127                         } else {
 128                                 var r3_r = i1_r - i3_i
 129                                 var r3_i = i1_i + i3_r
 130                         }
 131
 132                         output[2 * ((outputOffset) + (outputStride) * (i))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = r0_i
 133                         output[2 * ((outputOffset) + (outputStride) * (i + m1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (i + m1)) + 1] = r1_i
 134                         output[2 * ((outputOffset) + (outputStride) * (i + m2))] = r2_r, output[2 * ((outputOffset) + (outputStride) * (i + m2)) + 1] = r2_i
 135                         output[2 * ((outputOffset) + (outputStride) * (i + m3))] = r3_r, output[2 * ((outputOffset) + (outputStride) * (i + m3)) + 1] = r3_i
 136                 }
 137         }
 138
 139         function butterfly(output, outputOffset, outputStride, fStride, state, m, p) {
 140                 var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p)
 141
 142                 for (var u = 0; u < m; u++) {
 143                         for (var q1 = 0, k = u; q1 < p; q1++, k += m) {
 144                                 var x0_r = output[2 * ((outputOffset) + (outputStride) * (k))], x0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1]
 145                                 scratch[2 * (q1)] = x0_r, scratch[2 * (q1) + 1] = x0_i
 146                         }
 147
 148                         for (var q1 = 0, k = u; q1 < p; q1++, k += m) {
 149                                 var tOffset = 0
 150
 151                                 var x0_r = scratch[2 * (0)], x0_i = scratch[2 * (0) + 1]
 152                                 output[2 * ((outputOffset) + (outputStride) * (k))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = x0_i
 153
 154                                 for (var q = 1; q < p; q++) {
 155                                         tOffset = (tOffset + fStride * k) % n
 156
 157                                         var s0_r = output[2 * ((outputOffset) + (outputStride) * (k))], s0_i = output[2 * ((outputOffset) + (outputStride) * (k)) + 1]
 158
 159                                         var s1_r = scratch[2 * (q)], s1_i = scratch[2 * (q) + 1]
 160                                         var t1_r = t[2 * (tOffset)], t1_i = t[2 * (tOffset) + 1]
 161                                         var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r
 162
 163                                         var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i
 164                                         output[2 * ((outputOffset) + (outputStride) * (k))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (k)) + 1] = r0_i
 165                                 }
 166                         }
 167                 }
 168         }
 169
 170         function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) {
 171                 var p = factors.shift()
 172                 var m = factors.shift()
 173
 174                 if (m == 1) {
 175                         for (var i = 0; i < p * m; i++) {
 176                                 var x0_r = f[2 * ((fOffset) + (fStride * inputStride) * (i))], x0_i = f[2 * ((fOffset) + (fStride * inputStride) * (i)) + 1]
 177                                 output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i
 178                         }
 179                 } else {
 180                         for (var i = 0; i < p; i++) {
 181                                 work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state)
 182                         }
 183                 }
 184
 185                 switch (p) {
 186                         case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break
 187                         case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break
 188                         case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break
 189                         default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break
 190                 }
 191         }
 192
 193         var complex = function (n, inverse) {
 194                 if (arguments.length < 2) {
 195                         throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'")
 196                 }
 197
 198                 var n = ~~n, inverse = !!inverse
 199
 200                 if (n < 1) {
 201                         throw new RangeError("n is outside range, should be positive integer, was `" + n + "'")
 202                 }
 203
 204                 var state = {
 205                         n: n,
 206                         inverse: inverse,
 207
 208                         factors: [],
 209                         twiddle: new Float64Array(2 * n),
 210                         scratch: new Float64Array(2 * n)
 211                 }
 212
 213                 var t = state.twiddle, theta = 2 * Math.PI / n
 214
 215                 for (var i = 0; i < n; i++) {
 216                         if (inverse) {
 217                                 var phase =  theta * i
 218                         } else {
 219                                 var phase = -theta * i
 220                         }
 221
 222                         t[2 * (i)] = Math.cos(phase)
 223                         t[2 * (i) + 1] = Math.sin(phase)
 224                 }
 225
 226                 var p = 4, v = Math.floor(Math.sqrt(n))
 227
 228                 while (n > 1) {
 229                         while (n % p) {
 230                                 switch (p) {
 231                                         case 4: p = 2; break
 232                                         case 2: p = 3; break
 233                                         default: p += 2; break
 234                                 }
 235
 236                                 if (p > v) {
 237                                         p = n
 238                                 }
 239                         }
 240
 241                         n /= p
 242
 243                         state.factors.push(p)
 244                         state.factors.push(n)
 245                 }
 246
 247                 this.state = state
 248         }
 249
 250         complex.prototype.simple = function (output, input, t) {
 251                 this.process(output, 0, 1, input, 0, 1, t)
 252         }
 253
 254         complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) {
 255                 var outputStride = ~~outputStride, inputStride = ~~inputStride
 256
 257                 var type = t == 'real' ? t : 'complex'
 258
 259                 if (outputStride < 1) {
 260                         throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'")
 261                 }
 262
 263                 if (inputStride < 1) {
 264                         throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'")
 265                 }
 266
 267                 if (type == 'real') {
 268                         for (var i = 0; i < this.state.n; i++) {
 269                                 var x0_r = input[inputOffset + inputStride * i]
 270                                 var x0_i = 0.0
 271
 272                                 this.state.scratch[2 * (i)] = x0_r, this.state.scratch[2 * (i) + 1] = x0_i
 273                         }
 274
 275                         work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state)
 276                 } else {
 277                         if (input == output) {
 278                                 work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state)
 279
 280                                 for (var i = 0; i < this.state.n; i++) {
 281                                         var x0_r = this.state.scratch[2 * (i)], x0_i = this.state.scratch[2 * (i) + 1]
 282
 283                                         output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i
 284                                 }
 285                         } else {
 286                                 work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state)
 287                         }
 288                 }
 289         }
 290
 291         namespace.complex = complex
 292 }(FFT)
 293
 294 module.exports = FFT