982850fWoodblock strike simulator: Python reference model and interactive web appJeremy Magland 1// Anti-aliased decimation: Kaiser-windowed-sinc FIR lowpass at ~0.45 * fs_lo,
2// evaluated only at the retained output samples (polyphase, so the cost is
3// taps * T / R rather than taps * T). Replaces scipy.resample_poly(p, 1, R).
5function besselI0(x: number): number {
6 // Series expansion; converges quickly for the argument range used here.
7 let sum = 1.0;
8 let term = 1.0;
9 const half = x / 2.0;
10 for (let k = 1; k < 64; k++) {
11 term *= (half / k) * (half / k);
12 sum += term;
13 if (term < 1e-16 * sum) break;
14 }
15 return sum;
16}
18export function designDecimator(R: number, beta = 8.0): Float64Array {
19 // Cutoff at 0.45 of the output Nyquist-setting rate: fc/fs_hi = 0.45/R.
20 const taps = 12 * R + 1; // odd length, linear phase
21 const h = new Float64Array(taps);
22 const M2 = (taps - 1) / 2;
23 const fcNorm = 0.45 / R; // cycles per input sample
24 const i0beta = besselI0(beta);
25 let sum = 0;
26 for (let i = 0; i < taps; i++) {
27 const m = i - M2;
28 const sinc = m === 0 ? 2 * fcNorm : Math.sin(2 * Math.PI * fcNorm * m) / (Math.PI * m);
29 const r = m / M2;
30 const w = besselI0(beta * Math.sqrt(Math.max(0, 1 - r * r))) / i0beta;
31 h[i] = sinc * w;
32 sum += h[i];
33 }
34 for (let i = 0; i < taps; i++) h[i] /= sum; // unity DC gain
35 return h;
36}
38export function decimate(p: Float64Array, R: number): Float64Array {
39 if (R === 1) return p.slice();
40 const h = designDecimator(R);
41 const taps = h.length;
42 const M2 = (taps - 1) / 2;
43 const T = p.length;
44 const outLen = Math.floor(T / R);
45 const y = new Float64Array(outLen);
46 for (let m = 0; m < outLen; m++) {
47 const centre = m * R;
48 let acc = 0;
49 // y[m] = sum_i h[i] * p[centre + M2 - i], zero outside the buffer
50 const iLo = Math.max(0, centre + M2 - (T - 1));
51 const iHi = Math.min(taps - 1, centre + M2);
52 for (let i = iLo; i <= iHi; i++) acc += h[i] * p[centre + M2 - i];
53 y[m] = acc;
54 }
55 return y;
56}