// On-device landmark (Shazam/Dejavu-style) audio fingerprinter — pure JS.
//
// A faithful port of the app's Swift `LandmarkFingerprint` (Sources/LicenseCore/
// LandmarkFingerprint.swift) and its Python twin (EjePlayerCodes/extension/
// fp_ref.py). It MUST stay parameter-identical to them so the hashes it emits
// collide with the bundled index (data/fp/), which was built by the Swift code.
// Validated bit-for-bit against fp_ref.py on identical PCM (tests/fp.engine.test.mjs).
//
// Pipeline: mono 8 kHz Float32 → STFT(1024/256, Hann) → log(1+|X|^2) spectrogram
// (512 bins, bin 0 packs DC^2+Nyquist^2 like vDSP) → separable sliding-max local
// peaks → 16-band adaptive threshold (mean+0.5σ) → density cap (6 peaks / 24
// frames) → pair hashes (fanOut 6, dt 1..60), code = (f1&0x1FF)<<16 |
// (f2&0x1FF)<<7 | (dt&0x7F), anchored at the first peak's frame.
//
// Decoding/resampling to 8 kHz mono is the caller's job (Web Audio in the
// browser; ffmpeg in tests) — this module is DOM-free and takes raw samples.

export const FP = {
  sr: 8000, nfft: 1024, hop: 256, bins: 512,
  peakFreqRadius: 10, peakTimeRadius: 10, peakThreshSigma: 0.5,
  freqBands: 16, densityWindow: 24, peaksPerWindow: 6,
  fanOut: 6, dtMin: 1, dtMax: 60,
};

// ---- radix-2 FFT (double precision, matches numpy's fft to ~1e-10) ----------
// Precompute bit-reversal + twiddles for the fixed nfft once.
const N = FP.nfft;
const _rev = new Uint16Array(N);
for (let i = 0; i < N; i++) {
  let r = 0;
  for (let b = 0, x = i; b < 10; b++, x >>= 1) r = (r << 1) | (x & 1);
  _rev[i] = r;
}
const _cos = new Float64Array(N / 2);
const _sin = new Float64Array(N / 2);
for (let i = 0; i < N / 2; i++) {
  _cos[i] = Math.cos(-2 * Math.PI * i / N);
  _sin[i] = Math.sin(-2 * Math.PI * i / N);
}
const _re = new Float64Array(N);
const _im = new Float64Array(N);

// In-place forward FFT of a real signal already loaded (bit-reversed) into _re/_im.
function fftReal(frame) {
  for (let i = 0; i < N; i++) { _re[i] = frame[_rev[i]]; _im[i] = 0; }
  for (let len = 2; len <= N; len <<= 1) {
    const half = len >> 1;
    const step = N / len;
    for (let i = 0; i < N; i += len) {
      for (let j = 0, k = 0; j < half; j++, k += step) {
        const c = _cos[k], s = _sin[k];
        const a = i + j, b = a + half;
        const tr = _re[b] * c - _im[b] * s;
        const ti = _re[b] * s + _im[b] * c;
        _re[b] = _re[a] - tr; _im[b] = _im[a] - ti;
        _re[a] += tr;        _im[a] += ti;
      }
    }
  }
}

// ---- Hann window (symmetric, == np.hanning == vDSP Hann shape) ---------------
const _win = new Float64Array(N);
for (let n = 0; n < N; n++) _win[n] = 0.5 - 0.5 * Math.cos((2 * Math.PI * n) / (N - 1));

// ---- spectrogram: (frames × bins) of log(1+|X|^2) ---------------------------
function spectrogram(x) {
  const { hop, bins } = FP;
  if (x.length < N) return { mags: new Float64Array(0), frames: 0 };
  const frames = 1 + ((x.length - N) / hop | 0);
  if (frames < 2) return { mags: new Float64Array(0), frames: 0 };
  const mags = new Float64Array(frames * bins);
  const frame = new Float64Array(N);
  for (let f = 0; f < frames; f++) {
    const start = f * hop;
    for (let n = 0; n < N; n++) frame[n] = x[start + n] * _win[n];
    fftReal(frame);
    const base = f * bins;
    // bin 0 packs DC^2 + Nyquist^2 (vDSP real-FFT layout)
    mags[base] = Math.log1p(_re[0] * _re[0] + _re[bins] * _re[bins]);
    for (let b = 1; b < bins; b++) {
      mags[base + b] = Math.log1p(_re[b] * _re[b] + _im[b] * _im[b]);
    }
  }
  return { mags, frames };
}

// ---- clamped rectangular sliding maximum ------------------------------------
// alongFreq=true → window over bins (columns); false → over frames (rows).
function slidingMax(a, frames, bins, radius, alongFreq) {
  const out = new Float64Array(a.length);
  if (alongFreq) {
    for (let r = 0; r < frames; r++) {
      const base = r * bins;
      for (let c = 0; c < bins; c++) {
        let m = -Infinity;
        const lo = Math.max(0, c - radius), hi = Math.min(bins - 1, c + radius);
        for (let k = lo; k <= hi; k++) { const v = a[base + k]; if (v > m) m = v; }
        out[base + c] = m;
      }
    }
  } else {
    for (let c = 0; c < bins; c++) {
      for (let r = 0; r < frames; r++) {
        let m = -Infinity;
        const lo = Math.max(0, r - radius), hi = Math.min(frames - 1, r + radius);
        for (let k = lo; k <= hi; k++) { const v = a[k * bins + c]; if (v > m) m = v; }
        out[r * bins + c] = m;
      }
    }
  }
  return out;
}

// ---- peak constellation -----------------------------------------------------
function constellation(x) {
  const { bins, peakFreqRadius, peakTimeRadius, peakThreshSigma, freqBands,
          densityWindow, peaksPerWindow } = FP;
  const { mags, frames } = spectrogram(x);
  if (frames < 2) return [];

  const fMax = slidingMax(mags, frames, bins, peakFreqRadius, true);
  const tfMax = slidingMax(fMax, frames, bins, peakTimeRadius, false);

  // adaptive per-band thresholds: mean + σ·std over each band's cells
  const bandSize = Math.max(1, (bins / freqBands) | 0);   // 32
  const bandThresh = new Float64Array(freqBands);
  for (let b = 0; b < freqBands; b++) {
    const lo = b * bandSize, hi = Math.min(bins, lo + bandSize);
    if (lo >= hi) continue;
    let sum = 0, sumsq = 0, cnt = 0;
    for (let t = 0; t < frames; t++) {
      const rb = t * bins;
      for (let f = lo; f < hi; f++) { const v = mags[rb + f]; sum += v; sumsq += v * v; cnt++; }
    }
    const mean = sum / cnt;
    bandThresh[b] = mean + peakThreshSigma * Math.sqrt(Math.max(0, sumsq / cnt - mean * mean));
  }

  // candidate peaks in (t asc, f asc) order — matches fp_ref np.where row-major
  const cand = [];   // {t, f, mag}
  for (let t = 0; t < frames; t++) {
    const base = t * bins;
    for (let f = 0; f < bins; f++) {
      const v = mags[base + f];
      if (v === tfMax[base + f] && v > bandThresh[Math.min(freqBands - 1, (f / bandSize) | 0)]) {
        cand.push({ t, f, mag: v });
      }
    }
  }
  if (!cand.length) return [];

  // density cap: per densityWindow keep the loudest peaksPerWindow (stable)
  const out = [];
  let i = 0;
  while (i < cand.length) {
    const wStart = cand[i].t - (cand[i].t % densityWindow);
    let j = i;
    const window = [];
    while (j < cand.length && cand[j].t < wStart + densityWindow) { window.push(cand[j]); j++; }
    window.sort((a, b) => b.mag - a.mag);   // stable in modern JS engines
    for (const p of window.slice(0, peaksPerWindow)) out.push({ t: p.t, f: p.f });
    i = j;
  }
  return out;
}

// Linear-resample `samples` by `ratio` (output length = len*ratio). ratio>1
// stretches/slows (lowers pitch); ratio<1 compresses/speeds (raises pitch).
// Used to fingerprint a query at several speeds: old tango transfers (and their
// YouTube re-uploads) are commonly digitized a few % fast/slow vs the library
// copy, which shifts every spectral peak into different bins AND scales the
// inter-peak Δt — so the only way the landmark hashes collide again is to undo
// the speed. Resampling the query by the right ratio restores both. (Validated
// against real YouTube audio: a 5.5%-slow upload recovers at ratio ≈ 1.055.)
export function resample(samples, ratio) {
  if (ratio === 1) return samples;
  const n = Math.max(1, Math.round(samples.length * ratio));
  const out = new Float32Array(n);
  const inv = 1 / ratio;
  const last = samples.length - 1;
  for (let i = 0; i < n; i++) {
    const x = i * inv;
    const x0 = x | 0;
    if (x0 >= last) { out[i] = samples[last]; break; }
    const f = x - x0;
    out[i] = samples[x0] * (1 - f) + samples[x0 + 1] * f;
  }
  return out;
}

// ---- pair hashing -----------------------------------------------------------
// Returns parallel arrays { codes: Uint32Array, ts: Uint32Array }.
export function fingerprint(samples) {
  const peaks = constellation(samples);
  const { fanOut, dtMin, dtMax } = FP;
  if (peaks.length < 2) return { codes: new Uint32Array(0), ts: new Uint32Array(0) };
  peaks.sort((a, b) => a.t - b.t);
  const codes = [], ts = [];
  const n = peaks.length;
  for (let i = 0; i < n; i++) {
    let j = i + 1, cnt = 0;
    while (j < n && cnt < fanOut) {
      const dt = peaks[j].t - peaks[i].t;
      if (dt > dtMax) break;
      if (dt >= dtMin) {
        const f1 = peaks[i].f & 0x1FF, f2 = peaks[j].f & 0x1FF;
        codes.push(((f1 << 16) | (f2 << 7) | (dt & 0x7F)) >>> 0);
        ts.push(peaks[i].t);
        cnt++;
      }
      j++;
    }
  }
  return { codes: Uint32Array.from(codes), ts: Uint32Array.from(ts) };
}
