Chris@25: 'use strict';
Chris@25: 
Chris@25: !function(exports, complex_array) {
Chris@25: 
Chris@25:   var
Chris@25:     ComplexArray = complex_array.ComplexArray,
Chris@25:     // Math constants and functions we need.
Chris@25:     PI = Math.PI,
Chris@25:     SQRT1_2 = Math.SQRT1_2,
Chris@25:     sqrt = Math.sqrt,
Chris@25:     cos = Math.cos,
Chris@25:     sin = Math.sin
Chris@25: 
Chris@25:   ComplexArray.prototype.FFT = function() {
Chris@25:     return FFT(this, false)
Chris@25:   }
Chris@25: 
Chris@25:   exports.FFT = function(input) {
Chris@25:     return ensureComplexArray(input).FFT()
Chris@25:   }
Chris@25: 
Chris@25:   ComplexArray.prototype.InvFFT = function() {
Chris@25:     return FFT(this, true)
Chris@25:   }
Chris@25: 
Chris@25:   exports.InvFFT = function(input) {
Chris@25:     return ensureComplexArray(input).InvFFT()
Chris@25:   }
Chris@25: 
Chris@25:   // Applies a frequency-space filter to input, and returns the real-space
Chris@25:   // filtered input.
Chris@25:   // filterer accepts freq, i, n and modifies freq.real and freq.imag.
Chris@25:   ComplexArray.prototype.frequencyMap = function(filterer) {
Chris@25:     return this.FFT().map(filterer).InvFFT()
Chris@25:   }
Chris@25: 
Chris@25:   exports.frequencyMap = function(input, filterer) {
Chris@25:     return ensureComplexArray(input).frequencyMap(filterer)
Chris@25:   }
Chris@25: 
Chris@25:   function ensureComplexArray(input) {
Chris@25:     return complex_array.isComplexArray(input) && input ||
Chris@25:         new ComplexArray(input)
Chris@25:   }
Chris@25: 
Chris@25:   function FFT(input, inverse) {
Chris@25:     var n = input.length
Chris@25: 
Chris@25:     if (n & (n - 1)) {
Chris@25:       return FFT_Recursive(input, inverse)
Chris@25:     } else {
Chris@25:       return FFT_2_Iterative(input, inverse)
Chris@25:     }
Chris@25:   }
Chris@25: 
Chris@25:   function FFT_Recursive(input, inverse) {
Chris@25:     var
Chris@25:       n = input.length,
Chris@25:       // Counters.
Chris@25:       i, j,
Chris@25:       output,
Chris@25:       // Complex multiplier and its delta.
Chris@25:       f_r, f_i, del_f_r, del_f_i,
Chris@25:       // Lowest divisor and remainder.
Chris@25:       p, m,
Chris@25:       normalisation,
Chris@25:       recursive_result,
Chris@25:       _swap, _real, _imag
Chris@25: 
Chris@25:     if (n === 1) {
Chris@25:       return input
Chris@25:     }
Chris@25: 
Chris@25:     output = new ComplexArray(n, input.ArrayType)
Chris@25: 
Chris@25:     // Use the lowest odd factor, so we are able to use FFT_2_Iterative in the
Chris@25:     // recursive transforms optimally.
Chris@25:     p = LowestOddFactor(n)
Chris@25:     m = n / p
Chris@25:     normalisation = 1 / sqrt(p)
Chris@25:     recursive_result = new ComplexArray(m, input.ArrayType)
Chris@25: 
Chris@25:     // Loops go like O(n Σ p_i), where p_i are the prime factors of n.
Chris@25:     // for a power of a prime, p, this reduces to O(n p log_p n)
Chris@25:     for(j = 0; j < p; j++) {
Chris@25:       for(i = 0; i < m; i++) {
Chris@25:         recursive_result.real[i] = input.real[i * p + j]
Chris@25:         recursive_result.imag[i] = input.imag[i * p + j]
Chris@25:       }
Chris@25:       // Don't go deeper unless necessary to save allocs.
Chris@25:       if (m > 1) {
Chris@25:         recursive_result = FFT(recursive_result, inverse)
Chris@25:       }
Chris@25: 
Chris@25:       del_f_r = cos(2*PI*j/n)
Chris@25:       del_f_i = (inverse ? -1 : 1) * sin(2*PI*j/n)
Chris@25:       f_r = 1
Chris@25:       f_i = 0
Chris@25: 
Chris@25:       for(i = 0; i < n; i++) {
Chris@25:         _real = recursive_result.real[i % m]
Chris@25:         _imag = recursive_result.imag[i % m]
Chris@25: 
Chris@25:         output.real[i] += f_r * _real - f_i * _imag
Chris@25:         output.imag[i] += f_r * _imag + f_i * _real
Chris@25: 
Chris@25:         _swap = f_r * del_f_r - f_i * del_f_i
Chris@25:         f_i = f_r * del_f_i + f_i * del_f_r
Chris@25:         f_r = _swap
Chris@25:       }
Chris@25:     }
Chris@25: 
Chris@25:     // Copy back to input to match FFT_2_Iterative in-placeness
Chris@25:     // TODO: faster way of making this in-place?
Chris@25:     for(i = 0; i < n; i++) {
Chris@25:       input.real[i] = normalisation * output.real[i]
Chris@25:       input.imag[i] = normalisation * output.imag[i]
Chris@25:     }
Chris@25: 
Chris@25:     return input
Chris@25:   }
Chris@25: 
Chris@25:   function FFT_2_Iterative(input, inverse) {
Chris@25:     var
Chris@25:       n = input.length,
Chris@25:       // Counters.
Chris@25:       i, j,
Chris@25:       output, output_r, output_i,
Chris@25:       // Complex multiplier and its delta.
Chris@25:       f_r, f_i, del_f_r, del_f_i, temp,
Chris@25:       // Temporary loop variables.
Chris@25:       l_index, r_index,
Chris@25:       left_r, left_i, right_r, right_i,
Chris@25:       // width of each sub-array for which we're iteratively calculating FFT.
Chris@25:       width
Chris@25: 
Chris@25:     output = BitReverseComplexArray(input)
Chris@25:     output_r = output.real
Chris@25:     output_i = output.imag
Chris@25:     // Loops go like O(n log n):
Chris@25:     //   width ~ log n; i,j ~ n
Chris@25:     width = 1
Chris@25:     while (width < n) {
Chris@25:       del_f_r = cos(PI/width)
Chris@25:       del_f_i = (inverse ? -1 : 1) * sin(PI/width)
Chris@25:       for (i = 0; i < n/(2*width); i++) {
Chris@25:         f_r = 1
Chris@25:         f_i = 0
Chris@25:         for (j = 0; j < width; j++) {
Chris@25:           l_index = 2*i*width + j
Chris@25:           r_index = l_index + width
Chris@25: 
Chris@25:           left_r = output_r[l_index]
Chris@25:           left_i = output_i[l_index]
Chris@25:           right_r = f_r * output_r[r_index] - f_i * output_i[r_index]
Chris@25:           right_i = f_i * output_r[r_index] + f_r * output_i[r_index]
Chris@25: 
Chris@25:           output_r[l_index] = SQRT1_2 * (left_r + right_r)
Chris@25:           output_i[l_index] = SQRT1_2 * (left_i + right_i)
Chris@25:           output_r[r_index] = SQRT1_2 * (left_r - right_r)
Chris@25:           output_i[r_index] = SQRT1_2 * (left_i - right_i)
Chris@25:           temp = f_r * del_f_r - f_i * del_f_i
Chris@25:           f_i = f_r * del_f_i + f_i * del_f_r
Chris@25:           f_r = temp
Chris@25:         }
Chris@25:       }
Chris@25:       width <<= 1
Chris@25:     }
Chris@25: 
Chris@25:     return output
Chris@25:   }
Chris@25: 
Chris@25:   function BitReverseIndex(index, n) {
Chris@25:     var bitreversed_index = 0
Chris@25: 
Chris@25:     while (n > 1) {
Chris@25:       bitreversed_index <<= 1
Chris@25:       bitreversed_index += index & 1
Chris@25:       index >>= 1
Chris@25:       n >>= 1
Chris@25:     }
Chris@25:     return bitreversed_index
Chris@25:   }
Chris@25: 
Chris@25:   function BitReverseComplexArray(array) {
Chris@25:     var n = array.length,
Chris@25:         flips = {},
Chris@25:         swap,
Chris@25:         i
Chris@25: 
Chris@25:     for(i = 0; i < n; i++) {
Chris@25:       var r_i = BitReverseIndex(i, n)
Chris@25: 
Chris@25:       if (flips.hasOwnProperty(i) || flips.hasOwnProperty(r_i)) continue
Chris@25: 
Chris@25:       swap = array.real[r_i]
Chris@25:       array.real[r_i] = array.real[i]
Chris@25:       array.real[i] = swap
Chris@25: 
Chris@25:       swap = array.imag[r_i]
Chris@25:       array.imag[r_i] = array.imag[i]
Chris@25:       array.imag[i] = swap
Chris@25: 
Chris@25:       flips[i] = flips[r_i] = true
Chris@25:     }
Chris@25: 
Chris@25:     return array
Chris@25:   }
Chris@25: 
Chris@25:   function LowestOddFactor(n) {
Chris@25:     var factor = 3,
Chris@25:         sqrt_n = sqrt(n)
Chris@25: 
Chris@25:     while(factor <= sqrt_n) {
Chris@25:       if (n % factor === 0) return factor
Chris@25:       factor = factor + 2
Chris@25:     }
Chris@25:     return n
Chris@25:   }
Chris@25: 
Chris@25: }(
Chris@25:   typeof exports === 'undefined' && (this.fft = {}) || exports,
Chris@25:   typeof require === 'undefined' && (this.complex_array) ||
Chris@25:     require('./complex_array')
Chris@25: )