util.js

// ── Math shorthand ──

export let { cos, sin, abs, exp, sqrt, PI, cosh, acosh, acos, pow, log10 } = Math
export let PI2 = 2 * PI

// ── Internal helpers (used by window functions) ──

/** Cosine-sum: w(i) = Σ (-1)^k aₖ cos(2πki/(N-1)). Used by blackman, nuttall, etc. */
export function cosineSum (i, N, a) {
	let f = PI2 * i / (N - 1), v = a[0]
	for (let k = 1; k < a.length; k++) v += (k % 2 ? -1 : 1) * a[k] * cos(k * f)
	return v
}

/** Modified Bessel function of the first kind, order 0. Used by kaiser, kaiserBesselDerived. */
export function i0 (x) {
	let s = 1, t = 1
	for (let k = 1; k <= 25; k++) { t *= (x / (2 * k)) * (x / (2 * k)); s += t; if (t < 1e-15 * s) break }
	return s
}

/** Gegenbauer (ultraspherical) polynomial C_n^mu(x) via recurrence. Used by ultraspherical. */
export function gegen (n, mu, x) {
	if (n === 0) return 1
	if (n === 1) return 2 * mu * x
	let c0 = 1, c1 = 2 * mu * x
	for (let k = 2; k <= n; k++) {
		let c2 = (2 * x * (k + mu - 1) * c1 - (k + 2 * mu - 2) * c0) / k
		c0 = c1; c1 = c2
	}
	return c1
}

/**
 * Window from its spectrum sampled at N points, p[k] for frequencies πk/N: w = IDFT(p), centred, peak 1.
 * Even N shifts half a sample first, so both lengths centre alike. As SciPy's chebwin. Used by dolphChebyshev, ultraspherical.
 */
export function fromSpectrum (p) {
	let N = p.length, odd = N % 2, n = odd ? (N + 1) / 2 : N / 2 + 1, w = new Float64Array(N)
	for (let j = odd ? 0 : 1; j < n; j++) {
		let s = 0
		for (let k = 0; k < N; k++) s += p[k] * cos(PI * k * (odd ? 2 * j : 2 * j - 1) / N)
		if (odd) w[n - 1 + j] = w[n - 1 - j] = s
		else w[n - 2 + j] = w[n - 1 - j] = s
	}
	return normalize(w)
}

/** Chebyshev polynomial T_n(x) for any real x. Used by dolphChebyshev. */
export function chebyshev (n, x) {
	if (abs(x) <= 1) return cos(n * acos(x))
	let v = cosh(n * acosh(abs(x)))
	return x < 0 && n % 2 ? -v : v
}

/** Normalize array to peak absolute value of 1. Used by array-computed windows. */
export function normalize (w) {
	let peak = 0
	for (let i = 0; i < w.length; i++) if (abs(w[i]) > peak) peak = abs(w[i])
	if (peak > 0) for (let i = 0; i < w.length; i++) w[i] /= peak
	return w
}

// ── Public utilities ──

/** Generate a full window as Float64Array. */
export function generate (fn, N, ...params) {
	let w = new Float64Array(N)
	for (let i = 0; i < N; i++) w[i] = fn(i, N, ...params)
	return w
}

/** Apply window to a signal in-place. */
export function apply (signal, fn, ...params) {
	for (let i = 0, N = signal.length; i < N; i++) signal[i] *= fn(i, N, ...params)
	return signal
}

/** Equivalent noise bandwidth in frequency bins. Rectangular = 1.0, Hann ≈ 1.5. */
export function enbw (fn, N, ...params) {
	let s = 0, s2 = 0
	for (let i = 0; i < N; i++) { let v = fn(i, N, ...params); s += v; s2 += v * v }
	return N * s2 / (s * s)
}

/** Worst-case amplitude error in dB when a tone falls between DFT bins. */
export function scallopLoss (fn, N, ...params) {
	let s = 0, re = 0, im = 0
	for (let i = 0; i < N; i++) { let v = fn(i, N, ...params); s += v; re += v * cos(PI * i / N); im -= v * sin(PI * i / N) }
	return s === 0 ? Infinity : -20 * log10(sqrt(re * re + im * im) / abs(s))
}

/** COLA deviation. Returns max relative deviation from constant overlap-add sum; 0 = perfect. */
export function cola (fn, N, hop, ...params) {
	let win = generate(fn, N, ...params)
	let sums = new Float64Array(hop)
	for (let t = 0; t < hop; t++) for (let k = t; k < N; k += hop) sums[t] += win[k]
	let mean = 0
	for (let t = 0; t < hop; t++) mean += sums[t]
	mean /= hop
	if (mean === 0) return Infinity
	let maxDev = 0
	for (let t = 0; t < hop; t++) { let d = abs(sums[t] - mean) / mean; if (d > maxDev) maxDev = d }
	return maxDev
}