/**
 * Harmonic Analysis Module
 * Implements Harmonic Product Spectrum (HPS) and harmonic matching algorithms
 */

import { FFT, WindowFunction, SpectralAnalysis } from './dsp';

/**
 * Harmonic analysis configuration
 */
export interface HarmonicAnalysisConfig {
  sampleRate?: number;
  fftSize?: number;
  windowType?: 'hann' | 'hamming' | 'blackman' | 'rectangular';
  maxHarmonics?: number;
  harmonicTolerance?: number;
  minHarmonicStrength?: number;
  useHPS?: boolean;
  hpsDecimation?: number;
}

/**
 * Harmonic analysis result
 */
export interface HarmonicResult {
  fundamental: number;
  confidence: number;
  harmonics: Array<{
    frequency: number;
    amplitude: number;
    harmonicNumber: number;
    strength: number;
  }>;
  harmonicity: number;
  spectralCentroid: number;
  spectralRolloff: number;
}

/**
 * Harmonic Product Spectrum (HPS) analyzer
 */
export class HarmonicAnalyzer {
  private config: Required<HarmonicAnalysisConfig>;
  private fft: FFT;
  private window: Float32Array;

  constructor(config: HarmonicAnalysisConfig = {}) {
    this.config = {
      sampleRate: config.sampleRate ?? 44100,
      fftSize: config.fftSize ?? 2048,
      windowType: config.windowType ?? 'hann',
      maxHarmonics: config.maxHarmonics ?? 8,
      harmonicTolerance: config.harmonicTolerance ?? 0.05,
      minHarmonicStrength: config.minHarmonicStrength ?? 0.1,
      useHPS: config.useHPS ?? true,
      hpsDecimation: config.hpsDecimation ?? 4,
    };

    this.fft = new FFT(this.config.fftSize);
    this.window = WindowFunction.create(this.config.fftSize, this.config.windowType);
  }

  /**
   * Analyze harmonics in a signal
   */
  analyzeHarmonics(signal: Float32Array): HarmonicResult | null {
    if (signal.length !== this.config.fftSize) {
      throw new Error(`Signal length must be ${this.config.fftSize}`);
    }

    // Apply window
    const windowedSignal = WindowFunction.apply(signal, this.window);

    // Compute FFT
    const fftResult = this.fft.forward(windowedSignal);
    const magnitude = new Float32Array(this.config.fftSize / 2 + 1);
    
    // Calculate magnitude from complex FFT result
    for (let i = 0; i < magnitude.length; i++) {
      const real = fftResult[i].real;
      const imag = fftResult[i].imag;
      magnitude[i] = Math.sqrt(real * real + imag * imag);
    }

    // Find fundamental using HPS or spectral analysis
    const fundamental = this.config.useHPS 
      ? this.findFundamentalHPS(magnitude)
      : this.findFundamentalSpectral(magnitude);

    if (!fundamental) {
      return null;
    }

    // Find harmonics
    const harmonics = this.findHarmonics(magnitude, fundamental.frequency);

    // Calculate harmonicity
    const harmonicity = this.calculateHarmonicity(harmonics);

    // Calculate spectral features
    const spectralCentroid = this.calculateSpectralCentroid(magnitude);
    const spectralRolloff = this.calculateSpectralRolloff(magnitude);

    return {
      fundamental: fundamental.frequency,
      confidence: fundamental.confidence,
      harmonics,
      harmonicity,
      spectralCentroid,
      spectralRolloff,
    };
  }

  /**
   * Find fundamental frequency using Harmonic Product Spectrum
   */
  private findFundamentalHPS(magnitude: Float32Array): { frequency: number; confidence: number } | null {
    const hps = this.computeHPS(magnitude);
    
    // Find peak in HPS
    const peak = this.findPeakInHPS(hps);
    if (!peak) {
      return null;
    }

    const frequency = (peak.bin * this.config.sampleRate) / this.config.fftSize;
    const confidence = Math.min(1.0, peak.strength);

    return { frequency, confidence };
  }

  /**
   * Compute Harmonic Product Spectrum
   */
  private computeHPS(magnitude: Float32Array): Float32Array {
    const hps = new Float32Array(magnitude.length);
    const decimation = this.config.hpsDecimation;

    for (let i = 0; i < magnitude.length; i++) {
      let product = magnitude[i];
      
      // Multiply by downsampled versions
      for (let harmonic = 2; harmonic <= decimation; harmonic++) {
        const downsampledIndex = Math.floor(i / harmonic);
        if (downsampledIndex < magnitude.length && downsampledIndex >= 0) {
          product *= magnitude[downsampledIndex];
        }
      }
      
      hps[i] = product;
    }

    // Normalize HPS
    const maxValue = Math.max(...hps);
    if (maxValue > 0) {
      for (let i = 0; i < hps.length; i++) {
        hps[i] /= maxValue;
      }
    }

    return hps;
  }

  /**
   * Find peak in HPS
   */
  private findPeakInHPS(hps: Float32Array): { bin: number; strength: number } | null {
    const minBin = Math.floor((50 * this.config.fftSize) / this.config.sampleRate);
    const maxBin = Math.floor((2000 * this.config.fftSize) / this.config.sampleRate);
    
    let maxStrength = 0;
    let bestBin = -1;

    for (let bin = minBin; bin <= maxBin && bin < hps.length; bin++) {
      if (hps[bin] > maxStrength) {
        maxStrength = hps[bin];
        bestBin = bin;
      }
    }

    if (bestBin === -1 || maxStrength < this.config.minHarmonicStrength) {
      return null;
    }

    return { bin: bestBin, strength: maxStrength };
  }

  /**
   * Find fundamental using spectral analysis
   */
  private findFundamentalSpectral(magnitude: Float32Array): { frequency: number; confidence: number } | null {
    // Find spectral peaks
    const peaks = this.findSpectralPeaks(magnitude);
    if (peaks.length === 0) {
      return null;
    }

    // Find fundamental by looking for lowest frequency with strong harmonics
    const fundamental = this.findFundamentalFromPeaks(peaks);
    if (!fundamental) {
      return null;
    }

    return {
      frequency: fundamental.frequency,
      confidence: fundamental.strength,
    };
  }

  /**
   * Find spectral peaks
   */
  private findSpectralPeaks(magnitude: Float32Array): Array<{ frequency: number; strength: number }> {
    const peaks: Array<{ frequency: number; strength: number }> = [];
    const minBin = Math.floor((50 * this.config.fftSize) / this.config.sampleRate);
    const maxBin = Math.floor((2000 * this.config.fftSize) / this.config.sampleRate);

    for (let bin = minBin; bin <= maxBin && bin < magnitude.length - 1; bin++) {
      if (this.isLocalMaximum(magnitude, bin)) {
        const frequency = (bin * this.config.sampleRate) / this.config.fftSize;
        const strength = magnitude[bin];
        
        if (strength > this.config.minHarmonicStrength) {
          peaks.push({ frequency, strength });
        }
      }
    }

    return peaks.sort((a, b) => b.strength - a.strength);
  }

  /**
   * Check if bin is a local maximum
   */
  private isLocalMaximum(magnitude: Float32Array, bin: number): boolean {
    if (bin <= 0 || bin >= magnitude.length - 1) {
      return false;
    }

    const center = magnitude[bin];
    const left = magnitude[bin - 1];
    const right = magnitude[bin + 1];

    return center > left && center > right;
  }

  /**
   * Find fundamental from spectral peaks
   */
  private findFundamentalFromPeaks(peaks: Array<{ frequency: number; strength: number }>): { frequency: number; strength: number } | null {
    if (peaks.length === 0) {
      return null;
    }

    // Try to find fundamental by looking for harmonic relationships
    for (const peak of peaks) {
      const harmonicCount = this.countHarmonics(peaks, peak.frequency);
      
      if (harmonicCount >= 2) {
        return peak;
      }
    }

    // Fall back to strongest peak
    return peaks[0];
  }

  /**
   * Count harmonics of a fundamental frequency
   */
  private countHarmonics(peaks: Array<{ frequency: number; strength: number }>, fundamental: number): number {
    let count = 0;
    const tolerance = this.config.harmonicTolerance;

    for (let harmonic = 2; harmonic <= this.config.maxHarmonics; harmonic++) {
      const expectedFreq = fundamental * harmonic;
      
      for (const peak of peaks) {
        if (Math.abs(peak.frequency - expectedFreq) / expectedFreq < tolerance) {
          count++;
          break;
        }
      }
    }

    return count;
  }

  /**
   * Find harmonics of a fundamental frequency
   */
  private findHarmonics(magnitude: Float32Array, fundamental: number): Array<{
    frequency: number;
    amplitude: number;
    harmonicNumber: number;
    strength: number;
  }> {
    const harmonics: Array<{
      frequency: number;
      amplitude: number;
      harmonicNumber: number;
      strength: number;
    }> = [];

    for (let harmonic = 1; harmonic <= this.config.maxHarmonics; harmonic++) {
      const expectedFreq = fundamental * harmonic;
      const bin = Math.round((expectedFreq * this.config.fftSize) / this.config.sampleRate);
      
      if (bin < magnitude.length) {
        const amplitude = magnitude[bin];
        const strength = amplitude / magnitude[0]; // Normalize by DC component
        
        if (strength > this.config.minHarmonicStrength) {
          harmonics.push({
            frequency: expectedFreq,
            amplitude,
            harmonicNumber: harmonic,
            strength,
          });
        }
      }
    }

    return harmonics;
  }

  /**
   * Calculate harmonicity (ratio of harmonic to inharmonic energy)
   */
  private calculateHarmonicity(harmonics: Array<{
    frequency: number;
    amplitude: number;
    harmonicNumber: number;
    strength: number;
  }>): number {
    if (harmonics.length === 0) {
      return 0;
    }

    const totalHarmonicEnergy = harmonics.reduce((sum, h) => sum + h.amplitude * h.amplitude, 0);
    const totalEnergy = harmonics.reduce((sum, h) => sum + h.amplitude, 0);
    
    return totalHarmonicEnergy / (totalEnergy * totalEnergy + 1e-8);
  }

  /**
   * Update configuration
   */
  updateConfig(config: Partial<HarmonicAnalysisConfig>): void {
    this.config = { ...this.config, ...config };

    // Recreate FFT and window if size changed
    if (config.fftSize && config.fftSize !== this.config.fftSize) {
      this.fft = new FFT(this.config.fftSize);
      this.window = WindowFunction.create(this.config.fftSize, this.config.windowType);
    } else if (config.windowType && config.windowType !== this.config.windowType) {
      this.window = WindowFunction.create(this.config.fftSize, this.config.windowType);
    }
  }

  /**
   * Calculate spectral centroid
   */
  private calculateSpectralCentroid(magnitude: Float32Array): number {
    let weightedSum = 0;
    let magnitudeSum = 0;

    for (let i = 0; i < magnitude.length; i++) {
      const frequency = (i * this.config.sampleRate) / (2 * magnitude.length);
      weightedSum += frequency * magnitude[i];
      magnitudeSum += magnitude[i];
    }

    return magnitudeSum > 0 ? weightedSum / magnitudeSum : 0;
  }

  /**
   * Calculate spectral rolloff
   */
  private calculateSpectralRolloff(magnitude: Float32Array): number {
    const totalEnergy = magnitude.reduce((sum, val) => sum + val * val, 0);
    const threshold = 0.85 * totalEnergy; // 85% of total energy
    
    let cumulativeEnergy = 0;
    for (let i = 0; i < magnitude.length; i++) {
      cumulativeEnergy += magnitude[i] * magnitude[i];
      if (cumulativeEnergy >= threshold) {
        return (i * this.config.sampleRate) / (2 * magnitude.length);
      }
    }

    return (magnitude.length * this.config.sampleRate) / (2 * magnitude.length);
  }

  /**
   * Get current configuration
   */
  getConfig(): Readonly<Required<HarmonicAnalysisConfig>> {
    return { ...this.config };
  }
}

/**
 * Harmonic matching utilities
 */
export const HarmonicMatching = {
  /**
   * Match harmonics to a fundamental frequency
   */
  matchHarmonics(
    peaks: Array<{ frequency: number; strength: number }>,
    fundamental: number,
    tolerance: number = 0.05
  ): Array<{ frequency: number; strength: number; harmonicNumber: number }> {
    const matched: Array<{ frequency: number; strength: number; harmonicNumber: number }> = [];

    for (let harmonic = 1; harmonic <= 8; harmonic++) {
      const expectedFreq = fundamental * harmonic;
      
      for (const peak of peaks) {
        if (Math.abs(peak.frequency - expectedFreq) / expectedFreq < tolerance) {
          matched.push({
            ...peak,
            harmonicNumber: harmonic,
          });
          break;
        }
      }
    }

    return matched;
  },

  /**
   * Find best fundamental from harmonic matches
   */
  findBestFundamental(
    peaks: Array<{ frequency: number; strength: number }>,
    candidates: number[]
  ): { frequency: number; score: number } | null {
    let bestScore = 0;
    let bestFundamental: { frequency: number; score: number } | null = null;

    for (const candidate of candidates) {
      const matches = this.matchHarmonics(peaks, candidate);
      const score = this.calculateHarmonicScore(matches);
      
      if (score > bestScore) {
        bestScore = score;
        bestFundamental = { frequency: candidate, score };
      }
    }

    return bestFundamental;
  },

  /**
   * Calculate harmonic score
   */
  calculateHarmonicScore(matches: Array<{ frequency: number; strength: number; harmonicNumber: number }>): number {
    if (matches.length === 0) {
      return 0;
    }

    // Weight by harmonic number (lower harmonics are more important)
    const weightedScore = matches.reduce((sum, match) => {
      const weight = 1 / match.harmonicNumber;
      return sum + match.strength * weight;
    }, 0);

    return weightedScore / matches.length;
  },

  /**
   * Detect octave errors
   */
  detectOctaveError(
    detected: number,
    expected: number,
    tolerance: number = 0.1
  ): { isOctaveError: boolean; correctedFrequency: number } {
    const ratio = detected / expected;
    const octaveRatio = Math.round(ratio);
    
    if (Math.abs(ratio - octaveRatio) < tolerance && octaveRatio >= 2) {
      return {
        isOctaveError: true,
        correctedFrequency: detected / octaveRatio,
      };
    }

    return {
      isOctaveError: false,
      correctedFrequency: detected,
    };
  },
};

/**
 * Quick utility: Analyze harmonics
 */
export function analyzeHarmonics(
  signal: Float32Array,
  sampleRate: number = 44100,
  options?: HarmonicAnalysisConfig
): HarmonicResult | null {
  const analyzer = new HarmonicAnalyzer({
    sampleRate,
    ...options,
  });

  return analyzer.analyzeHarmonics(signal);
}

/**
 * Harmonic analysis utilities
 */
export const HarmonicUtils = {
  /**
   * Calculate harmonic-to-noise ratio
   */
  calculateHNR(harmonics: Array<{ amplitude: number }>, noiseFloor: number): number {
    const harmonicEnergy = harmonics.reduce((sum, h) => sum + h.amplitude * h.amplitude, 0);
    const noiseEnergy = noiseFloor * noiseFloor;
    
    return 10 * Math.log10(harmonicEnergy / (noiseEnergy + 1e-8));
  },

  /**
   * Calculate spectral irregularity
   */
  calculateSpectralIrregularity(magnitude: Float32Array): number {
    let irregularity = 0;
    
    for (let i = 1; i < magnitude.length - 1; i++) {
      const diff = magnitude[i] - (magnitude[i - 1] + magnitude[i + 1]) / 2;
      irregularity += diff * diff;
    }
    
    return irregularity / magnitude.length;
  },

  /**
   * Calculate spectral flatness
   */
  calculateSpectralFlatness(magnitude: Float32Array): number {
    let geometricMean = 1;
    let arithmeticMean = 0;
    
    for (const mag of magnitude) {
      if (mag > 0) {
        geometricMean *= Math.pow(mag, 1 / magnitude.length);
        arithmeticMean += mag;
      }
    }
    
    arithmeticMean /= magnitude.length;
    
    return geometricMean / (arithmeticMean + 1e-8);
  },

  /**
   * Extract harmonic series
   */
  extractHarmonicSeries(
    fundamental: number,
    maxHarmonics: number = 8
  ): number[] {
    const series: number[] = [];
    
    for (let i = 1; i <= maxHarmonics; i++) {
      series.push(fundamental * i);
    }
    
    return series;
  },
};