nordic-skating / commit
feat(tools): test the film's chirps against the model, and find nothing to test
tools/analyse-video.mjs | 115 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 113 insertions(+), 2 deletions(-) diff --git a/tools/analyse-video.mjs b/tools/analyse-video.mjs index 3afae95..2c84140 100644 --- a/tools/analyse-video.mjs +++ b/tools/analyse-video.mjs @@ -1,6 +1,6 @@ /** * Measure the ice in the reference video: the note it sings, how often a skate strikes it, - * and how long it rings. + * how long it rings, and whether its chirps fall the way the model says they should. * * node tools/analyse-video.mjs [video.mp4] * @@ -10,7 +10,7 @@ */ import { spawnSync } from 'node:child_process' import { readdirSync } from 'node:fs' -import { thicknessFromFrequency } from '../src/ice/acoustics.ts' +import { groupSpeed, thicknessFromFrequency, wavenumber } from '../src/ice/acoustics.ts' const RATE = 22050 /** The film's skating ends here, seconds; after it there is music, which is also tonal. */ @@ -164,3 +164,114 @@ for (let t = 0; t < SKATING_ENDS; t += 10) { `its amplitude falls to 1/e in ${(8.686 / fade).toFixed(2)} s`) } } + +// --- The chirps ----------------------------------------------------------------------- +// A crack some way from the microphone arrives high frequencies first. The model says by +// how much: a frequency f arrives r / groupSpeed(f) after the crack, r metres away. On +// deep water under thin ice the group speed goes as f to the power 0.6, so the arrival +// time should fall as f to the power -0.6, whatever the distance and whatever the ice. +// That power is what is tested here, since it is the one thing about a chirp that does +// not depend on knowing where the crack was. +// +// The sound is split into bands from just above the note to 2.4 kHz, and for each sudden burst the moment +// each band comes up out of the quiet is its arrival. A burst counts as a chirp when most bands have a +// clear arrival, the arrivals are spread over more than a few hundredths of a second, +// and higher bands arrive earlier. +{ + const size = 512 + const hop = 16 + const hz = RATE / size + const perSecond = RATE / hop + const BANDS = 10 + // From just above the note, where nearly all of a chirp's sound is, up to 2.4 kHz. + const lowest = note * 1.08 + const edge = (b) => lowest * (2400 / lowest) ** (b / BANDS) + const middle = Array.from({ length: BANDS }, (_, b) => Math.sqrt(edge(b) * edge(b + 1))) + const frames = Math.floor((Math.min(audio.length, SKATING_ENDS * RATE) - size) / hop) + const loud = Array.from({ length: BANDS }, () => new Float32Array(frames)) + const total = new Float32Array(frames) + for (let n = 0; n < frames; n++) { + const p = spectrum(n * hop, size) + for (let b = 0; b < BANDS; b++) { + let sum = 0 + for (let i = Math.round(edge(b) / hz); i < Math.round(edge(b + 1) / hz); i++) sum += p[i] + loud[b][n] = sum + total[n] += sum + } + } + const before = Math.round(0.08 * perSecond) + const after = Math.round(0.4 * perSecond) + const lead = Math.round(0.015 * perSecond) + const mean = (a, from, to) => { let s = 0; for (let i = from; i < to; i++) s += a[i]; return s / (to - from) } + /** How well a straight line fits, 0 to 1, and its slope. */ + const line = (xs, ts) => { + const mx = xs.reduce((a, b) => a + b) / xs.length + const mt = ts.reduce((a, b) => a + b) / ts.length + let sxx = 0, sxt = 0, stt = 0 + for (let i = 0; i < xs.length; i++) { + sxx += (xs[i] - mx) ** 2 + sxt += (xs[i] - mx) * (ts[i] - mt) + stt += (ts[i] - mt) ** 2 + } + return { fit: sxt ** 2 / (sxx * stt), slope: sxt / sxx } + } + const chirps = [] + let bursts = 0 + let last = -Infinity + for (let n = before; n < frames - after; n++) { + if (total[n] < 20 * mean(total, n - before, n - 4) || n - last < 0.25 * perSecond) continue + last = n + bursts++ + const fs = [] + const ts = [] + for (let b = 0; b < BANDS; b++) { + const quiet = mean(loud[b], n - before, n - 4) + // The band's arrival is when it first comes up out of the quiet and stays up for + // a hundredth of a second: its onset, which later sounds in the band cannot move. + let onset = -1 + const stays = Math.round(0.01 * perSecond) + for (let i = n - lead; i < n + after && onset < 0; i++) { + if (loud[b][i] > 8 * quiet && mean(loud[b], i, i + stays) > 8 * quiet) onset = i + } + if (onset < 0) continue + fs.push(middle[b]) + ts.push((onset - n) / perSecond) + } + if (fs.length < 7) continue + const spread = Math.max(...ts) - Math.min(...ts) + if (spread < 0.03) continue + // Higher earlier: of every pair of bands, how many are in that order. + let ordered = 0 + let pairs = 0 + for (let i = 0; i < fs.length; i++) for (let j = i + 1; j < fs.length; j++) { pairs++; if (ts[j] < ts[i]) ordered++ } + if (ordered / pairs < 0.8) continue + // The power of f the arrivals fall as: whichever makes them lie straightest. + let best = { fit: 0, power: 0 } + for (let power = 0.1; power <= 2.001; power += 0.05) { + const { fit } = line(fs.map((f) => f ** -power), ts) + if (fit > best.fit) best = { fit, power } + } + const model = line(fs.map((f) => 1 / groupSpeed(wavenumber(f, thickness), thickness)), ts) + chirps.push({ t: (n * hop) / RATE, bands: fs.length, spread, power: best.power, powerFit: best.fit, r: model.slope, fit: model.fit }) + } + console.log(`\nTHE CHIRPS (${bursts} bursts; ${chirps.length} are chirps: most bands arriving, over 30 ms or more, higher first)`) + for (const c of chirps) { + console.log(` ${c.t.toFixed(2).padStart(7)} s ${String(c.bands).padStart(2)} bands over ${(c.spread * 1000).toFixed(0).padStart(3)} ms ` + + `falls as f^-${c.power.toFixed(2)} (fit ${c.powerFit.toFixed(2)}) by the model: ${c.r.toFixed(0).padStart(4)} m away (fit ${c.fit.toFixed(2)})`) + } + const good = chirps.filter((c) => c.fit > 0.8) + if (good.length < 5) { + // Which is how it comes out for the film. In nearly every burst every band arrives + // within 30 ms of every other: the microphones were close to the skater, and a sound + // from close by has had no distance to spread out over. The film has the note and + // not the chirp, and the chirp's law stands on the paper it is from (README). + console.log(` Too few to test anything with: ${bursts - chirps.length} of the ${bursts} bursts show no spread to measure, ` + + `and only ${good.length} of the rest fit the model well.`) + } else { + const powers = good.map((c) => c.power).sort((a, b) => a - b) + const far = good.map((c) => c.r).sort((a, b) => a - b) + console.log(` ${good.length} fit the model well (over 0.8). Among them the arrivals fall as f to the power ` + + `-${quantile(powers, 0.5).toFixed(2)} (middle half -${quantile(powers, 0.25).toFixed(2)} to -${quantile(powers, 0.75).toFixed(2)}); the model says -0.60.`) + console.log(` Their cracks were a median ${quantile(far, 0.5).toFixed(0)} m from the microphone (${far[0].toFixed(0)} to ${far[far.length - 1].toFixed(0)} m).`) + } +}