/**
 * Digital Signal Processing utilities
 * FFT, windowing functions, and spectral analysis
 */

/**
 * Complex number representation
 */
export interface Complex {
  real: number;
  imag: number;
}

/**
 * Window function types
 */
export type WindowType = 'hann' | 'hamming' | 'blackman' | 'rectangular' | 'bartlett';

/**
 * FFT result with magnitude and phase
 */
export interface FFTResult {
  magnitude: Float32Array;
  phase: Float32Array;
  frequencies: Float32Array;
}

/**
 * Fast Fourier Transform (Cooley-Tukey radix-2 decimation-in-time)
 */
export class FFT {
  private size: number;
  private reversedBits: Uint32Array;
  private cosTable: Float32Array;
  private sinTable: Float32Array;

  constructor(size: number) {
    if (!this.isPowerOfTwo(size)) {
      throw new Error('FFT size must be a power of 2');
    }

    this.size = size;
    this.reversedBits = new Uint32Array(size);
    this.cosTable = new Float32Array(size / 2);
    this.sinTable = new Float32Array(size / 2);

    // Precompute bit-reversed indices
    const bits = Math.log2(size);
    for (let i = 0; i < size; i++) {
      this.reversedBits[i] = this.reverseBits(i, bits);
    }

    // Precompute twiddle factors
    for (let i = 0; i < size / 2; i++) {
      const angle = -2 * Math.PI * i / size;
      this.cosTable[i] = Math.cos(angle);
      this.sinTable[i] = Math.sin(angle);
    }
  }

  /**
   * Forward FFT (time domain to frequency domain)
   */
  forward(input: Float32Array): Complex[] {
    if (input.length !== this.size) {
      throw new Error(`Input size must be ${this.size}`);
    }

    const output: Complex[] = new Array(this.size);

    // Bit-reverse copy
    for (let i = 0; i < this.size; i++) {
      output[this.reversedBits[i]] = {
        real: input[i],
        imag: 0,
      };
    }

    // Cooley-Tukey FFT
    for (let len = 2; len <= this.size; len *= 2) {
      const halfLen = len / 2;
      const tableStep = this.size / len;

      for (let i = 0; i < this.size; i += len) {
        for (let j = 0; j < halfLen; j++) {
          const k = j * tableStep;
          const cos = this.cosTable[k];
          const sin = this.sinTable[k];

          const even = output[i + j];
          const odd = output[i + j + halfLen];

          const tReal = odd.real * cos - odd.imag * sin;
          const tImag = odd.real * sin + odd.imag * cos;

          output[i + j] = {
            real: even.real + tReal,
            imag: even.imag + tImag,
          };

          output[i + j + halfLen] = {
            real: even.real - tReal,
            imag: even.imag - tImag,
          };
        }
      }
    }

    return output;
  }

  /**
   * Get magnitude spectrum
   */
  getMagnitude(complex: Complex[]): Float32Array {
    const magnitude = new Float32Array(complex.length);
    for (let i = 0; i < complex.length; i++) {
      magnitude[i] = Math.sqrt(
        complex[i].real * complex[i].real + 
        complex[i].imag * complex[i].imag
      );
    }
    return magnitude;
  }

  /**
   * Get phase spectrum
   */
  getPhase(complex: Complex[]): Float32Array {
    const phase = new Float32Array(complex.length);
    for (let i = 0; i < complex.length; i++) {
      phase[i] = Math.atan2(complex[i].imag, complex[i].real);
    }
    return phase;
  }

  /**
   * Get power spectrum (magnitude squared)
   */
  getPowerSpectrum(complex: Complex[]): Float32Array {
    const power = new Float32Array(complex.length);
    for (let i = 0; i < complex.length; i++) {
      power[i] = 
        complex[i].real * complex[i].real + 
        complex[i].imag * complex[i].imag;
    }
    return power;
  }

  /**
   * Compute FFT and return magnitude, phase, and frequencies
   */
  compute(input: Float32Array, sampleRate: number): FFTResult {
    const complex = this.forward(input);
    const magnitude = this.getMagnitude(complex);
    const phase = this.getPhase(complex);
    const frequencies = this.getFrequencies(sampleRate);

    return { magnitude, phase, frequencies };
  }

  /**
   * Get frequency bins
   */
  getFrequencies(sampleRate: number): Float32Array {
    const frequencies = new Float32Array(this.size);
    const binWidth = sampleRate / this.size;
    
    for (let i = 0; i < this.size; i++) {
      frequencies[i] = i * binWidth;
    }
    
    return frequencies;
  }

  private isPowerOfTwo(n: number): boolean {
    return n > 0 && (n & (n - 1)) === 0;
  }

  private reverseBits(x: number, bits: number): number {
    let result = 0;
    for (let i = 0; i < bits; i++) {
      result = (result << 1) | (x & 1);
      x >>= 1;
    }
    return result;
  }
}

/**
 * Window functions for reducing spectral leakage
 */
export class WindowFunction {
  /**
   * Create a window function of specified type and size
   */
  static create(size: number, type: WindowType = 'hann'): Float32Array {
    switch (type) {
      case 'hann':
        return WindowFunction.hann(size);
      case 'hamming':
        return WindowFunction.hamming(size);
      case 'blackman':
        return WindowFunction.blackman(size);
      case 'bartlett':
        return WindowFunction.bartlett(size);
      case 'rectangular':
        return WindowFunction.rectangular(size);
      default:
        throw new Error(`Unknown window type: ${type}`);
    }
  }

  /**
   * Hann window (raised cosine)
   */
  static hann(size: number): Float32Array {
    const window = new Float32Array(size);
    for (let i = 0; i < size; i++) {
      window[i] = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (size - 1)));
    }
    return window;
  }

  /**
   * Hamming window
   */
  static hamming(size: number): Float32Array {
    const window = new Float32Array(size);
    const alpha = 0.54;
    const beta = 1 - alpha;
    
    for (let i = 0; i < size; i++) {
      window[i] = alpha - beta * Math.cos((2 * Math.PI * i) / (size - 1));
    }
    return window;
  }

  /**
   * Blackman window
   */
  static blackman(size: number): Float32Array {
    const window = new Float32Array(size);
    const a0 = 0.42;
    const a1 = 0.5;
    const a2 = 0.08;
    
    for (let i = 0; i < size; i++) {
      const n = i / (size - 1);
      window[i] = a0 - 
                   a1 * Math.cos(2 * Math.PI * n) + 
                   a2 * Math.cos(4 * Math.PI * n);
    }
    return window;
  }

  /**
   * Bartlett window (triangular)
   */
  static bartlett(size: number): Float32Array {
    const window = new Float32Array(size);
    const half = (size - 1) / 2;
    
    for (let i = 0; i < size; i++) {
      window[i] = 1 - Math.abs((i - half) / half);
    }
    return window;
  }

  /**
   * Rectangular window (no windowing)
   */
  static rectangular(size: number): Float32Array {
    const window = new Float32Array(size);
    window.fill(1);
    return window;
  }

  /**
   * Apply window to signal
   */
  static apply(signal: Float32Array, window: Float32Array): Float32Array {
    if (signal.length !== window.length) {
      throw new Error('Signal and window must have same length');
    }

    const windowed = new Float32Array(signal.length);
    for (let i = 0; i < signal.length; i++) {
      windowed[i] = signal[i] * window[i];
    }
    return windowed;
  }
}

/**
 * Spectral analysis utilities
 */
export class SpectralAnalysis {
  /**
   * Find spectral peaks
   */
  static findPeaks(
    spectrum: Float32Array,
    threshold: number = 0.1,
    minDistance: number = 1
  ): number[] {
    const peaks: number[] = [];
    const maxValue = Math.max(...spectrum);
    const thresholdValue = maxValue * threshold;

    for (let i = minDistance; i < spectrum.length - minDistance; i++) {
      // Check if this is a local maximum
      let isPeak = spectrum[i] >= thresholdValue;
      
      for (let j = i - minDistance; j <= i + minDistance && isPeak; j++) {
        if (j !== i && spectrum[j] >= spectrum[i]) {
          isPeak = false;
        }
      }

      if (isPeak) {
        peaks.push(i);
      }
    }

    return peaks;
  }

  /**
   * Parabolic interpolation for more accurate peak frequency
   */
  static parabolicInterpolation(
    spectrum: Float32Array,
    peakIndex: number
  ): number {
    if (peakIndex <= 0 || peakIndex >= spectrum.length - 1) {
      return peakIndex;
    }

    const alpha = spectrum[peakIndex - 1];
    const beta = spectrum[peakIndex];
    const gamma = spectrum[peakIndex + 1];

    const p = 0.5 * (alpha - gamma) / (alpha - 2 * beta + gamma);
    
    return peakIndex + p;
  }

  /**
   * Convert bin index to frequency
   */
  static binToFrequency(bin: number, sampleRate: number, fftSize: number): number {
    return (bin * sampleRate) / fftSize;
  }

  /**
   * Convert frequency to bin index
   */
  static frequencyToBin(frequency: number, sampleRate: number, fftSize: number): number {
    return Math.round((frequency * fftSize) / sampleRate);
  }

  /**
   * Compute spectral centroid (center of mass of spectrum)
   */
  static spectralCentroid(magnitude: Float32Array, frequencies: Float32Array): number {
    let numerator = 0;
    let denominator = 0;

    for (let i = 0; i < magnitude.length; i++) {
      numerator += magnitude[i] * frequencies[i];
      denominator += magnitude[i];
    }

    return denominator > 0 ? numerator / denominator : 0;
  }

  /**
   * Compute spectral flux (change in spectrum over time)
   */
  static spectralFlux(
    currentSpectrum: Float32Array,
    previousSpectrum: Float32Array
  ): number {
    if (currentSpectrum.length !== previousSpectrum.length) {
      throw new Error('Spectra must have same length');
    }

    let flux = 0;
    for (let i = 0; i < currentSpectrum.length; i++) {
      const diff = currentSpectrum[i] - previousSpectrum[i];
      flux += diff * diff;
    }

    return Math.sqrt(flux);
  }

  /**
   * Zero-pad signal to next power of 2
   */
  static zeroPad(signal: Float32Array, targetSize?: number): Float32Array {
    if (!targetSize) {
      targetSize = 1;
      while (targetSize < signal.length) {
        targetSize *= 2;
      }
    }

    const padded = new Float32Array(targetSize);
    padded.set(signal);
    return padded;
  }

  /**
   * Compute Harmonic Product Spectrum (for pitch detection)
   */
  static harmonicProductSpectrum(
    spectrum: Float32Array,
    harmonics: number = 4
  ): Float32Array {
    const hps = new Float32Array(spectrum.length);
    hps.set(spectrum);

    // Multiply spectrum with its downsampled versions
    for (let h = 2; h <= harmonics; h++) {
      for (let i = 0; i < Math.floor(spectrum.length / h); i++) {
        hps[i] *= spectrum[i * h];
      }
    }

    return hps;
  }

  /**
   * Autocorrelation using FFT (faster for large signals)
   */
  static autocorrelation(signal: Float32Array): Float32Array {
    const fftSize = 1;
    let size = fftSize;
    while (size < signal.length * 2) {
      size *= 2;
    }

    const fft = new FFT(size);
    
    // Zero-pad signal
    const padded = new Float32Array(size);
    padded.set(signal);

    // FFT of signal
    const complex = fft.forward(padded);
    
    // Multiply by complex conjugate (power spectrum)
    for (let i = 0; i < complex.length; i++) {
      const mag = complex[i].real * complex[i].real + complex[i].imag * complex[i].imag;
      complex[i].real = mag;
      complex[i].imag = 0;
    }

    // IFFT to get autocorrelation (not implemented here, use direct method for now)
    // For now, return direct autocorrelation
    return this.directAutocorrelation(signal);
  }

  /**
   * Direct autocorrelation (simpler but slower)
   */
  private static directAutocorrelation(signal: Float32Array): Float32Array {
    const result = new Float32Array(signal.length);
    
    for (let lag = 0; lag < signal.length; lag++) {
      let sum = 0;
      for (let i = 0; i < signal.length - lag; i++) {
        sum += signal[i] * signal[i + lag];
      }
      result[lag] = sum;
    }

    return result;
  }
}

/**
 * Helper: Get next power of 2
 */
export function nextPowerOfTwo(n: number): number {
  return Math.pow(2, Math.ceil(Math.log2(n)));
}

/**
 * Helper: Check if number is power of 2
 */
export function isPowerOfTwo(n: number): boolean {
  return n > 0 && (n & (n - 1)) === 0;
}