2// Checks the MATLAB layer against the numbers the chapters claim, through the
3// numbl CLI, with no browser involved.
4//
5// npm run test:matlab (NUMBL_DIR defaults to ~/src/numbl)
6//
7// Every assertion here is a statement made somewhere in the prose, so if one
8// of these fails either the code or the text is wrong, and the text is not
9// allowed to win by default.
10//
11// Almost everything measured in this tutorial is deterministic, so most
12// tolerances are tight: a quoted number is checked to the precision it is
13// quoted to. The exceptions are wall-clock times, which are never asserted
14// as values, only as orderings, and the reverberation times of Chapter 8,
15// which depend on one source-receiver pair and are quoted to two figures.
16import { join, dirname } from 'node:path'
17import { fileURLToPath } from 'node:url'
18import { runLab } from './runlab.mjs'
20const root = join(dirname(fileURLToPath(import.meta.url)), '..')
21const src = (p) => join(root, 'src', 'chapters', p)
22const ONLY = process.argv.includes('--only') ? process.argv[process.argv.indexOf('--only') + 1] : null
24let failures = 0
25const check = (name, ok, detail = '') => {
26 console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
27 if (!ok) failures++
28}
29const near = (a, b, tol) => Number.isFinite(a) && Math.abs(a - b) <= tol
30/** within `frac` of b, relatively */
31const rel = (a, b, frac) => Number.isFinite(a) && Math.abs(a - b) <= frac * Math.abs(b)
32const note = (m) => 440 * 2 ** ((m - 69) / 12)
33const maxAbs = (a) => Math.max(...a.filter((v) => v != null).map(Math.abs))
34const section = (id, title) => {
35 if (ONLY && ONLY !== id) return false
36 console.log(`\n${title}`)
37 return true
38}
40// ---------------------------------------------------------------- Chapter 1
41if (section('ch1', 'Chapter 1: an exact string, and aliasing')) {
42 const D = src('ch01/driver.m')
43 const S = src('ch01/pluck.m')
44 const pl = (where) =>
45 runLab(D, S, { mode: 'pluck', fs: 44100, f0: 110, beta: 0.2, where, xp: 0.12, T: 2.5, T60: 3, fmax: 6000 }).out
46 const b = pl('bridge')
47 check('200 partials below Nyquist at 110 Hz', b.nyq === 200, String(b.nyq))
48 check('12 of the first 60 missing at beta = 0.2', b.missing === 12, String(b.missing))
49 check('band-limited vs dAlembert, bridge: 5.0e-2', rel(b.reldiff, 0.05, 0.02), b.reldiff.toExponential(3))
50 const p = pl('pickup')
51 check('band-limited vs dAlembert, pickup: 3.5e-4', rel(p.reldiff, 3.5e-4, 0.03), p.reldiff.toExponential(3))
52 const table = {
53 bridge: { 45: -26.0, 57: -22.0, 69: -18.9, 81: -15.9, 93: -12.7, 100: -11.4 },
54 pickup: { 45: -69.2, 57: -60.0, 69: -51.0, 81: -42.0, 93: -32.0, 100: -28.6 },
55 }
56 for (const where of ['bridge', 'pickup']) {
57 for (const [m, want] of Object.entries(table[where])) {
58 const { out } = runLab(D, S, { mode: 'alias', fs: 44100, f0: note(Number(m)), beta: 0.2, where, xp: 0.12, T: 1.5, T60: 1.2 })
59 check(`${where}, note ${m}: folded energy ${want} dB`, near(out.measured, want, 0.06), out.measured.toFixed(2))
60 if (Number(m) !== 45) {
61 check(` and the prediction agrees to 0.03 dB`, near(out.measured, out.predicted, 0.035), (out.measured - out.predicted).toFixed(3))
62 } else {
63 check(` and at A2 they differ by about 1 dB (bridge) or less`, Math.abs(out.measured - out.predicted) < 1.2, (out.measured - out.predicted).toFixed(2))
64 }
65 }
66 }
67}
69// ---------------------------------------------------------------- Chapter 2
70if (section('ch2', 'Chapter 2: leapfrog and the magic time step')) {
71 const D = src('ch02/driver.m')
72 const S = src('ch02/leapfrog.m')
73 const run = (N) => runLab(D, S, { mode: 'run', fs: 44100, f0: 44100 / 400, N, T: 2.5, where: 'bridge', beta: 0.2, xp: 0.12 }).out
74 const r200 = run(200)
75 check('N = 200 is exact to ~1e-12', r200.relerr < 1e-11 && r200.relerr > 1e-14, r200.relerr.toExponential(2))
76 const r100 = run(100)
77 check('N = 100: relative error 0.72', near(r100.relerr, 0.72, 0.005), r100.relerr.toFixed(3))
78 check('N = 100: partial 1 is 0.05 cents flat', near(r100.pred1, -0.053, 0.002), r100.pred1.toFixed(4))
79 check('N = 100: partial 10 is 5.4 cents flat', near(r100.pred10, -5.36, 0.01), r100.pred10.toFixed(3))
80 const r50 = run(50)
81 check('N = 50: error 1.04, partial 10 is 27 cents flat', near(r50.relerr, 1.04, 0.01) && near(r50.pred10, -26.8, 0.1),
82 `${r50.relerr.toFixed(3)}, ${r50.pred10.toFixed(2)}`)
83 const r201 = run(201)
84 check('N = 201 blows up after 69 steps, 1.6 ms', r201.unstable && near(r201.blowup, (69 / 44100) * 1000, 1e-9), String(r201.blowup))
86 const { out: c } = runLab(D, S, { mode: 'cents', f0: 44100 / 400, N: 100, lams: [0.25, 0.5, 0.75, 0.9, 1] })
87 check('measured and predicted cents agree below 1e-7', c.agree < 1e-7, c.agree.toExponential(2))
88 check('at lambda = 1 every error is below 1e-7 cents', maxAbs(c.meas[4]) < 1e-7, maxAbs(c.meas[4]).toExponential(2))
89 // where each predicted curve crosses -5 cents, in points per wavelength
90 const cross = (k) => {
91 const pred = c.pred[k]
92 for (let n = 1; n < pred.length; n++) if (pred[n] < -5) return (2 * c.N) / (n + 0.5)
93 return NaN
94 }
95 for (const [k, want] of [[0, 23], [1, 21], [2, 16], [3, 11]]) {
96 check(`5 cents at lambda = ${c.lams[k]} needs about ${want} points per wavelength`, near(cross(k), want, 1.2), cross(k).toFixed(1))
97 }
98 check('near two points per wavelength the error is several hundred cents',
99 c.pred[1][c.pred[1].length - 1] < -300, c.pred[1][c.pred[1].length - 1].toFixed(0))
101 const click = (lam) => runLab(D, S, { mode: 'click', fs: 44100, c: 343, L: 17.15, lam, T: 1.2, sigt: 3e-5 }).out
102 const k1 = click(1)
103 check('the click at lambda = 1 is exact to ~3e-11', k1.rel < 1e-10, k1.rel.toExponential(2))
104 const k5 = click(0.5)
105 check('at lambda = 0.5 the grid carries up to 7.35 kHz', near(k5.fcut, 7350, 10), k5.fcut.toFixed(0))
106 check('and the waveform error is 1.7', near(k5.rel, 1.7, 0.05), k5.rel.toFixed(3))
107 const k25 = click(0.25)
108 check('at lambda = 0.25, up to 3.55 kHz', near(k25.fcut, 3550, 10), k25.fcut.toFixed(0))
109 // the group velocity quoted in the text
110 const w = 2 * Math.PI * 5000 / 44100
111 const kh = 2 * Math.asin(Math.sin(w / 2) / 0.5)
112 const vg = Math.cos(kh / 2) / Math.cos(w / 2)
113 check('group velocity at 5 kHz is about three quarters of c', near(vg, 0.765, 0.01), vg.toFixed(3))
114 check('so after 334 m it is about 0.3 s late', near((334 / 343) * (1 / vg - 1), 0.3, 0.01), ((334 / 343) * (1 / vg - 1)).toFixed(3))
115}
117// ---------------------------------------------------------------- Chapter 3
118if (section('ch3', 'Chapter 3: the digital waveguide')) {
119 const D = src('ch03/driver.m')
120 const S = src('ch03/ks.m')
121 const { out: e } = runLab(D, S, { mode: 'equiv', fs: 44100, T: 1, beta: 0.2, Ns: [25, 100, 400, 1600, 3200] })
122 check('the loop and leapfrog agree to ~1e-11', Math.max(...e.diff) < 1e-10, Math.max(...e.diff).toExponential(2))
123 check('leapfrog time grows with N; the loop does not',
124 e.tl[4] > 10 * e.tl[1] && e.tw[4] < 5 * Math.max(e.tw[1], 0.5), `${e.tl.map((v) => v.toFixed(0))} / ${e.tw.map((v) => v.toFixed(1))}`)
125 check('at N = 3200 the loop is several hundred times faster', e.tl[4] / e.tw[4] > 200, (e.tl[4] / e.tw[4]).toFixed(0))
127 const notes = Array.from({ length: 25 }, (_, i) => 36 + 3 * i)
128 const { out: k } = runLab(D, S, { mode: 'keys', fs: 44100, beta: 0.2, notes, Tn: 0.5, scale: [84, 86], step: 0.22, ring: 0.9 })
129 const worst = (row) => row.reduce((m, v) => (Math.abs(v) > Math.abs(m) ? v : m), 0)
130 const lowNotes = notes.map((m, i) => (note(m) < 262 ? Math.abs(k.cents[0][i]) : 0))
131 check('whole-sample loop within 2 cents up to about 260 Hz', Math.max(...lowNotes) < 2, Math.max(...lowNotes).toFixed(2))
132 const iw = k.cents[0].findIndex((v) => v === worst(k.cents[0]))
133 check('whole-sample loop worst 47 cents, at F#7', near(worst(k.cents[0]), 46.9, 0.2) && notes[iw] === 102, `${worst(k.cents[0]).toFixed(2)} at ${notes[iw]}`)
134 check('linear interpolation worst 0.14 cents', near(Math.abs(worst(k.cents[1])), 0.14, 0.01), worst(k.cents[1]).toFixed(3))
135 check('allpass worst 0.63 cents', near(Math.abs(worst(k.cents[2])), 0.63, 0.01), worst(k.cents[2]).toFixed(3))
136 const t60 = (m) => k.t60[0][notes.indexOf(m)]
137 check('C2 rings about 9700 s', rel(t60(36), 9740, 0.02), t60(36).toFixed(0))
138 check('A4 32 s, C6 2.4 s, C8 36 ms', rel(t60(69), 32.2, 0.02) && rel(t60(84), 2.43, 0.02) && rel(t60(108), 0.0362, 0.03),
139 `${t60(69).toFixed(1)}, ${t60(84).toFixed(2)}, ${t60(108).toFixed(4)}`)
140 const pred = (f) => -3 / (f * Math.log10(Math.cos((Math.PI * f) / 44100)))
141 check('and each agrees with the loop-gain formula to a few per cent', rel(t60(69), pred(note(69)), 0.02) && rel(t60(108), pred(note(108)), 0.05),
142 `${pred(note(69)).toFixed(1)}, ${pred(note(108)).toFixed(4)}`)
143 const ratios = k.t60[1].map((v, i) => v / k.t60[0][i])
144 check('linear interpolation cuts decay time by up to about a half', Math.min(...ratios) > 0.44 && Math.min(...ratios) < 0.52, Math.min(...ratios).toFixed(2))
145}
147// ------------------------------------------------------------ Chapters 4, 5
148const MEM = src('membrane/driver.m')
149const DRUM = { c: 150 * Math.SQRT2, kmax: 90, width: 0.03, xr: 0.77, yr: 0.63, fmax: 3000, T: 1.0, xs: 0.31, ys: 0.42 }
150if (section('ch4', 'Chapter 4: a square drum')) {
151 const S = src('ch04/drum.m')
152 const st = (N, lam) => runLab(MEM, S, { mode: 'strike', ...DRUM, N, lam, scheme: 'fd2' }).out
153 const s48 = st(48, 0.7)
154 check('N = 48: 14.5 kHz, error 0.59, about 2600 modes', near(s48.fs, 14546, 1) && near(s48.rel, 0.59, 0.005) && Math.abs(s48.nmodes - 2600) < 50,
155 `${s48.fs.toFixed(0)}, ${s48.rel.toFixed(3)}, ${s48.nmodes}`)
156 check('N = 64: 0.54; N = 24: 0.80', near(st(64, 0.7).rel, 0.54, 0.006) && near(st(24, 0.7).rel, 0.80, 0.006))
157 const { out: c } = runLab(MEM, S, { mode: 'cents', N: 64, lams: [0.3, 0.5, Math.SQRT1_2], rs: [8, 16, 32], nang: 9 })
158 check('measured and predicted agree to better than 1e-6 cents', c.agree < 1e-6, c.agree.toExponential(2))
159 const at = (m, n, k) => c.meas[k][c.m.findIndex((v, i) => v === m && c.n[i] === n)]
160 check('16 ppw, axis: -5.4 at 1/sqrt2, -10.0 at 0.3', near(at(8, 1, 2), -5.35, 0.06) && near(at(8, 1, 0), -9.96, 0.06),
161 `${at(8, 1, 2).toFixed(2)}, ${at(8, 1, 0).toFixed(2)}`)
162 check('16 ppw, diagonal: 0 and -5.1', near(at(6, 6, 2), 0, 1e-6) && near(at(6, 6, 0), -5.14, 0.06), `${at(6, 6, 2).toExponential(1)}, ${at(6, 6, 0).toFixed(2)}`)
163 check('8 ppw axis -23, 4 ppw axis -101', near(at(16, 1, 2), -22.7, 0.2) && near(at(32, 2, 2), -100.7, 0.5), `${at(16, 1, 2).toFixed(1)}, ${at(32, 2, 2).toFixed(1)}`)
164 // the leading-order formula, -2847 (cos^4 + sin^4 - lam^2) / P^2, at each
165 // mode's own angle and wavelength; the spread between an axis mode and a
166 // diagonal one then comes out the same at both Courant numbers
167 const lead = (m, n, lam) => {
168 const th = Math.atan2(n, m)
169 const P = (2 * 64) / Math.hypot(m, n)
170 return (-2847 * (Math.cos(th) ** 4 + Math.sin(th) ** 4 - lam * lam)) / (P * P)
171 }
172 for (const [k, lam] of [[0, 0.3], [2, Math.SQRT1_2]]) {
173 const got = at(6, 6, k) - at(8, 1, k)
174 const want = lead(6, 6, lam) - lead(8, 1, lam)
175 check(`the direction spread at lambda = ${lam.toFixed(3)} follows the leading-order formula`, rel(got, want, 0.03), `${got.toFixed(2)} vs ${want.toFixed(2)}`)
176 }
177 const P5 = Math.sqrt((2847 * 0.5) / 5)
178 check('5 cents along the axes needs about 17 points per wavelength', near(P5, 17, 0.5), P5.toFixed(1))
179}
181if (section('ch5', 'Chapter 5: fewer points per wavelength')) {
182 const S = src('ch05/drum.m')
183 const { out: o } = runLab(MEM, S, {
184 mode: 'orders', Ns: [16, 24, 32, 48, 64, 96], schemes: ['fd2', 'fd4', 'fd44'], lams: [0.7, 0.6, 0.7], modes: [4, 3, 9, 2, 14, 1],
185 })
186 const w = o.worst
187 const last = 5
188 check('second order: 409 cents at 16, 7.7 at 96, order 2.03', near(w[0][0], 409, 1) && near(w[0][last], 7.69, 0.02) && near(o.orders[0][last], 2.03, 0.02),
189 `${w[0][0].toFixed(0)}, ${w[0][last].toFixed(2)}, ${o.orders[0][last].toFixed(3)}`)
190 check('fourth order in both: 245 to 0.19, order 4.0', near(w[2][0], 245, 1) && near(w[2][last], 0.187, 0.003) && o.orders[2].slice(1).every((q) => near(q, 4.0, 0.1)),
191 `${w[2][0].toFixed(0)}, ${w[2][last].toFixed(3)}, ${o.orders[2].slice(1).map((q) => q.toFixed(2))}`)
192 check('fourth order in both at 48: 3.0 cents, below any second-order grid', near(w[2][3], 2.97, 0.03) && Math.min(...w[0]) > 3,
193 w[2][3].toFixed(2))
194 check('space only: 23 cents at 24 against 146; 5.1 at 96', near(w[1][1], 23.3, 0.2) && near(w[0][1], 146, 1) && near(w[1][last], 5.10, 0.03),
195 `${w[1][1].toFixed(1)}, ${w[0][1].toFixed(0)}, ${w[1][last].toFixed(2)}`)
196 check('space only ends near second order', near(o.orders[1][last], 1.77, 0.05), o.orders[1][last].toFixed(2))
197 const nm = 3
198 check('at 96 the lowest mode is 0.7 cents sharp (space only), 0.1 flat (second order)',
199 near(o.cents[1 * nm][last], 0.69, 0.02) && near(o.cents[0][last], -0.095, 0.005), `${o.cents[nm][last].toFixed(3)}, ${o.cents[0][last].toFixed(3)}`)
200 check('fourth order in both costs about 4.5 times as much per step', o.ms[2][last] / o.ms[0][last] > 3 && o.ms[2][last] / o.ms[0][last] < 6.5,
201 (o.ms[2][last] / o.ms[0][last]).toFixed(2))
202 check('space only costs about 2.5 times the plain scheme at 96', o.ms[1][last] / o.ms[0][last] > 1.8 && o.ms[1][last] / o.ms[0][last] < 3.5,
203 (o.ms[1][last] / o.ms[0][last]).toFixed(2))
204 const st = (s, N) => runLab(MEM, S, { mode: 'strike', ...DRUM, N, lam: { fd2: 0.7, fd4: 0.6, fd44: 0.7 }[s], scheme: s }).out
205 const a = st('fd2', 32).rel
206 const b = st('fd4', 32).rel
207 const cc = st('fd44', 32).rel
208 check('N = 32: 0.69, 0.81, 0.37', near(a, 0.69, 0.006) && near(b, 0.81, 0.006) && near(cc, 0.37, 0.006), `${a.toFixed(3)} ${b.toFixed(3)} ${cc.toFixed(3)}`)
209 const d24 = st('fd44', 24)
210 const f48 = st('fd2', 48)
211 check('fourth order at 24 beats second order at 48, for less work', near(d24.rel, 0.46, 0.006) && near(f48.rel, 0.59, 0.006) && d24.ms < f48.ms,
212 `${d24.rel.toFixed(3)} in ${d24.ms.toFixed(0)} ms, ${f48.rel.toFixed(3)} in ${f48.ms.toFixed(0)} ms`)
213}
215// ---------------------------------------------------------------- Chapter 6
216if (section('ch6', 'Chapter 6: a round drum')) {
217 const D = src('ch06/driver.m')
218 const S = src('ch06/disk.m')
219 const { out: h } = runLab(D, S, {
220 mode: 'hear', fs: 22050, f01: 150, T: 1.5, width: 0.03, mmax: 40, kmax: 14, fmax: 4000, r0: 0.4, rr: 0.55, thr: 2.1, tau: 0.5,
221 })
222 const want = [1, 1.593, 2.136, 2.295, 2.653, 2.917]
223 check('the lowest ratios are 1, 1.593, 2.136, 2.295, 2.653, 2.917', want.every((v, i) => near(h.tratio[i], v, 0.0006)), h.tratio.slice(0, 6).map((v) => v.toFixed(4)).join(' '))
224 const { out: c } = runLab(D, S, { mode: 'hear', fs: 22050, f01: 150, T: 0.2, width: 0.03, mmax: 40, kmax: 14, fmax: 4000, r0: 0, rr: 0.55, thr: 2.1, tau: 0.5 })
225 check('struck at the centre, only m = 0 sounds', c.m.every((m) => m === 0), `${c.m.length} partials`)
226 check('at ratios 1, 2.295, 3.598, 4.903', [1, 2.295, 3.598, 4.903].every((v, i) => near(c.ratio[i], v, 0.001)), c.ratio.slice(0, 4).map((v) => v.toFixed(3)).join(' '))
228 const { out: g } = runLab(D, S, { mode: 'grid', Ns: [16, 24, 32, 40, 48, 64], methods: ['staircase', 'sw'], modes: [0, 1, 1, 1], periods: 10 })
229 const st = g.worst[0]
230 const sw = g.worst[1]
231 check('staircase: 54 cents at 16, 17 at 64', near(st[0], 53.9, 0.3) && near(st[5], 17.2, 0.2), `${st[0].toFixed(1)}, ${st[5].toFixed(1)}`)
232 check('staircase: fitted order about 0.8, and it rises from 24 to 32', near(Math.log(st[0] / st[5]) / Math.log(4), 0.82, 0.05) && st[2] > st[1],
233 (Math.log(st[0] / st[5]) / Math.log(4)).toFixed(2))
234 check('staircase errors are all flat', g.cents.slice(0, 2).flat().every((v) => v < 0))
235 check('Shortley-Weller: 7.4 at 16, 0.69 at 64, order 1.7', near(sw[0], 7.43, 0.05) && near(sw[5], 0.688, 0.01) && near(Math.log(sw[0] / sw[5]) / Math.log(4), 1.72, 0.03),
236 `${sw[0].toFixed(2)}, ${sw[5].toFixed(3)}`)
237 const ratio = st.map((v, i) => v / sw[i])
238 check('better by a factor of 7 to 25 on every grid', Math.min(...ratio) > 7 && Math.max(...ratio) < 26, `${Math.min(...ratio).toFixed(1)} to ${Math.max(...ratio).toFixed(1)}`)
239 check('below 5 cents by N = 24', sw[1] < 5)
240 check('its Courant number ranges 0.26 to 0.62; the staircase keeps 0.67',
241 near(Math.min(...g.lam[1]), 0.257, 0.005) && near(Math.max(...g.lam[1]), 0.615, 0.005) && g.lam[0].every((l) => near(l, 0.672, 0.001)))
242 check('Shortley-Weller at 16 beats the staircase at 64 at a small fraction of the time', sw[0] < st[5] && g.ms[1][0] < g.ms[0][5] / 15,
243 `${g.ms[1][0].toFixed(0)} ms against ${g.ms[0][5].toFixed(0)} ms`)
244}
246// ---------------------------------------------------------------- Chapter 7
247if (section('ch7', 'Chapter 7: letting sound out')) {
248 const D = src('ch07/driver.m')
249 const S = src('ch07/openbox.m')
250 const REGION = { N: 100, L: 2, lam: 0.7, sigma: 0.06, xs: 0.7, ys: 1.1, T: 0.0085 }
251 const { out: o } = runLab(D, S, { mode: 'sweep', ...REGION, methods: ['rigid', 'mur', 'sponge', 'pml'], Ws: [5, 10, 20, 40] })
252 const r = o.refl
253 check('rigid +2.4 dB, Mur -22.7 dB', near(r[0][0], 2.44, 0.05) && near(r[1][0], -22.74, 0.05), `${r[0][0].toFixed(2)}, ${r[1][0].toFixed(2)}`)
254 const sponge = [-4.1, -8.8, -14.8, -21.3]
255 const pml = [-62.1, -80.2, -95.2, -107.3]
256 check('damping layer -4.1, -8.8, -14.8, -21.3', sponge.every((v, j) => near(r[2][j], v, 0.06)), r[2].map((v) => v.toFixed(2)).join(' '))
257 check('PML -62.1, -80.2, -95.2, -107.3', pml.every((v, j) => near(r[3][j], v, 0.1)), r[3].map((v) => v.toFixed(2)).join(' '))
258 check('the widest damping layer does no better than Mur', r[2][3] > r[1][0])
259 check('a five-node PML is about 40 dB better than anything else', r[1][0] - r[3][0] > 38)
260 // a timing, so only its rough size is checked
261 check('a five-node PML costs about twice Mur, or a little more', o.ms[3][0] / o.ms[1][0] > 1.4 && o.ms[3][0] / o.ms[1][0] < 3.5, (o.ms[3][0] / o.ms[1][0]).toFixed(2))
262 // the snapshot of the default lab
263 const { out: s } = runLab(D, S, { mode: 'snap', ...REGION, tsnap: 0.0042, method: 'sponge', W: 10 })
264 check('the snapshot lab reproduces the sweep', near(s.refl, r[2][1], 1e-9), s.refl.toFixed(2))
265 // Mur's oblique reflection, as quoted
266 const mur = (a) => (1 - Math.cos(a)) / (1 + Math.cos(a))
267 check('Mur reflects 0.17 at 45 degrees and a third at 60', near(mur(Math.PI / 4), 0.17, 0.005) && near(mur(Math.PI / 3), 1 / 3, 1e-9))
268}
270// ---------------------------------------------------------------- Chapter 8
271if (section('ch8', 'Chapter 8: a room, by images')) {
272 const D = src('ch08/driver.m')
273 const S = src('ch08/images.m')
274 const ROOM = { L: [6, 4.5, 3], xs: [1.5, 1.2, 1.4], xr: [4.3, 3.1, 1.6] }
275 const { out: d } = runLab(D, S, { mode: 'decay', ...ROOM, fs: 16000, alphas: [0.05, 0.1, 0.2, 0.3, 0.5, 0.7], tcap: 1.5, taps: 8 })
276 const meas = [2.47, 1.28, 0.59, 0.37, 0.21, 0.12]
277 const lat = [2.55, 1.24, 0.59, 0.37, 0.19, 0.11]
278 const sab = [2.23, 1.12, 0.56, 0.37, 0.22, 0.16]
279 const eyr = [2.17, 1.06, 0.5, 0.31, 0.16, 0.09]
280 const two = (a, b) => a.every((v, i) => Math.abs(b[i] - v) <= 0.006 + 0.005 * v)
281 check('measured T60 as tabulated', two(meas, d.rt), d.rt.map((v) => v.toFixed(3)).join(' '))
282 check('lattice as tabulated', two(lat, d.lattice), d.lattice.map((v) => v.toFixed(3)).join(' '))
283 check('Sabine as tabulated', two(sab, d.sabine), d.sabine.map((v) => v.toFixed(3)).join(' '))
284 check('Eyring as tabulated', two(eyr, d.eyring), d.eyring.map((v) => v.toFixed(3)).join(' '))
285 const dev = d.rt.map((v, i) => v / d.lattice[i] - 1)
286 check('measured within about 3 % of the lattice up to 0.3', dev.slice(0, 4).every((v) => Math.abs(v) < 0.035), dev.map((v) => v.toFixed(3)).join(' '))
287 check('and 10 % and 14 % long at 0.5 and 0.7', near(dev[4], 0.1, 0.015) && near(dev[5], 0.14, 0.015))
288 const fac = d.lattice.map((v, i) => v / d.eyring[i])
289 check('the lattice is a fixed multiple of Eyring, 1.17', fac.every((v) => near(v, 1.175, 0.01)), fac.map((v) => v.toFixed(3)).join(' '))
290 const sl = d.sabine.map((v, i) => v / d.lattice[i])
291 check('Sabine short of the lattice by 13 % at 0.05 and long by 46 % at 0.7', near(sl[0], 0.87, 0.01) && near(sl[5], 1.46, 0.01), `${sl[0].toFixed(3)}, ${sl[5].toFixed(3)}`)
292 const room = (alpha, taps) => runLab(D, S, { mode: 'room', ...ROOM, fs: 22050, alpha, taps, tcap: 1.2, notes: [60, 64], step: 0.24 }).out
293 const r = room(0.2, 8)
294 check('the default room: 0.59 measured, 0.78 broadband, direct at 9.9 ms', near(r.rt, 0.59, 0.006) && near(r.rtbroad, 0.78, 0.006) && near(r.direct, 9.9, 0.05),
295 `${r.rt.toFixed(3)}, ${r.rtbroad.toFixed(3)}, ${r.direct.toFixed(2)}`)
296 const r0 = room(0.2, 0)
297 check('nearest-sample placement changes the decay by about one per cent', Math.abs(r0.rt / r.rt - 1) < 0.02, (r0.rt / r.rt - 1).toFixed(4))
298 check('alpha 0.1: about 1.3 s; alpha 0.6: about 0.15 s', near(room(0.1, 8).rt, 1.28, 0.03) && near(room(0.6, 8).rt, 0.155, 0.01))
299 const schroeder = 2000 * Math.sqrt(r.rt / 81)
300 check('the Schroeder frequency is about 170 Hz', near(schroeder, 170, 5), schroeder.toFixed(0))
301}
303// ---------------------------------------------------------------- Chapter 9
304if (section('ch9', 'Chapter 9: a room, by finite differences')) {
305 const D = src('ch09/driver.m')
306 const S = src('ch09/room3d.m')
307 const BOX = { L: [3.2, 2.4, 2.0], lam: 0.55, sigma: 0.0005, xs: [0.8, 0.6, 0.8], xr: [2.4, 1.6, 1.2] }
308 const { out: o } = runLab(D, S, { mode: 'converge', ...BOX, T: 0.03, Ns: [16, 24, 32, 48, 64] })
309 const want = [1.29, 0.92, 0.48, 0.21, 0.11]
310 check('errors 1.29, 0.92, 0.48, 0.21, 0.11', want.every((v, i) => near(o.err[i], v, 0.006)), o.err.map((v) => v.toFixed(3)).join(' '))
311 check('observed orders 0.84, 2.26, 2.06, 2.15', [0.84, 2.26, 2.06, 2.15].every((v, i) => near(o.orders[i + 1], v, 0.02)))
312 check('points per wavelength 1.5 to 6.0 at 1.1 kHz', near(o.fhi, 1136, 1) && near(o.ppw[0], 1.5, 0.02) && near(o.ppw[4], 6.0, 0.05))
313 const p1 = 6 * Math.sqrt(o.err[4] / 0.01)
314 check('about 20 points per wavelength for 1 per cent', near(p1, 20, 1), p1.toFixed(1))
315 check('about 1.6e10 nodes for the full band', rel(o.fullcells, 1.6e10, 0.02), o.fullcells.toExponential(2))
316 check('the time per node settles at a few tens of ns', o.rate > 10 && o.rate < 80, o.rate.toFixed(1))
317 const yearsAt38 = (o.fullband / o.rate) * 38 / 3.156e7
318 const hoursAt38 = (o.lowband / o.rate) * 38 / 3600
319 check('at 38 ns: about seven years, and about six hours to 2 kHz', near(yearsAt38, 7, 0.3) && near(hoursAt38, 6.1, 0.3), `${yearsAt38.toFixed(2)} y, ${hoursAt38.toFixed(2)} h`)
320 const { out: h } = runLab(D, S, { mode: 'hear', ...BOX, T: 0.3, fs: 22050, Ns: [16, 32] })
321 check('the listening grids run at 3.1 and 6.2 kHz', near(h.fsk[0], 3118, 1) && near(h.fsk[1], 6236, 1), h.fsk.map((v) => v.toFixed(0)).join(', '))
322 check('and their clips were raised to a playable rate', h.clips.every((c) => c.fs >= 8000), h.clips.map((c) => c.fs.toFixed(0)).join(', '))
323}
325console.log(failures === 0 ? '\nall checks passed' : `\n${failures} check(s) failed`)
326process.exit(failures === 0 ? 0 : 1)