1// The modal synthesis: a single mode struck by an impulse gives the closed
2// form damped sinusoid; the radiation gain of a uniform piston equals the
3// analytic value; the attack and the tail join; the sliced parts sum to the
4// whole.
5import { test } from "node:test";
6import assert from "node:assert/strict";
7import { makeGrid, defaultParams, type SimParams, listenerTerms, RHO_AIR, SPEED_OF_SOUND, nodeWeight } from "../src/params.ts";
8import { uniformPlate } from "../src/model/reference.ts";
9import { assemble } from "../src/model/assemble.ts";
10import { computeModes, localRunner } from "../src/modal/modes.ts";
11import { LN1000 } from "../src/loss.ts";
12import {
13 modalModel, modalRates, strikeGains, simulateContact, listenerWeights, radiationGains, tailStart, nodalAccelerationTable, attackPressure, modalTail, movieFrames, synthesize, sliceModel, sliceContact, type ModalModel,
14} from "../src/modal/synth.ts";
16const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
18async function smallModel(): Promise<{ model: ModalModel; pl: ReturnType<typeof uniformPlate> }> {
19 const g = makeGrid(0.2, 0.15, 0.0125);
20 const pl = uniformPlate(g, steel, 1e-3);
21 const sys = assemble(pl);
22 const modes = await computeModes(sys, 4000, localRunner(sys), { sliceSize: 20 });
23 return { model: modalModel(pl, sys, modes), pl };
24}
26test("a single mode under an impulse follows the closed-form damped sinusoid", async () => {
27 const { model, pl } = await smallModel();
28 const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.1, fps: 100 };
29 // Keep only mode 0 by slicing the model.
30 const one = sliceModel(model, 0, 1);
31 // Give the mode a decay time of exactly 1 s through the material term alone.
32 one.sigmaMat[0] = LN1000; one.sigmaTe[0] = 0; one.sigmaRad[0] = 0;
33 const sigma = modalRates(one);
34 const g = strikeGains(one, pl, p);
35 const c = simulateContact(one, sigma, g, p, p.fps);
36 assert.ok(c.peakForce > 0 && c.contactTime > 0);
37 // After tEnd the mode is q(t) = Re[A exp(s (t - tEnd))]; check the movie frames against it.
38 const nFrames = 10;
39 const frames = movieFrames(one, sigma, c, p.fps, nFrames);
40 const w = 2 * Math.PI * one.freqs[0], s = sigma[0], wd = Math.sqrt(w * w - s * s);
41 for (let f = 1; f < nFrames; f++) {
42 const t = f / p.fps - c.tEnd;
43 if (t < 0) continue;
44 const q = Math.exp(-s * t) * (c.aRe[0] * Math.cos(wd * t) - c.aIm[0] * Math.sin(wd * t));
45 // Compare at a node away from the edge.
46 const k = Math.floor(pl.grid.ny / 2) * pl.grid.nx + Math.floor(pl.grid.nx / 3);
47 const phi = one.wShapes[one.wRow[k] * one.count];
48 assert.ok(Math.abs(frames[f * one.nGrid + k] - phi * q) < 1e-9 + 1e-6 * Math.abs(phi * q));
49 }
50 assert.ok(Math.abs(sigma[0] - LN1000 / 1.0) < 1e-9, "decay rate from T60 = 1 s");
51});
53test("the radiation gain of a uniform piston equals the analytic Rayleigh sum", async () => {
54 const { model, pl } = await smallModel();
55 const p: SimParams = { ...defaultParams, listenerDist: 1.0, listenerAngle: 0 };
56 // A fake model with one mode of uniform shape (mass-normalized to 1/sqrt(M)).
57 const uniform: ModalModel = { ...sliceModel(model, 0, 1), freqs: new Float64Array([1000]) };
58 let M = 0;
59 for (let k = 0; k < uniform.nGrid; k++) M += uniform.masses[k];
60 const phi = 1 / Math.sqrt(M);
61 for (let r = 0; r < uniform.nRows; r++) uniform.wShapes[r] = phi;
62 const { weight, delay } = listenerWeights(pl, p);
63 const sigma = new Float64Array([0]);
64 const gains = radiationGains(uniform, sigma, weight, delay);
65 // Analytic: G = sum_k rho dA / (2 pi r_k) phi exp(-i omega tau_k).
66 const lt = listenerTerms(pl, p);
67 let dmin = Infinity;
68 for (let k = 0; k < lt.weight.length; k++) if (lt.weight[k] > 0) dmin = Math.min(dmin, lt.delay[k]);
69 let re = 0, im = 0;
70 const w = 2 * Math.PI * 1000;
71 for (let k = 0; k < lt.weight.length; k++) {
72 if (lt.weight[k] === 0) continue;
73 const tau = lt.delay[k] - dmin;
74 re += lt.weight[k] * phi * Math.cos(w * tau);
75 im -= lt.weight[k] * phi * Math.sin(w * tau);
76 }
77 // The shapes are float32, so agreement is to about 1e-7.
78 assert.ok(Math.abs(gains.re[0] - re) < 1e-6 * Math.abs(re) && Math.abs(gains.im[0] - im) < 1e-6 * Math.max(Math.abs(im), 1e-6));
79 // On axis at 1 m the weights sum to rho A / (2 pi r) within the delay spread.
80 let area = 0;
81 for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) area += nodeWeight(pl.grid, i, j) * pl.grid.hx * pl.grid.hy;
82 const approx = (RHO_AIR * area) / (2 * Math.PI * 1.0) * phi;
83 assert.ok(Math.abs(Math.hypot(re, im) / approx - 1) < 0.05);
84 void SPEED_OF_SOUND;
85});
87test("the attack and the tail join, and the sliced parts sum to the whole", async () => {
88 const { model, pl } = await smallModel();
89 const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.2, fps: 120, strikeVelocity: 2 };
90 const whole = synthesize(model, pl, p, true);
91 const nSamples = Math.round(p.duration * p.fs);
92 assert.equal(whole.audio.length, nSamples);
93 // No discontinuity at the tail start beyond the local sample-to-sample variation.
94 const { maxDelay, weight, delay } = listenerWeights(pl, p);
95 const n0 = Math.ceil(tailStart(whole.contact, maxDelay) * p.fs - 1e-9);
96 let maxStep = 0, jump = Math.abs(whole.audio[n0] - whole.audio[n0 - 1]);
97 for (let n = n0 + 2; n < n0 + 200; n++) maxStep = Math.max(maxStep, Math.abs(whole.audio[n] - whole.audio[n - 1]));
98 assert.ok(jump < 3 * maxStep, `jump ${jump} vs typical step ${maxStep}`);
99 // The tail decays: the last 20 ms is quieter than the first 20 ms.
100 const rms = (a: Float32Array) => Math.sqrt(a.reduce((s, v) => s + v * v, 0) / a.length);
101 assert.ok(rms(whole.audio.subarray(nSamples - 960)) < rms(whole.audio.subarray(0, 960)));
102 // Split the modes into three slices and sum the parts.
103 const sigma = modalRates(model);
104 const g = strikeGains(model, pl, p);
105 const c = simulateContact(model, sigma, g, p, p.fps);
106 const tStart = tailStart(c, maxDelay);
107 const len = Math.min(c.tableLen, Math.ceil(tStart * c.tableRate) + 2);
108 const audio = new Float32Array(nSamples);
109 const table = new Float32Array(model.nRows * len);
110 const nFrames = Math.round(p.duration * p.fps);
111 const frames = new Float32Array(nFrames * model.nGrid);
112 const cuts = [0, Math.floor(model.count / 3), Math.floor((2 * model.count) / 3), model.count];
113 for (let s = 0; s < 3; s++) {
114 const sub = sliceModel(model, cuts[s], cuts[s + 1]);
115 const cc = sliceContact(c, cuts[s], cuts[s + 1]);
116 const sg = sigma.slice(cuts[s], cuts[s + 1]);
117 nodalAccelerationTable(sub, cc, len, table);
118 modalTail(sub, sg, cc, radiationGains(sub, sg, weight, delay), p.fs, tStart, nSamples, audio);
119 movieFrames(sub, sg, cc, p.fps, nFrames, frames);
120 }
121 attackPressure(model, table, len, c.tableRate, weight, delay, p.fs, tStart, audio);
122 let err = 0, nrm = 0;
123 for (let n = 0; n < nSamples; n++) { err += (audio[n] - whole.audio[n]) ** 2; nrm += whole.audio[n] ** 2; }
124 assert.ok(Math.sqrt(err / nrm) < 1e-5, `audio split error ${Math.sqrt(err / nrm)}`);
125 err = 0; nrm = 0;
126 for (let i = 0; i < frames.length; i++) { err += (frames[i] - whole.frames[i]) ** 2; nrm += whole.frames[i] ** 2; }
127 assert.ok(Math.sqrt(err / nrm) < 1e-5, `frame split error ${Math.sqrt(err / nrm)}`);
128});
130test("a lighter hit gives a proportionally smaller signal in the linear model", async () => {
131 const { model, pl } = await smallModel();
132 const p1: SimParams = { ...defaultParams, fs: 48000, duration: 0.05, strikeVelocity: 0.2 };
133 const p2: SimParams = { ...p1, strikeVelocity: 0.1 };
134 const a1 = synthesize(model, pl, p1, false).audio, a2 = synthesize(model, pl, p2, false).audio;
135 const peak = (a: Float32Array) => a.reduce((m, v) => Math.max(m, Math.abs(v)), 0);
136 // The Hertz contact is nonlinear (a faster hit is shorter and brighter), so the ratio exceeds 2 somewhat.
137 const ratio = peak(a1) / peak(a2);
138 assert.ok(ratio > 1.8 && ratio < 3.5, `ratio ${ratio}`);
139});