/**
 * 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.
 */

export type FFTArrayTypes = number[] | Float32Array | Float64Array;

export class FFT {
  private csize: number;
  private width: number;
  private table: Float64Array;
  private bitrev: Uint32Array;

  private data: FFTArrayTypes | null;
  private out: FFTArrayTypes | null;
  private inv: boolean;

  constructor(readonly size: number) {
    if (size <= 1 || (size & (size - 1)) !== 0)
      throw new Error('FFT size must be a power of two and bigger than 1');

    this.csize = size << 1;

    const table = new Float64Array(size * 2);
    for (var i = 0; i < table.length; i += 2) {
      const angle = Math.PI * i / size;
      table[i] = Math.cos(angle);
      table[i + 1] = -Math.sin(angle);
    }
    this.table = table;

    // Find size's power of two
    var power = 0;
    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
    this.width = power % 2 === 0 ? power - 1 : power;

    // Pre-compute bit-reversal patterns
    const bitrev = new Uint32Array(1 << this.width);
    for (var j = 0; j < bitrev.length; j++) {
      bitrev[j] = 0;
      for (var shift = 0; shift < this.width; shift += 2) {
        var revShift = this.width - shift - 2;
        bitrev[j] |= ((j >>> shift) & 3) << revShift;
      }
    }
    this.bitrev = bitrev;

    this.out = null;
    this.data = null;
    this.inv = false;
  }

  fromComplexArray(complex: FFTArrayTypes, storage?: FFTArrayTypes) {
    const res: FFTArrayTypes = storage || new Float64Array(complex.length >>> 1);
    for (var i = 0; i < complex.length; i += 2)
      res[i >>> 1] = complex[i];
    return res;
  }

  createComplexArray() {
    return new Float64Array(this.csize);
  }

  toComplexArray(input: FFTArrayTypes, storage?: FFTArrayTypes) {
    const res = storage || this.createComplexArray();
    for (var i = 0; i < res.length; i += 2) {
      res[i] = input[i >>> 1];
      res[i + 1] = 0;
    }
    return res;
  }

  completeSpectrum(spectrum: FFTArrayTypes) {
    const size = this.csize;
    const half = size >>> 1;
    for (var i = 2; i < half; i += 2) {
      spectrum[size - i] = spectrum[i];
      spectrum[size - i + 1] = -spectrum[i + 1];
    }
  }

  transform(out: FFTArrayTypes, data: FFTArrayTypes) {
    if (out === data)
      throw new Error('Input and output buffers must be different');
    if (data.length < this.csize)
      throw new Error('Input buffer too small');
    if (out.length < this.csize)
      throw new Error('Output buffer too small');

    this.out = out;
    this.data = data;
    this.inv = false;
    this.transform4();
    this.out = null;
    this.data = null;
  }

  realTransform(out: FFTArrayTypes, data: FFTArrayTypes) {
    if (out === data)
      throw new Error('Input and output buffers must be different');
    if (data.length < this.size)
      throw new Error('Input buffer too small');
    if (out.length < this.csize)
      throw new Error('Output buffer too small');

    this.out = out;
    this.data = data;
    this.inv = false;
    this.realTransform4();
    this.out = null;
    this.data = null;
  }

  inverseTransform(out: FFTArrayTypes, data: FFTArrayTypes) {
    if (out === data)
      throw new Error('Input and output buffers must be different');
    if (data.length < this.csize)
      throw new Error('Input buffer too small');
    if (out.length < this.csize)
      throw new Error('Output buffer too small');

    this.out = out;
    this.data = data;
    this.inv = true;
    this.transform4();
    for (var i = 0; i < out.length; i++)
      out[i] /= this.size;
    this.out = null;
    this.data = null;
  }

  // radix-4 implementation
  private transform4() {
    const out = this.out!;
    const size = this.csize;

    // Initial step (permute and transform)
    var width = this.width;
    var step = 1 << width;
    var len = (size / step) << 1;

    var outOff;
    var t;
    const bitrev = this.bitrev;
    if (len === 4) {
      for (outOff = 0, t = 0; outOff < size; outOff += len, t++) {
        const off = bitrev[t];
        this.singleTransform2(outOff, off, step);
      }
    } else {
      // len === 8
      for (outOff = 0, t = 0; outOff < size; outOff += len, t++) {
        const off = bitrev[t];
        this.singleTransform4(outOff, off, step);
      }
    }

    // Loop through steps in decreasing order
    var inv = this.inv ? -1 : 1;
    var table = this.table;
    for (step >>= 2; step >= 2; step >>= 2) {
      len = (size / step) << 1;
      var quarterLen = len >>> 2;

      // Loop through offsets in the data
      for (outOff = 0; outOff < size; outOff += len) {
        // Full case
        var limit = outOff + quarterLen;
        for (var i = outOff, k = 0; i < limit; i += 2, k += step) {
          const A = i;
          const B = A + quarterLen;
          const C = B + quarterLen;
          const D = C + quarterLen;

          // Original values
          const Ar = out[A];
          const Ai = out[A + 1];
          const Br = out[B];
          const Bi = out[B + 1];
          const Cr = out[C];
          const Ci = out[C + 1];
          const Dr = out[D];
          const Di = out[D + 1];

          // Middle values
          const MAr = Ar;
          const MAi = Ai;

          const tableBr = table[k];
          const tableBi = inv * table[k + 1];
          const MBr = Br * tableBr - Bi * tableBi;
          const MBi = Br * tableBi + Bi * tableBr;

          const tableCr = table[2 * k];
          const tableCi = inv * table[2 * k + 1];
          const MCr = Cr * tableCr - Ci * tableCi;
          const MCi = Cr * tableCi + Ci * tableCr;

          const tableDr = table[3 * k];
          const tableDi = inv * table[3 * k + 1];
          const MDr = Dr * tableDr - Di * tableDi;
          const MDi = Dr * tableDi + Di * tableDr;

          // Pre-Final values
          const T0r = MAr + MCr;
          const T0i = MAi + MCi;
          const T1r = MAr - MCr;
          const T1i = MAi - MCi;
          const T2r = MBr + MDr;
          const T2i = MBi + MDi;
          const T3r = inv * (MBr - MDr);
          const T3i = inv * (MBi - MDi);

          // Final values
          const FAr = T0r + T2r;
          const FAi = T0i + T2i;

          const FCr = T0r - T2r;
          const FCi = T0i - T2i;

          const FBr = T1r + T3i;
          const FBi = T1i - T3r;

          const FDr = T1r - T3i;
          const FDi = T1i + T3r;

          out[A] = FAr;
          out[A + 1] = FAi;
          out[B] = FBr;
          out[B + 1] = FBi;
          out[C] = FCr;
          out[C + 1] = FCi;
          out[D] = FDr;
          out[D + 1] = FDi;
        }
      }
    }
  }

  // radix-2 implementation
  //
  // NOTE: Only called for len=4
  private singleTransform2(outOff: number, off: number, step: number) {
    const out = this.out!;
    const data = this.data!;

    const evenR = data[off];
    const evenI = data[off + 1];
    const oddR = data[off + step];
    const oddI = data[off + step + 1];

    const leftR = evenR + oddR;
    const leftI = evenI + oddI;
    const rightR = evenR - oddR;
    const rightI = evenI - oddI;

    out[outOff] = leftR;
    out[outOff + 1] = leftI;
    out[outOff + 2] = rightR;
    out[outOff + 3] = rightI;
  }

  // radix-4
  //
  // NOTE: Only called for len=8
  private singleTransform4(outOff: number, off: number, step: number) {
    const out = this.out!;
    const data = this.data!;
    const inv = this.inv ? -1 : 1;
    const step2 = step * 2;
    const step3 = step * 3;

    // Original values
    const Ar = data[off];
    const Ai = data[off + 1];
    const Br = data[off + step];
    const Bi = data[off + step + 1];
    const Cr = data[off + step2];
    const Ci = data[off + step2 + 1];
    const Dr = data[off + step3];
    const Di = data[off + step3 + 1];

    // Pre-Final values
    const T0r = Ar + Cr;
    const T0i = Ai + Ci;
    const T1r = Ar - Cr;
    const T1i = Ai - Ci;
    const T2r = Br + Dr;
    const T2i = Bi + Di;
    const T3r = inv * (Br - Dr);
    const T3i = inv * (Bi - Di);

    // Final values
    const FAr = T0r + T2r;
    const FAi = T0i + T2i;

    const FBr = T1r + T3i;
    const FBi = T1i - T3r;

    const FCr = T0r - T2r;
    const FCi = T0i - T2i;

    const FDr = T1r - T3i;
    const FDi = T1i + T3r;

    out[outOff] = FAr;
    out[outOff + 1] = FAi;
    out[outOff + 2] = FBr;
    out[outOff + 3] = FBi;
    out[outOff + 4] = FCr;
    out[outOff + 5] = FCi;
    out[outOff + 6] = FDr;
    out[outOff + 7] = FDi;
  };

  // Real input radix-4 implementation
  private realTransform4() {
    const out = this.out!;
    var size = this.csize;

    // Initial step (permute and transform)
    var width = this.width;
    var step = 1 << width;
    var len = (size / step) << 1;

    const bitrev = this.bitrev;
    if (len === 4) {
      for (let outOff = 0, t = 0; outOff < size; outOff += len, t++) {
        const off = bitrev[t];
        this.singleRealTransform2(outOff, off >>> 1, step >>> 1);
      }
    } else {
      // len === 8
      for (let outOff = 0, t = 0; outOff < size; outOff += len, t++) {
        const off = bitrev[t];
        this.singleRealTransform4(outOff, off >>> 1, step >>> 1);
      }
    }

    // Loop through steps in decreasing order
    const inv = this.inv ? -1 : 1;
    var table = this.table;
    for (step >>= 2; step >= 2; step >>= 2) {
      len = (size / step) << 1;
      var halfLen = len >>> 1;
      var quarterLen = halfLen >>> 1;
      var hquarterLen = quarterLen >>> 1;

      // Loop through offsets in the data
      for (let outOff = 0; outOff < size; outOff += len) {
        for (var i = 0, k = 0; i <= hquarterLen; i += 2, k += step) {
          var A = outOff + i;
          var B = A + quarterLen;
          var C = B + quarterLen;
          var D = C + quarterLen;

          // Original values
          var Ar = out[A];
          var Ai = out[A + 1];
          var Br = out[B];
          var Bi = out[B + 1];
          var Cr = out[C];
          var Ci = out[C + 1];
          var Dr = out[D];
          var Di = out[D + 1];

          // Middle values
          var MAr = Ar;
          var MAi = Ai;

          var tableBr = table[k];
          var tableBi = inv * table[k + 1];
          var MBr = Br * tableBr - Bi * tableBi;
          var MBi = Br * tableBi + Bi * tableBr;

          var tableCr = table[2 * k];
          var tableCi = inv * table[2 * k + 1];
          var MCr = Cr * tableCr - Ci * tableCi;
          var MCi = Cr * tableCi + Ci * tableCr;

          var tableDr = table[3 * k];
          var tableDi = inv * table[3 * k + 1];
          var MDr = Dr * tableDr - Di * tableDi;
          var MDi = Dr * tableDi + Di * tableDr;

          // Pre-Final values
          var T0r = MAr + MCr;
          var T0i = MAi + MCi;
          var T1r = MAr - MCr;
          var T1i = MAi - MCi;
          var T2r = MBr + MDr;
          var T2i = MBi + MDi;
          var T3r = inv * (MBr - MDr);
          var T3i = inv * (MBi - MDi);

          // Final values
          var FAr = T0r + T2r;
          var FAi = T0i + T2i;

          var FBr = T1r + T3i;
          var FBi = T1i - T3r;

          out[A] = FAr;
          out[A + 1] = FAi;
          out[B] = FBr;
          out[B + 1] = FBi;

          // Output final middle point
          if (i === 0) {
            var FCr = T0r - T2r;
            var FCi = T0i - T2i;
            out[C] = FCr;
            out[C + 1] = FCi;
            continue;
          }

          // Do not overwrite ourselves
          if (i === hquarterLen)
            continue;

          // In the flipped case:
          // MAi = -MAi
          // MBr=-MBi, MBi=-MBr
          // MCr=-MCr
          // MDr=MDi, MDi=MDr
          var ST0r = T1r;
          var ST0i = -T1i;
          var ST1r = T0r;
          var ST1i = -T0i;
          var ST2r = -inv * T3i;
          var ST2i = -inv * T3r;
          var ST3r = -inv * T2i;
          var ST3i = -inv * T2r;

          var SFAr = ST0r + ST2r;
          var SFAi = ST0i + ST2i;

          var SFBr = ST1r + ST3i;
          var SFBi = ST1i - ST3r;

          var SA = outOff + quarterLen - i;
          var SB = outOff + halfLen - i;

          out[SA] = SFAr;
          out[SA + 1] = SFAi;
          out[SB] = SFBr;
          out[SB + 1] = SFBi;
        }
      }
    }
  }

  // radix-2 implementation
  //
  // NOTE: Only called for len=4
  private singleRealTransform2(outOff: number, off: number, step: number) {
    const out = this.out!;
    const data = this.data!;

    const evenR = data[off];
    const oddR = data[off + step];

    const leftR = evenR + oddR;
    const rightR = evenR - oddR;

    out[outOff] = leftR;
    out[outOff + 1] = 0;
    out[outOff + 2] = rightR;
    out[outOff + 3] = 0;
  }

  // radix-4
  //
  // NOTE: Only called for len=8
  private singleRealTransform4(outOff: number, off: number, step: number) {
    const out = this.out!;
    const data = this.data!;
    const inv = this.inv ? -1 : 1;
    const step2 = step * 2;
    const step3 = step * 3;

    // Original values
    const Ar = data[off];
    const Br = data[off + step];
    const Cr = data[off + step2];
    const Dr = data[off + step3];

    // Pre-Final values
    const T0r = Ar + Cr;
    const T1r = Ar - Cr;
    const T2r = Br + Dr;
    const T3r = inv * (Br - Dr);

    // Final values
    const FAr = T0r + T2r;

    const FBr = T1r;
    const FBi = -T3r;

    const FCr = T0r - T2r;

    const FDr = T1r;
    const FDi = T3r;

    out[outOff] = FAr;
    out[outOff + 1] = 0;
    out[outOff + 2] = FBr;
    out[outOff + 3] = FBi;
    out[outOff + 4] = FCr;
    out[outOff + 5] = 0;
    out[outOff + 6] = FDr;
    out[outOff + 7] = FDi;
  }
}