jevstrudel.git/packages/superdough/fft.js

jevstrudel.git / packages / superdough / fft.js

fft.jssource488 lines · 13.8 KB · raw

1'use strict';

sourced from https://github.com/indutny/fft.js/ LICENSE This software is licensed under the MIT License. Copyright Fedor Indutny, 2017. 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.

10export default class FFT {
11  constructor(size) {
12    this.size = size | 0;
13    if (this.size <= 1 || (this.size & (this.size - 1)) !== 0)
14      throw new Error('FFT size must be a power of two and bigger than 1');
15
16    this._csize = size << 1;

NOTE: Use of var is intentional for old V8 versions

19    var table = new Array(this.size * 2);
20    for (var i = 0; i < table.length; i += 2) {
21      const angle = (Math.PI * i) / this.size;
22      table[i] = Math.cos(angle);
23      table[i + 1] = -Math.sin(angle);
24    }
25    this.table = table;

Find size's power of two

28    var power = 0;
29    for (var t = 1; this.size > t; t <<= 1) power++;

Calculate initial step's width:

If we are full radix-4 - it is 2x smaller to give inital len=8
Otherwise it is the same as power to give len=4

34    this._width = power % 2 === 0 ? power - 1 : power;

Pre-compute bit-reversal patterns

37    this._bitrev = new Array(1 << this._width);
38    for (var j = 0; j < this._bitrev.length; j++) {
39      this._bitrev[j] = 0;
40      for (var shift = 0; shift < this._width; shift += 2) {
41        var revShift = this._width - shift - 2;
42        this._bitrev[j] |= ((j >>> shift) & 3) << revShift;
43      }
44    }

46    this._out = null;
47    this._data = null;
48    this._inv = 0;
49  }
50  fromComplexArray(complex, storage) {
51    var res = storage || new Array(complex.length >>> 1);
52    for (var i = 0; i < complex.length; i += 2) res[i >>> 1] = complex[i];
53    return res;
54  }
55  createComplexArray() {
56    const res = new Array(this._csize);
57    for (var i = 0; i < res.length; i++) res[i] = 0;
58    return res;
59  }
60  toComplexArray(input, storage) {
61    var res = storage || this.createComplexArray();
62    for (var i = 0; i < res.length; i += 2) {
63      res[i] = input[i >>> 1];
64      res[i + 1] = 0;
65    }
66    return res;
67  }
68  completeSpectrum(spectrum) {
69    var size = this._csize;
70    var half = size >>> 1;
71    for (var i = 2; i < half; i += 2) {
72      spectrum[size - i] = spectrum[i];
73      spectrum[size - i + 1] = -spectrum[i + 1];
74    }
75  }
76  transform(out, data) {
77    if (out === data) throw new Error('Input and output buffers must be different');
78
79    this._out = out;
80    this._data = data;
81    this._inv = 0;
82    this._transform4();
83    this._out = null;
84    this._data = null;
85  }
86  realTransform(out, data) {
87    if (out === data) throw new Error('Input and output buffers must be different');
88
89    this._out = out;
90    this._data = data;
91    this._inv = 0;
92    this._realTransform4();
93    this._out = null;
94    this._data = null;
95  }
96  inverseTransform(out, data) {
97    if (out === data) throw new Error('Input and output buffers must be different');
98
99    this._out = out;
100    this._data = data;
101    this._inv = 1;
102    this._transform4();
103    for (var i = 0; i < out.length; i++) out[i] /= this.size;
104    this._out = null;
105    this._data = null;
106  }
107  // radix-4 implementation
108  //
109  // NOTE: Uses of `var` are intentional for older V8 version that do not
110  // support both `let compound assignments` and `const phi`
111  _transform4() {
112    var out = this._out;
113    var size = this._csize;

Initial step (permute and transform)

116    var width = this._width;
117    var step = 1 << width;
118    var len = (size / step) << 1;

120    var outOff;
121    var t;
122    var bitrev = this._bitrev;
123    if (len === 4) {
124      for (outOff = 0, t = 0; outOff < size; outOff += len, t++) {
125        const off = bitrev[t];
126        this._singleTransform2(outOff, off, step);
127      }
128    } else {
129      // len === 8
130      for (outOff = 0, t = 0; outOff < size; outOff += len, t++) {
131        const off = bitrev[t];
132        this._singleTransform4(outOff, off, step);
133      }
134    }

Loop through steps in decreasing order

137    var inv = this._inv ? -1 : 1;
138    var table = this.table;
139    for (step >>= 2; step >= 2; step >>= 2) {
140      len = (size / step) << 1;
141      var quarterLen = len >>> 2;

Loop through offsets in the data

144      for (outOff = 0; outOff < size; outOff += len) {
145        // Full case
146        var limit = outOff + quarterLen;
147        for (var i = outOff, k = 0; i < limit; i += 2, k += step) {
148          const A = i;
149          const B = A + quarterLen;
150          const C = B + quarterLen;
151          const D = C + quarterLen;

Original values

154          const Ar = out[A];
155          const Ai = out[A + 1];
156          const Br = out[B];
157          const Bi = out[B + 1];
158          const Cr = out[C];
159          const Ci = out[C + 1];
160          const Dr = out[D];
161          const Di = out[D + 1];

Middle values

164          const MAr = Ar;
165          const MAi = Ai;

167          const tableBr = table[k];
168          const tableBi = inv * table[k + 1];
169          const MBr = Br * tableBr - Bi * tableBi;
170          const MBi = Br * tableBi + Bi * tableBr;
171
172          const tableCr = table[2 * k];
173          const tableCi = inv * table[2 * k + 1];
174          const MCr = Cr * tableCr - Ci * tableCi;
175          const MCi = Cr * tableCi + Ci * tableCr;
176
177          const tableDr = table[3 * k];
178          const tableDi = inv * table[3 * k + 1];
179          const MDr = Dr * tableDr - Di * tableDi;
180          const MDi = Dr * tableDi + Di * tableDr;

Pre-Final values

183          const T0r = MAr + MCr;
184          const T0i = MAi + MCi;
185          const T1r = MAr - MCr;
186          const T1i = MAi - MCi;
187          const T2r = MBr + MDr;
188          const T2i = MBi + MDi;
189          const T3r = inv * (MBr - MDr);
190          const T3i = inv * (MBi - MDi);

Final values

193          const FAr = T0r + T2r;
194          const FAi = T0i + T2i;

196          const FCr = T0r - T2r;
197          const FCi = T0i - T2i;
198
199          const FBr = T1r + T3i;
200          const FBi = T1i - T3r;
201
202          const FDr = T1r - T3i;
203          const FDi = T1i + T3r;
204
205          out[A] = FAr;
206          out[A + 1] = FAi;
207          out[B] = FBr;
208          out[B + 1] = FBi;
209          out[C] = FCr;
210          out[C + 1] = FCi;
211          out[D] = FDr;
212          out[D + 1] = FDi;
213        }
214      }
215    }
216  }
217  // radix-2 implementation
218  //
219  // NOTE: Only called for len=4
220  _singleTransform2(outOff, off, step) {
221    const out = this._out;
222    const data = this._data;
223
224    const evenR = data[off];
225    const evenI = data[off + 1];
226    const oddR = data[off + step];
227    const oddI = data[off + step + 1];
228
229    const leftR = evenR + oddR;
230    const leftI = evenI + oddI;
231    const rightR = evenR - oddR;
232    const rightI = evenI - oddI;
233
234    out[outOff] = leftR;
235    out[outOff + 1] = leftI;
236    out[outOff + 2] = rightR;
237    out[outOff + 3] = rightI;
238  }
239  // radix-4
240  //
241  // NOTE: Only called for len=8
242  _singleTransform4(outOff, off, step) {
243    const out = this._out;
244    const data = this._data;
245    const inv = this._inv ? -1 : 1;
246    const step2 = step * 2;
247    const step3 = step * 3;

Original values

250    const Ar = data[off];
251    const Ai = data[off + 1];
252    const Br = data[off + step];
253    const Bi = data[off + step + 1];
254    const Cr = data[off + step2];
255    const Ci = data[off + step2 + 1];
256    const Dr = data[off + step3];
257    const Di = data[off + step3 + 1];

Pre-Final values

260    const T0r = Ar + Cr;
261    const T0i = Ai + Ci;
262    const T1r = Ar - Cr;
263    const T1i = Ai - Ci;
264    const T2r = Br + Dr;
265    const T2i = Bi + Di;
266    const T3r = inv * (Br - Dr);
267    const T3i = inv * (Bi - Di);

Final values

270    const FAr = T0r + T2r;
271    const FAi = T0i + T2i;

273    const FBr = T1r + T3i;
274    const FBi = T1i - T3r;
275
276    const FCr = T0r - T2r;
277    const FCi = T0i - T2i;
278
279    const FDr = T1r - T3i;
280    const FDi = T1i + T3r;
281
282    out[outOff] = FAr;
283    out[outOff + 1] = FAi;
284    out[outOff + 2] = FBr;
285    out[outOff + 3] = FBi;
286    out[outOff + 4] = FCr;
287    out[outOff + 5] = FCi;
288    out[outOff + 6] = FDr;
289    out[outOff + 7] = FDi;
290  }
291  // Real input radix-4 implementation
292  _realTransform4() {
293    var out = this._out;
294    var size = this._csize;

Initial step (permute and transform)

297    var width = this._width;
298    var step = 1 << width;
299    var len = (size / step) << 1;

301    var outOff;
302    var t;
303    var bitrev = this._bitrev;
304    if (len === 4) {
305      for (outOff = 0, t = 0; outOff < size; outOff += len, t++) {
306        const off = bitrev[t];
307        this._singleRealTransform2(outOff, off >>> 1, step >>> 1);
308      }
309    } else {
310      // len === 8
311      for (outOff = 0, t = 0; outOff < size; outOff += len, t++) {
312        const off = bitrev[t];
313        this._singleRealTransform4(outOff, off >>> 1, step >>> 1);
314      }
315    }

Loop through steps in decreasing order

318    var inv = this._inv ? -1 : 1;
319    var table = this.table;
320    for (step >>= 2; step >= 2; step >>= 2) {
321      len = (size / step) << 1;
322      var halfLen = len >>> 1;
323      var quarterLen = halfLen >>> 1;
324      var hquarterLen = quarterLen >>> 1;

Loop through offsets in the data

327      for (outOff = 0; outOff < size; outOff += len) {
328        for (var i = 0, k = 0; i <= hquarterLen; i += 2, k += step) {
329          var A = outOff + i;
330          var B = A + quarterLen;
331          var C = B + quarterLen;
332          var D = C + quarterLen;

Original values

335          var Ar = out[A];
336          var Ai = out[A + 1];
337          var Br = out[B];
338          var Bi = out[B + 1];
339          var Cr = out[C];
340          var Ci = out[C + 1];
341          var Dr = out[D];
342          var Di = out[D + 1];

Middle values

345          var MAr = Ar;
346          var MAi = Ai;

348          var tableBr = table[k];
349          var tableBi = inv * table[k + 1];
350          var MBr = Br * tableBr - Bi * tableBi;
351          var MBi = Br * tableBi + Bi * tableBr;
352
353          var tableCr = table[2 * k];
354          var tableCi = inv * table[2 * k + 1];
355          var MCr = Cr * tableCr - Ci * tableCi;
356          var MCi = Cr * tableCi + Ci * tableCr;
357
358          var tableDr = table[3 * k];
359          var tableDi = inv * table[3 * k + 1];
360          var MDr = Dr * tableDr - Di * tableDi;
361          var MDi = Dr * tableDi + Di * tableDr;

Pre-Final values

364          var T0r = MAr + MCr;
365          var T0i = MAi + MCi;
366          var T1r = MAr - MCr;
367          var T1i = MAi - MCi;
368          var T2r = MBr + MDr;
369          var T2i = MBi + MDi;
370          var T3r = inv * (MBr - MDr);
371          var T3i = inv * (MBi - MDi);

Final values

374          var FAr = T0r + T2r;
375          var FAi = T0i + T2i;

377          var FBr = T1r + T3i;
378          var FBi = T1i - T3r;
379
380          out[A] = FAr;
381          out[A + 1] = FAi;
382          out[B] = FBr;
383          out[B + 1] = FBi;

Output final middle point

386          if (i === 0) {
387            var FCr = T0r - T2r;
388            var FCi = T0i - T2i;
389            out[C] = FCr;
390            out[C + 1] = FCi;
391            continue;
392          }

Do not overwrite ourselves

395          if (i === hquarterLen) continue;

In the flipped case: MAi = -MAi MBr=-MBi, MBi=-MBr MCr=-MCr MDr=MDi, MDi=MDr

402          var ST0r = T1r;
403          var ST0i = -T1i;
404          var ST1r = T0r;
405          var ST1i = -T0i;
406          var ST2r = -inv * T3i;
407          var ST2i = -inv * T3r;
408          var ST3r = -inv * T2i;
409          var ST3i = -inv * T2r;

411          var SFAr = ST0r + ST2r;
412          var SFAi = ST0i + ST2i;
413
414          var SFBr = ST1r + ST3i;
415          var SFBi = ST1i - ST3r;
416
417          var SA = outOff + quarterLen - i;
418          var SB = outOff + halfLen - i;
419
420          out[SA] = SFAr;
421          out[SA + 1] = SFAi;
422          out[SB] = SFBr;
423          out[SB + 1] = SFBi;
424        }
425      }
426    }
427  }
428  // radix-2 implementation
429  //
430  // NOTE: Only called for len=4
431  _singleRealTransform2(outOff, off, step) {
432    const out = this._out;
433    const data = this._data;
434
435    const evenR = data[off];
436    const oddR = data[off + step];
437
438    const leftR = evenR + oddR;
439    const rightR = evenR - oddR;
440
441    out[outOff] = leftR;
442    out[outOff + 1] = 0;
443    out[outOff + 2] = rightR;
444    out[outOff + 3] = 0;
445  }
446  // radix-4
447  //
448  // NOTE: Only called for len=8
449  _singleRealTransform4(outOff, off, step) {
450    const out = this._out;
451    const data = this._data;
452    const inv = this._inv ? -1 : 1;
453    const step2 = step * 2;
454    const step3 = step * 3;

Original values

457    const Ar = data[off];
458    const Br = data[off + step];
459    const Cr = data[off + step2];
460    const Dr = data[off + step3];

Pre-Final values

463    const T0r = Ar + Cr;
464    const T1r = Ar - Cr;
465    const T2r = Br + Dr;
466    const T3r = inv * (Br - Dr);

Final values

469    const FAr = T0r + T2r;

471    const FBr = T1r;
472    const FBi = -T3r;
473
474    const FCr = T0r - T2r;
475
476    const FDr = T1r;
477    const FDi = T3r;
478
479    out[outOff] = FAr;
480    out[outOff + 1] = 0;
481    out[outOff + 2] = FBr;
482    out[outOff + 3] = FBi;
483    out[outOff + 4] = FCr;
484    out[outOff + 5] = 0;
485    out[outOff + 6] = FDr;
486    out[outOff + 7] = FDi;
487  }
488}