// Anti-aliased decimation: Kaiser-windowed-sinc FIR lowpass at ~0.45 * fs_lo,
// evaluated only at the retained output samples (polyphase, so the cost is
// taps * T / R rather than taps * T). Replaces scipy.resample_poly(p, 1, R).

function besselI0(x: number): number {
  // Series expansion; converges quickly for the argument range used here.
  let sum = 1.0;
  let term = 1.0;
  const half = x / 2.0;
  for (let k = 1; k < 64; k++) {
    term *= (half / k) * (half / k);
    sum += term;
    if (term < 1e-16 * sum) break;
  }
  return sum;
}

export function designDecimator(R: number, beta = 8.0): Float64Array {
  // Cutoff at 0.45 of the output Nyquist-setting rate: fc/fs_hi = 0.45/R.
  const taps = 12 * R + 1; // odd length, linear phase
  const h = new Float64Array(taps);
  const M2 = (taps - 1) / 2;
  const fcNorm = 0.45 / R; // cycles per input sample
  const i0beta = besselI0(beta);
  let sum = 0;
  for (let i = 0; i < taps; i++) {
    const m = i - M2;
    const sinc = m === 0 ? 2 * fcNorm : Math.sin(2 * Math.PI * fcNorm * m) / (Math.PI * m);
    const r = m / M2;
    const w = besselI0(beta * Math.sqrt(Math.max(0, 1 - r * r))) / i0beta;
    h[i] = sinc * w;
    sum += h[i];
  }
  for (let i = 0; i < taps; i++) h[i] /= sum; // unity DC gain
  return h;
}

export function decimate(p: Float64Array, R: number): Float64Array {
  if (R === 1) return p.slice();
  const h = designDecimator(R);
  const taps = h.length;
  const M2 = (taps - 1) / 2;
  const T = p.length;
  const outLen = Math.floor(T / R);
  const y = new Float64Array(outLen);
  for (let m = 0; m < outLen; m++) {
    const centre = m * R;
    let acc = 0;
    // y[m] = sum_i h[i] * p[centre + M2 - i], zero outside the buffer
    const iLo = Math.max(0, centre + M2 - (T - 1));
    const iHi = Math.min(taps - 1, centre + M2);
    for (let i = iLo; i <= iHi; i++) acc += h[i] * p[centre + M2 - i];
    y[m] = acc;
  }
  return y;
}