1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// Rayleigh-integral listener with the rim-diffraction back term, as in cymsim/radiation.py, but
13f4dbdexport .cymsim fileJeremy Magland 2// with the retarded-time sum evaluated at the audio rate: the normal acceleration of every node
3// is decimated to the audio rate inside the time loop (the same Kaiser-windowed sinc kernel as
4// the Python package, applied through DECIMATION_SPAN accumulator slots per node), the element
5// means of the finished samples are kept in a short history, and each element's history is read
6// at its fractional delay with a windowed-sinc interpolator of DELAY_TAPS taps. The per-element
7// coefficients, the delay taps and the decimation taps are computed here. See THEORY.md.
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 8import type { Listener } from "../config/schema.ts";
9import type { CymbalProfile } from "./geometry.ts";
10import type { Mesh } from "./mesh.ts";
13f4dbdexport .cymsim fileJeremy Magland 12/** Audio samples spanned by the decimation kernel (accumulator slots per node); as in cymsim/decimation.py. */
13export const DECIMATION_SPAN = 32;
14export const KAISER_BETA = 8.0;
15export const CUTOFF = 0.46;
16/** Taps of the fractional-delay interpolator, centred: lags i0 - DELAY_TAPS/2 + 1 ... i0 + DELAY_TAPS/2. */
17export const DELAY_TAPS = 16;
19/** Modified Bessel function I0 by its power series (converges quickly for x < 30). */
20export function besselI0(x: number): number {
21 let sum = 1, term = 1;
22 const y = (x * x) / 4;
23 for (let k = 1; k < 60; k++) { term *= y / (k * k); sum += term; if (term < 1e-17 * sum) break; }
24 return sum;
25}
27/** h(d) for d = -W m / 2 ... W m / 2 - 1, unit DC gain; identical to cymsim.decimation.decimation_kernel. */
28export function decimationKernel(m: number, W = DECIMATION_SPAN, beta = KAISER_BETA, cutoff = CUTOFF): Float64Array {
29 const L = W * m, half = L / 2, fc = cutoff / m, h = new Float64Array(L), i0 = besselI0(beta);
30 let sum = 0;
31 for (let i = 0; i < L; i++) {
32 const d = i - half, x = 2 * Math.PI * fc * d, sinc = d === 0 ? 1 : Math.sin(x) / x;
33 const u = d / half, win = besselI0(beta * Math.sqrt(Math.max(1 - u * u, 0))) / i0;
34 h[i] = sinc * win; sum += h[i];
35 }
36 for (let i = 0; i < L; i++) h[i] /= sum;
37 return h;
38}
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 39
40/**
13f4dbdexport .cymsim fileJeremy Magland 41 * Decimation taps in polyphase layout: taps[w * m + kk] is the weight with which the sub-step of
42 * phase kk (0..m-1) within audio sample j0 contributes to audio sample j0 + w - W/2 + 1, for slot
43 * w = 0..W-1. This is the (W, m) layout stored in a .cymsim file.
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 44 */
45export function decimationTaps(m: number, W = DECIMATION_SPAN): Float64Array {
13f4dbdexport .cymsim fileJeremy Magland 46 const h = decimationKernel(m, W), taps = new Float64Array(W * m);
47 for (let w = 0; w < W; w++) for (let kk = 0; kk < m; kk++) taps[w * m + kk] = h[kk + (W / 2 - 1 - w) * m + (W * m) / 2];
48 return taps;
49}
51/** Windowed-sinc taps that read x(j - i0 - fr) from x(j - i0 - t + T/2 - 1), t = 0..T-1 (Kaiser, beta 6, unit sum). */
52export function fractionalDelayTaps(fr: number, T = DELAY_TAPS): Float64Array {
53 const taps = new Float64Array(T), half = T / 2, i0 = besselI0(6);
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 54 let sum = 0;
56 const k = t - half + 1, x = k - fr; // lag offset relative to i0 (k = -T/2 + 1 ... T/2)
57 const sinc = x === 0 ? 1 : Math.sin(Math.PI * x) / (Math.PI * x);
58 const u = x / half, win = besselI0(6 * Math.sqrt(Math.max(1 - u * u, 0))) / i0;
59 taps[t] = sinc * win; sum += taps[t];
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 60 }
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 62 return taps;
63}
65export interface Radiation {
66 coefF: Float64Array; coefB: Float64Array; // rho A / (2 pi r), front and back
13f4dbdexport .cymsim fileJeremy Magland 67 i0F: Uint32Array; i0B: Uint32Array; // integer part of the delay in audio samples
68 frF: Float64Array; frB: Float64Array; // fractional part
69 tapsF: Float64Array; tapsB: Float64Array; // (E, DELAY_TAPS) interpolator weights
70 ringLength: number; // element history ring length in audio samples
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 71 taps: Float64Array; // (DECIMATION_SPAN, substeps) decimation weights
72 rFront: Float64Array; rBack: Float64Array;
73}
75export function buildRadiation(mesh: Mesh, profile: CymbalProfile, z0Node: Float64Array, listener: Listener, fs: number, substeps: number, nRim = 720): Radiation {
76 const E = mesh.e, [Lx, Ly, Lz] = listener.position;
77 const cx = new Float64Array(E), cy = new Float64Array(E), cz = new Float64Array(E);
78 for (let e = 0; e < E; e++) {
79 cx[e] = mesh.centroid[2 * e]; cy[e] = mesh.centroid[2 * e + 1];
80 cz[e] = (z0Node[mesh.tri[3 * e]] + z0Node[mesh.tri[3 * e + 1]] + z0Node[mesh.tri[3 * e + 2]]) / 3;
81 }
82 let R = 0;
83 for (let n = 0; n < mesh.n; n++) R = Math.max(R, Math.hypot(mesh.xy[2 * n], mesh.xy[2 * n + 1]));
84 const rimX = new Float64Array(nRim), rimY = new Float64Array(nRim), rimZ = new Float64Array(nRim), d2 = new Float64Array(nRim);
85 for (let j = 0; j < nRim; j++) {
86 const phi = (2 * Math.PI * j) / nRim;
87 rimX[j] = R * Math.cos(phi); rimY[j] = R * Math.sin(phi); rimZ[j] = profile.z0(rimX[j], rimY[j]);
88 d2[j] = Math.hypot(rimX[j] - Lx, rimY[j] - Ly, rimZ[j] - Lz);
89 }
90 const rFront = new Float64Array(E), rBack = new Float64Array(E);
91 const coefF = new Float64Array(E), coefB = new Float64Array(E), i0F = new Uint32Array(E), frF = new Float64Array(E);
13f4dbdexport .cymsim fileJeremy Magland 92 const i0B = new Uint32Array(E), frB = new Float64Array(E), tapsF = new Float64Array(E * DELAY_TAPS), tapsB = new Float64Array(E * DELAY_TAPS);
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 93 const { c_air: c, rho_air: rho, back_diffraction_gain: g } = listener;
94 let maxI0 = 0;
95 for (let e = 0; e < E; e++) {
96 rFront[e] = Math.hypot(Lx - cx[e], Ly - cy[e], Lz - cz[e]);
97 let best = Infinity;
98 for (let j = 0; j < nRim; j++) best = Math.min(best, Math.hypot(cx[e] - rimX[j], cy[e] - rimY[j], cz[e] - rimZ[j]) + d2[j]);
99 rBack[e] = best;
100 coefF[e] = (rho * mesh.area[e]) / (2 * Math.PI * rFront[e]);
101 coefB[e] = (g * rho * mesh.area[e]) / (2 * Math.PI * rBack[e]);
102 const dF = (rFront[e] / c) * fs, dB = (rBack[e] / c) * fs;
103 i0F[e] = Math.floor(dF); frF[e] = dF - i0F[e]; i0B[e] = Math.floor(dB); frB[e] = dB - i0B[e];
13f4dbdexport .cymsim fileJeremy Magland 104 tapsF.set(fractionalDelayTaps(frF[e]), e * DELAY_TAPS); tapsB.set(fractionalDelayTaps(frB[e]), e * DELAY_TAPS);
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 105 maxI0 = Math.max(maxI0, i0F[e], i0B[e]);
106 }
13f4dbdexport .cymsim fileJeremy Magland 107 return { coefF, coefB, i0F, i0B, frF, frB, tapsF, tapsB, ringLength: maxI0 + DELAY_TAPS + 4, taps: decimationTaps(substeps), rFront, rBack };
108}
110/** The delayed, interpolated value of an element history ring at sample j: sum_t taps[t] x(j - i0 - t + T/2 - 1). */
111export function delayedValue(hist: ArrayLike<number>, stride: number, e: number, D: number, j: number, i0: number, taps: ArrayLike<number>, tapBase: number): number {
112 let s = 0;
113 for (let t = 0; t < DELAY_TAPS; t++) {
114 const lag = i0 + t - DELAY_TAPS / 2 + 1, jj = j - lag;
115 if (jj >= 0 && lag >= 0) s += taps[tapBase + t] * hist[(jj % D) * stride + e];
116 }
117 return s;