/*===========================================================================*\
 * Fast Fourier Transform (Cooley-Tukey Method)
 *
 * (c) Vail Systems. Joshua Jung and Ben Bryan. 2015
 *
 * This code is not designed to be highly optimized but as an educational
 * tool to understand the Fast Fourier Transform.
 *
 * Can be used with Vectors (typed arrays), Arrays
 *
 *
\*===========================================================================*/
var Vector = Require('numerics/vector');
var Matrix = Require('numerics/matrix');

//------------------------------------------------
// Note: Some of this code is not optimized and is
// primarily designed as an educational and testing
// tool.
// To get high performace would require transforming
// the recursive calls into a loop and then loop
// unrolling. All of this is best accomplished
// in C or assembly.
//-------------------------------------------------

//-------------------------------------------------
// The following code assumes a complex number is
// an array: [real, imaginary]
//-------------------------------------------------
var complex = {
  //-------------------------------------------------
  // Add two complex numbers
  //-------------------------------------------------
  add : function (a, b)
  {
      return [a[0] + b[0], a[1] + b[1]];
  },

  //-------------------------------------------------
  // Subtract two complex numbers
  //-------------------------------------------------
  subtract : function (a, b)
  {
      return [a[0] - b[0], a[1] - b[1]];
  },

  //-------------------------------------------------
  // Multiply two complex numbers
  //
  // (a + bi) * (c + di) = (ac - bd) + (ad + bc)i
  //-------------------------------------------------
  multiply : function (a, b) 
  {
      return [(a[0] * b[0] - a[1] * b[1]), 
              (a[0] * b[1] + a[1] * b[0])];
  },

  //-------------------------------------------------
  // Calculate |a + bi|
  //
  // sqrt(a*a + b*b)
  //-------------------------------------------------
  magnitude : function (c) 
  {
      return Math.sqrt(c[0]*c[0] + c[1]*c[1]); 
  },
  
  phase : function (c) 
  {
      return c[0]!=0?Math.atan(c[1]/c[0])*180/Math.PI:(c[1]>0?90:-90); 
  }

}
var mapExponent = {};
var fftUtil = {
  //-------------------------------------------------
  // By Eulers Formula:
  //
  // e^(i*x) = cos(x) + i*sin(x)
  //
  // and in DFT:
  //
  // x = -2*PI*(k/N)
  //-------------------------------------------------
  exponent : function (k, N) {
        var x = -2 * Math.PI * (k / N);

        mapExponent[N] = mapExponent[N] || {};
        mapExponent[N][k] = mapExponent[N][k] || [Math.cos(x), Math.sin(x)];// [Real, Imaginary]

        return mapExponent[N][k];
  },

  //-------------------------------------------------
  // Calculate FFT Magnitude for complex numbers.
  //-------------------------------------------------
  fftMag : function (fftBins) {
    if (isArray(fftBins)) {
      var ret = fftBins.map(complex.magnitude);
      return ret.slice(0, ret.length / 2);
    } else if (isVector(fftBins) || isMatrix(fftBins)) {
      // complex matrix (2 columns)
      if (fftBins.columns != 2) throw "fft.fftMag: Complex matrix columns != 2";
      var ret = Vector(fftBins.rows,{dtn:fftBins.dtn})
      return ret.eval(function (v,i) { 
       return complex.magnitude([fftBins.get(i,0),fftBins.get(i,1)]) }); 
    }
  },
  fftPha : function (fftBins) {
    if (isArray(fftBins)) {
      var ret = fftBins.map(complex.phase);
      return ret.slice(0, ret.length / 2);
    } else if (isVector(fftBins) || isMatrix(fftBins)) {
      // complex matrix (2 columns)
      if (fftBins.columns != 2) throw "fft.fftMag: Complex matrix columns != 2";
      var ret = Vector(fftBins.rows,{dtn:fftBins.dtn})
      return ret.eval(function (v,i) { 
       return complex.phase([fftBins.get(i,0),fftBins.get(i,1)]) }); 
    }
  },

  //-------------------------------------------------
  // Calculate Frequency Bins
  // 
  // Returns an array of the frequencies (in hertz) of
  // each FFT bin provided, assuming the sampleRate is
  // samples taken per second.
  //-------------------------------------------------
  fftFreq : function (fftBins, sampleRate) {
      var stepFreq = sampleRate / (fftBins.length);
      var ret = fftBins.slice(0, fftBins.length / 2);

      return ret.map(function (__, ix) {
          return ix * stepFreq;
      });
  }
}
    
// Bit-twiddle
var REVERSE_TABLE = new Array(256);
var INT_BITS = 32;  //Number of bits in an integer

(function(tab) {
  for(var i=0; i<256; ++i) {
    var v = i, r = i, s = 7;
    for (v >>>= 1; v; v >>>= 1) {
      r <<= 1;
      r |= v & 1;
      --s;
    }
    tab[i] = (r << s) & 0xff;
  }
})(REVERSE_TABLE);

var twiddle = {

  //Constants
  INT_BITS  : INT_BITS,
  INT_MAX   :  0x7fffffff,
  INT_MIN   : -1<<(INT_BITS-1),

  //Returns -1, 0, +1 depending on sign of x
  sign : function(v) {
    return (v > 0) - (v < 0);
  },

  //Computes absolute value of integer
  abs : function(v) {
    var mask = v >> (INT_BITS-1);
    return (v ^ mask) - mask;
  },

  //Computes minimum of integers x and y
  min : function(x, y) {
    return y ^ ((x ^ y) & -(x < y));
  },

  //Computes maximum of integers x and y
  max : function(x, y) {
    return x ^ ((x ^ y) & -(x < y));
  },

  //Checks if a number is a power of two
  isPow2 : function(v) {
    return !(v & (v-1)) && (!!v);
  },

  //Computes log base 2 of v
  log2 : function(v) {
    var r, shift;
    r =     (v > 0xFFFF) << 4; v >>>= r;
    shift = (v > 0xFF  ) << 3; v >>>= shift; r |= shift;
    shift = (v > 0xF   ) << 2; v >>>= shift; r |= shift;
    shift = (v > 0x3   ) << 1; v >>>= shift; r |= shift;
    return r | (v >> 1);
  },

  //Computes log base 10 of v
  log10 : function(v) {
    return  (v >= 1000000000) ? 9 : (v >= 100000000) ? 8 : (v >= 10000000) ? 7 :
            (v >= 1000000) ? 6 : (v >= 100000) ? 5 : (v >= 10000) ? 4 :
            (v >= 1000) ? 3 : (v >= 100) ? 2 : (v >= 10) ? 1 : 0;
  },

  //Counts number of bits
  popCount : function(v) {
    v = v - ((v >>> 1) & 0x55555555);
    v = (v & 0x33333333) + ((v >>> 2) & 0x33333333);
    return ((v + (v >>> 4) & 0xF0F0F0F) * 0x1010101) >>> 24;
  },

  //Counts number of trailing zeros
  countTrailingZeros : function (v) {
    var c = 32;
    v &= -v;
    if (v) c--;
    if (v & 0x0000FFFF) c -= 16;
    if (v & 0x00FF00FF) c -= 8;
    if (v & 0x0F0F0F0F) c -= 4;
    if (v & 0x33333333) c -= 2;
    if (v & 0x55555555) c -= 1;
    return c;
  },

  //Rounds to next power of 2
  nextPow2 : function(v) {
    v += v === 0;
    --v;
    v |= v >>> 1;
    v |= v >>> 2;
    v |= v >>> 4;
    v |= v >>> 8;
    v |= v >>> 16;
    return v + 1;
  },

  //Rounds down to previous power of 2
  prevPow2 : function(v) {
    v |= v >>> 1;
    v |= v >>> 2;
    v |= v >>> 4;
    v |= v >>> 8;
    v |= v >>> 16;
    return v - (v>>>1);
  },

  //Computes parity of word
  parity : function(v) {
    v ^= v >>> 16;
    v ^= v >>> 8;
    v ^= v >>> 4;
    v &= 0xf;
    return (0x6996 >>> v) & 1;
  },


  //Reverse bits in a 32 bit word
  reverse : function(v) {
    return  (REVERSE_TABLE[ v         & 0xff] << 24) |
            (REVERSE_TABLE[(v >>> 8)  & 0xff] << 16) |
            (REVERSE_TABLE[(v >>> 16) & 0xff] << 8)  |
             REVERSE_TABLE[(v >>> 24) & 0xff];
  },

  //Interleave bits of 2 coordinates with 16 bits.  Useful for fast quadtree codes
  interleave2 : function(x, y) {
    x &= 0xFFFF;
    x = (x | (x << 8)) & 0x00FF00FF;
    x = (x | (x << 4)) & 0x0F0F0F0F;
    x = (x | (x << 2)) & 0x33333333;
    x = (x | (x << 1)) & 0x55555555;

    y &= 0xFFFF;
    y = (y | (y << 8)) & 0x00FF00FF;
    y = (y | (y << 4)) & 0x0F0F0F0F;
    y = (y | (y << 2)) & 0x33333333;
    y = (y | (y << 1)) & 0x55555555;

    return x | (y << 1);
  },

  //Extracts the nth interleaved component
  deinterleave2 : function(v, n) {
    v = (v >>> n) & 0x55555555;
    v = (v | (v >>> 1))  & 0x33333333;
    v = (v | (v >>> 2))  & 0x0F0F0F0F;
    v = (v | (v >>> 4))  & 0x00FF00FF;
    v = (v | (v >>> 16)) & 0x000FFFF;
    return (v << 16) >> 16;
  },


  //Interleave bits of 3 coordinates, each with 10 bits.  Useful for fast octree codes
  interleave3 : function(x, y, z) {
    x &= 0x3FF;
    x  = (x | (x<<16)) & 4278190335;
    x  = (x | (x<<8))  & 251719695;
    x  = (x | (x<<4))  & 3272356035;
    x  = (x | (x<<2))  & 1227133513;

    y &= 0x3FF;
    y  = (y | (y<<16)) & 4278190335;
    y  = (y | (y<<8))  & 251719695;
    y  = (y | (y<<4))  & 3272356035;
    y  = (y | (y<<2))  & 1227133513;
    x |= (y << 1);

    z &= 0x3FF;
    z  = (z | (z<<16)) & 4278190335;
    z  = (z | (z<<8))  & 251719695;
    z  = (z | (z<<4))  & 3272356035;
    z  = (z | (z<<2))  & 1227133513;

    return x | (z << 2);
  },

  //Extracts nth interleaved component of a 3-tuple
  deinterleave3 : function(v, n) {
    v = (v >>> n)       & 1227133513;
    v = (v | (v>>>2))   & 3272356035;
    v = (v | (v>>>4))   & 251719695;
    v = (v | (v>>>8))   & 4278190335;
    v = (v | (v>>>16))  & 0x3FF;
    return (v<<22)>>22;
  },

  //Computes next combination in colexicographic order (this is mistakenly called nextPermutation on the bit twiddling hacks page)
  nextCombination : function(v) {
    var t = v | (v - 1);
    return (t + 1) | (((~t & -~t) - 1) >>> (twiddle.countTrailingZeros(v) + 1));
  }

}

function checkpow2(info,vector) {
  if (Math.floor(Math.log2(vector.length)) != Math.log2(vector.length))
    throw ('fft.'+info+' error: vector length must be a power of 2')
}

module.exports = {
  //-------------------------------------------------
  // Calculate FFT for vector where vector.length
  // is assumed to be a power of 2.
  //-------------------------------------------------
  fft: function fft(vector) {
    if (vector.data) vector=vector.data; // Matrix|Vector
    
    checkpow2('fft',vector);
    
    var X = [],
        N = vector.length;

    // Base case is X = x + 0i since our input is assumed to be real only.
    if (N == 1) {
      if (Array.isArray(vector[0])) //If input vector contains complex numbers
        return [[vector[0][0], vector[0][1]]];      
      else
        return [[vector[0], 0]];
    }

    // Recurse: all even samples
    var X_evens = fft(vector.filter(even)),

        // Recurse: all odd samples
        X_odds  = fft(vector.filter(odd));

    // Now, perform N/2 operations!
    for (var k = 0; k < N / 2; k++) {
      // t is a complex number!
      var t = X_evens[k],
          e = complex.multiply(fftUtil.exponent(k, N), X_odds[k]);

      X[k] = complex.add(t, e);
      X[k + (N / 2)] = complex.subtract(t, e);
    }

    function even(__, ix) {
      return ix % 2 == 0;
    }

    function odd(__, ix) {
      return ix % 2 == 1;
    }

    return X;
  },
  //-------------------------------------------------
  // Calculate FFT for vector where vector.length
  // is assumed to be a power of 2.  This is the in-
  // place implementation, to avoid the memory
  // footprint used by recursion.
  //-------------------------------------------------
  fftInPlace: function(vector) {
    if (vector.data) vector=vector.data; // Matrix|Vector
    checkpow2('fftInPlace',vector);
    
    var N = vector.length;

    var trailingZeros = twiddle.countTrailingZeros(N); //Once reversed, this will be leading zeros

    // Reverse bits
    for (var k = 0; k < N; k++) {
      var p = twiddle.reverse(k) >>> (twiddle.INT_BITS - trailingZeros);
      if (p > k) {
        var complexTemp = [vector[k], 0];
        vector[k] = vector[p];
        vector[p] = complexTemp;
      } else {
        vector[p] = [vector[p], 0];
      }
    }

    //Do the DIT now in-place
    for (var len = 2; len <= N; len += len) {
      for (var i = 0; i < len / 2; i++) {
        var w = fftUtil.exponent(i, len);
        for (var j = 0; j < N / len; j++) {
          var t = complex.multiply(w, vector[j * len + i + len / 2]);
          vector[j * len + i + len / 2] = complex.subtract(vector[j * len + i], t);
          vector[j * len + i] = complex.add(vector[j * len + i], t);
        }
      }
    }
  },
  ifft: function ifft(signal){
    if (signal.data) signal=signal.data; // Matrix
    checkpow2('ifft',signal);
    //Interchange real and imaginary parts
    var csignal=[];
    for(var i=0; i<signal.length; i++){
        csignal[i]=[signal[i][1], signal[i][0]];
    }

    //Apply fft
    var ps=module.exports.fft(csignal);

    //Interchange real and imaginary parts and normalize
    var res=[];
    for(var j=0; j<ps.length; j++){
        res[j]=[ps[j][1]/ps.length, ps[j][0]/ps.length];
    }
    return res;
  },
  fftMag: fftUtil.fftMag,
  fftPha: fftUtil.fftPha,
  fftFreq: fftUtil.fftFreq,
};