1// Losses given as decay times. The design script gives T60, the time in
2// which a partial decays by 60 dB, as a number, a table [f T60] or a
3// function of frequency; the design layer samples it on LOSS_FREQS, and this
4// module interpolates it, converts it to decay rates for the modes, and fits
5// the two-parameter model of the full simulation to it.
7import { type Loss, type Material, bendingConstant } from "./params.ts";
9/** Frequencies at which T60 is sampled: 128 points from 10 Hz to 40 kHz, log-spaced. */
10export const LOSS_FREQS: Float64Array = (() => {
11 const n = 128, f = new Float64Array(n);
12 for (let i = 0; i < n; i++) f[i] = 10 * Math.pow(4000, i / (n - 1));
13 return f;
14})();
16/** ln(1000): a rate sigma [1/s] and a decay time relate by sigma = LN1000 / T60. */
17export const LN1000 = Math.log(1000);
19/** The default losses: 4 s at low frequencies, falling to about 1 s at 5 kHz. */
20export function defaultT60(f: number): number {
21 return 4 / (1 + f / 2000);
22}
24export function lossFromFunction(fn: (f: number) => number): Loss {
25 const t60 = new Float64Array(LOSS_FREQS.length);
26 for (let i = 0; i < t60.length; i++) t60[i] = fn(LOSS_FREQS[i]);
27 return { freqs: LOSS_FREQS, t60 };
28}
30export function defaultLoss(): Loss {
31 return lossFromFunction(defaultT60);
32}
34export function constantLoss(t60: number): Loss {
35 return lossFromFunction(() => t60);
36}
38/** Log-log interpolation of a table of (frequency, T60) rows, clamped at the ends. */
39export function lossFromTable(freqs: ArrayLike<number>, t60s: ArrayLike<number>): Loss {
40 const n = freqs.length;
41 const idx = Array.from({ length: n }, (_, i) => i).sort((a, b) => freqs[a] - freqs[b]);
42 const fs = idx.map((i) => freqs[i]), ts = idx.map((i) => t60s[i]);
43 return lossFromFunction((f) => interpLogLog(fs, ts, f));
44}
46function interpLogLog(fs: number[], ts: number[], f: number): number {
47 const n = fs.length;
48 if (n === 1 || f <= fs[0]) return ts[0];
49 if (f >= fs[n - 1]) return ts[n - 1];
50 let i = 0;
51 while (i < n - 2 && fs[i + 1] < f) i++;
52 const t = (Math.log(f) - Math.log(fs[i])) / (Math.log(fs[i + 1]) - Math.log(fs[i]));
53 return Math.exp(Math.log(ts[i]) * (1 - t) + Math.log(ts[i + 1]) * t);
54}
56/** T60 at a frequency, by log-log interpolation of the sampled curve, clamped at its ends. */
57export function t60At(loss: Loss, f: number): number {
58 return interpLogLog(Array.from(loss.freqs), Array.from(loss.t60), f);
59}
61/** Decay rate sigma [1/s] of a partial at frequency f: its amplitude decays as exp(-sigma t). */
62export function decayRate(loss: Loss, f: number): number {
63 return LN1000 / t60At(loss, f);
64}
66/**
67 * The two-parameter loss model of the full simulation. A bending mode of
68 * angular frequency omega decays at the rate sigma0 + sigma1 k^2 with
69 * k^2 = omega / kappa; the fit is a non-negative least-squares fit of that
70 * line to the design's decay rates on log-spaced frequencies from 50 Hz to
71 * 10 kHz, at the mean thickness of the plate.
72 */
73export interface SigmaFit {
74 sigma0: number; // [1/s]
75 sigma1: number; // [m^2/s]
76 kappa: number; // [m^2/s] bending constant used for the fit
77}
79export function fitSigma(loss: Loss, mat: Material, hMean: number, fLo = 50, fHi = 10000, n = 40): SigmaFit {
80 const kappa = bendingConstant(mat, hMean);
81 const x: number[] = [], y: number[] = [];
82 for (let i = 0; i < n; i++) {
83 const f = fLo * Math.pow(fHi / fLo, i / (n - 1));
84 x.push((2 * Math.PI * f) / kappa); // k^2
85 y.push(decayRate(loss, f));
86 }
87 // Unconstrained least squares, then clamp: with two parameters the
88 // non-negative solution is the unconstrained one, or one parameter at
89 // zero with the other refit.
90 let sx = 0, sy = 0, sxx = 0, sxy = 0;
91 for (let i = 0; i < n; i++) { sx += x[i]; sy += y[i]; sxx += x[i] * x[i]; sxy += x[i] * y[i]; }
92 const det = n * sxx - sx * sx;
93 let s1 = (n * sxy - sx * sy) / det, s0 = (sy - s1 * sx) / n;
94 if (s1 < 0) { s1 = 0; s0 = sy / n; }
95 if (s0 < 0) { s0 = 0; s1 = sxy / sxx; }
96 return { sigma0: s0, sigma1: s1, kappa };
97}
99/** Decay rate of the fitted model at a frequency. */
100export function fittedRate(fit: SigmaFit, f: number): number {
101 return fit.sigma0 + (fit.sigma1 * 2 * Math.PI * f) / fit.kappa;
102}