ae42884Anti-aliased pickup: solve at R fs and low-pass filter the pickup down to fs (both solvers); CACHE_VERSION 3; chaos-check analysis scriptsJeremy Magland 1import { timeStep } from './solver';
2import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
4// Anti-aliased pickup. The solvers record the pickup once per output sample (point
5// sampling), so motion above fs/2 would fold back into the audible band: a 10 mm pluck of
6// the 5-bead string puts about 28% of its pickup power above 24 kHz. Instead, a run is
7// solved with an internal sample rate R fs, and the pickup is low-pass filtered (a Kaiser-
8// windowed sinc, passband to 0.45 fs, stopband from 0.5 fs, about 86 dB down) and kept
9// every R-th sample. R is chosen so that the internal Nyquist frequency is above twice the
10// system's highest linear frequency plus fs/2 (room for nonlinear combination tones),
11// capped at the solver's steps per sample, so the time step does not shrink.
13export function oversampling(sys: BeadSystem, exc: Excitation, opt: SimOptions): number {
14 const { omegaMax, stepsPerSample } = timeStep(sys, exc.mallet, opt.fs, opt.cfl);
15 const fmax = omegaMax / (2 * Math.PI);
16 return Math.max(1, Math.min(stepsPerSample, Math.ceil((2 * fmax + opt.fs / 2) / opt.fs)));
17}
19// the options the solver runs with: R times the sample rate, delays in internal samples
20export function oversampledOptions(opt: SimOptions, R: number): SimOptions {
21 return { ...opt, fs: opt.fs * R, pickupDelay: opt.pickupDelay?.map((d) => d * R) };
22}
24const BETA = 8.6; // Kaiser window: about 86 dB stopband
26// filter taps at the internal rate, symmetric about the middle, summing to 1
27const kernels = new Map<number, Float64Array>();
28export function lowpass(R: number): Float64Array {
29 let h = kernels.get(R);
30 if (h) return h;
31 const df = 0.05 / R; // transition width, cycles per internal sample (0.45 to 0.5 of fs)
32 const fc = 0.475 / R; // cutoff, cycles per internal sample
33 let N = Math.ceil((86 - 7.95) / (14.36 * df)) + 1;
34 if (N % 2 === 0) N++;
35 const M = (N - 1) / 2;
36 h = new Float64Array(N);
37 const i0b = besselI0(BETA);
38 let sum = 0;
39 for (let i = 0; i < N; i++) {
40 const x = i - M;
41 const sinc = x === 0 ? 2 * fc : Math.sin(2 * Math.PI * fc * x) / (Math.PI * x);
42 const r = x / M;
43 h[i] = sinc * (besselI0(BETA * Math.sqrt(Math.max(0, 1 - r * r))) / i0b);
44 sum += h[i];
45 }
46 for (let i = 0; i < N; i++) h[i] /= sum;
47 kernels.set(R, h);
48 return h;
49}
51function besselI0(x: number): number {
52 let s = 1;
53 let term = 1;
54 for (let k = 1; k < 50; k++) {
55 term *= (x / (2 * k)) ** 2;
56 s += term;
57 if (term < 1e-12 * s) break;
58 }
59 return s;
60}
62// Filter and decimate an internal-rate result to fs. complete: the whole run (the signal
63// is zero after its end); otherwise a partial result, of which only the output samples
64// whose filter window is covered are kept.
65export function decimateResult(r: SimResult, R: number, complete: boolean): SimResult {
66 if (R === 1) return { ...r, oversampling: 1 };
67 const h = lowpass(R);
68 const M = (h.length - 1) / 2;
69 const x = r.audio;
70 const nOut = complete ? Math.ceil(x.length / R) : Math.max(0, Math.floor((x.length - 1 - M) / R) + 1);
71 const audio = new Float32Array(nOut);
72 for (let i = 0; i < nOut; i++) {
73 const c = i * R; // centre, in internal samples; the signal is zero before t = 0
74 const a = Math.max(0, c - M);
75 const b = Math.min(x.length - 1, c + M);
76 let s = 0;
77 for (let j = a; j <= b; j++) s += h[j - c + M] * x[j];
78 audio[i] = s;
79 }
80 let contactForce: Float32Array | undefined;
81 if (r.contactForce) {
82 contactForce = new Float32Array(nOut);
83 for (let i = 0; i < nOut; i++) contactForce[i] = r.contactForce[Math.min(i * R, r.contactForce.length - 1)];
84 }
85 return { ...r, audio, contactForce, fs: r.fs / R, stepsPerSample: r.stepsPerSample * R, oversampling: R };
86}