fuzzyworld

nordic-skating / mistral-medium-35-littlejs / files / tools / analyse-video.mjs

Branch mistral-medium-35-littlejs · default branch · all branches

analyse-video.mjs

Blame: the commit that last changed each line · file · history

174b4a1a 2026-10-09 Pavel Safronov1/**
2 * Measure the ice in the reference video: the note it sings, how often a skate strikes it,
1fdf083d 2026-10-10 Pavel Safronov3 * how long it rings, and whether its chirps fall the way the model says they should.
174b4a1a 2026-10-09 Pavel Safronov4 *
5 * node tools/analyse-video.mjs [video.mp4]
6 *
7 * The video is National Geographic's and is not part of this repository (README,
8 * Attribution): put a copy beside this project's package.json, or pass its path. Nothing
9 * from it is kept, only the numbers printed here. Needs ffmpeg.
10 */
11import { spawnSync } from 'node:child_process'
12import { readdirSync } from 'node:fs'
1fdf083d 2026-10-10 Pavel Safronov13import { groupSpeed, thicknessFromFrequency, wavenumber } from '../src/ice/acoustics.ts'
174b4a1a 2026-10-09 Pavel Safronov14
15const RATE = 22050
16/** The film's skating ends here, seconds; after it there is music, which is also tonal. */
17const SKATING_ENDS = 170
18
19const video = process.argv[2] ?? readdirSync('.').find((f) => f.endsWith('.mp4'))
20if (!video) throw new Error('no video: pass its path, or put an .mp4 in this directory')
21const decoded = spawnSync('ffmpeg', ['-v', 'error', '-i', video, '-ac', '1', '-ar', String(RATE), '-f', 'f32le', '-'],
22 { maxBuffer: 1 << 28 })
23if (decoded.status !== 0) throw new Error(`ffmpeg: ${decoded.stderr}`)
24const audio = new Float32Array(decoded.stdout.buffer, decoded.stdout.byteOffset, decoded.stdout.byteLength >> 2)
25
26/** Power spectrum of `size` samples from `at`, Hann-windowed. `size` is a power of two. */
27function spectrum(at, size) {
28 const re = new Float64Array(size)
29 const im = new Float64Array(size)
30 for (let i = 0; i < size; i++) re[i] = (audio[at + i] ?? 0) * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / size))
31 for (let i = 1, j = 0; i < size; i++) {
32 let bit = size >> 1
33 for (; j & bit; bit >>= 1) j ^= bit
34 j ^= bit
35 if (i < j) [re[i], re[j]] = [re[j], re[i]]
36 }
37 for (let len = 2; len <= size; len <<= 1) {
38 const angle = (-2 * Math.PI) / len
39 for (let start = 0; start < size; start += len) {
40 for (let k = 0; k < len / 2; k++) {
41 const c = Math.cos(angle * k)
42 const s = Math.sin(angle * k)
43 const a = start + k
44 const b = a + len / 2
45 const xr = re[b] * c - im[b] * s
46 const xi = re[b] * s + im[b] * c
47 re[b] = re[a] - xr
48 im[b] = im[a] - xi
49 re[a] += xr
50 im[a] += xi
51 }
52 }
53 }
54 const power = new Float64Array(size / 2)
55 for (let i = 0; i < size / 2; i++) power[i] = re[i] ** 2 + im[i] ** 2
56 return power
57}
58
59const quantile = (sorted, q) => sorted[Math.min(sorted.length - 1, Math.floor(q * sorted.length))]
60
61// --- The note -------------------------------------------------------------------------
62// Long windows, for a fine reading of pitch. A frame counts when one narrow peak between
63// 300 and 1500 Hz stands well clear of the rest of that band.
64const tones = []
65{
66 const size = 4096
67 const hz = RATE / size
68 const lo = Math.round(300 / hz)
69 const hi = Math.round(1500 / hz)
70 for (let at = 0; at + size < Math.min(audio.length, SKATING_ENDS * RATE); at += size / 4) {
71 const p = spectrum(at, size)
72 let peak = lo
73 for (let i = lo; i <= hi; i++) if (p[i] > p[peak]) peak = i
74 let rest = 0
75 let n = 0
76 for (let i = lo; i <= hi; i++) if (Math.abs(i - peak) > 6) { rest += p[i]; n++ }
77 if (p[peak] < 25 * (rest / n)) continue
78 // The true peak lies between bins; a parabola through the three around it finds it.
79 const [a, b, c] = [p[peak - 1], p[peak], p[peak + 1]].map(Math.log)
80 tones.push({ t: at / RATE, f: (peak + (0.5 * (a - c)) / (a - 2 * b + c)) * hz })
81 }
82}
83const pitches = tones.map((x) => x.f).sort((a, b) => a - b)
84const note = quantile(pitches, 0.5)
85const thickness = thicknessFromFrequency(note)
86console.log(`THE NOTE (${tones.length} clear frames of ${((tones.length * 1024) / RATE).toFixed(0)} s)`)
87console.log(` median ${note.toFixed(0)} Hz, middle half ${quantile(pitches, 0.25).toFixed(0)}–${quantile(pitches, 0.75).toFixed(0)} Hz, ` +
88 `middle 90% ${quantile(pitches, 0.05).toFixed(0)}–${quantile(pitches, 0.95).toFixed(0)} Hz`)
89console.log(` ice that sings at ${note.toFixed(0)} Hz is ${(thickness * 100).toFixed(1)} cm thick ` +
90 `(${(thicknessFromFrequency(quantile(pitches, 0.95)) * 100).toFixed(1)}–${(thicknessFromFrequency(quantile(pitches, 0.05)) * 100).toFixed(1)} cm over the middle 90%)`)
91console.log(' by ten seconds of film:')
92for (let t = 0; t < SKATING_ENDS; t += 10) {
93 const here = tones.filter((x) => x.t >= t && x.t < t + 10).map((x) => x.f).sort((a, b) => a - b)
94 if (here.length >= 8) console.log(` ${String(t).padStart(3)} s ${quantile(here, 0.5).toFixed(0)} Hz (${here.length} frames)`)
95}
96
97// --- The strikes ----------------------------------------------------------------------
98// Short windows, for timing. Every push of a skate strikes the ice, and the ice rings at
99// its note and dies away. Followed here: the loudness of the band around the note.
100{
101 const size = 512
102 const hop = 128
103 const hz = RATE / size
104 const perSecond = RATE / hop
105 const lo = Math.round((note * 0.9) / hz)
106 const hi = Math.round((note * 1.25) / hz)
107 const frames = Math.floor((Math.min(audio.length, SKATING_ENDS * RATE) - size) / hop)
108 const loud = new Float64Array(frames)
109 for (let n = 0; n < frames; n++) {
110 const p = spectrum(n * hop, size)
111 for (let i = lo; i <= hi; i++) loud[n] += p[i]
112 }
113 // A strike: the band is suddenly ten times louder than over the 150 ms before it.
114 const before = Math.round(0.15 * perSecond)
115 const strikes = []
116 for (let n = before; n < frames; n++) {
117 let quiet = 0
118 for (let i = n - before; i < n - 3; i++) quiet += loud[i]
119 quiet /= before - 3
120 if (loud[n] > 10 * quiet && (strikes.length === 0 || n - strikes[strikes.length - 1] > 0.35 * perSecond)) strikes.push(n)
121 }
122 // Strikes that follow one another within three seconds are one run of skating.
123 const gaps = []
124 for (let i = 1; i < strikes.length; i++) {
125 const gap = (strikes[i] - strikes[i - 1]) / perSecond
126 if (gap < 3) gaps.push(gap)
127 }
128 gaps.sort((a, b) => a - b)
129 // How fast the ring dies: a straight line through its loudness in decibels, from its
130 // loudest moment until the next strike or for a second, whichever is sooner.
131 const fades = []
132 for (let i = 0; i < strikes.length; i++) {
133 const end = Math.min(strikes[i] + perSecond, (strikes[i + 1] ?? frames) - before, frames)
134 let top = strikes[i]
135 for (let n = strikes[i]; n < Math.min(end, strikes[i] + 0.1 * perSecond); n++) if (loud[n] > loud[top]) top = n
136 const ts = []
137 const dbs = []
138 for (let n = top; n < end; n++) {
139 ts.push((n - top) / perSecond)
140 dbs.push(10 * Math.log10(loud[n]))
141 }
142 if (ts.length < 0.3 * perSecond) continue
143 const mt = ts.reduce((a, b) => a + b) / ts.length
144 const md = dbs.reduce((a, b) => a + b) / dbs.length
145 let stt = 0
146 let std = 0
147 let sdd = 0
148 for (let k = 0; k < ts.length; k++) {
149 stt += (ts[k] - mt) ** 2
150 std += (ts[k] - mt) * (dbs[k] - md)
151 sdd += (dbs[k] - md) ** 2
152 }
153 const slope = std / stt
154 if (slope < 0 && std ** 2 / (stt * sdd) > 0.5) fades.push(-slope)
155 }
156 fades.sort((a, b) => a - b)
157 console.log(`\nTHE STRIKES (${strikes.length} found)`)
158 console.log(` time between strikes: median ${quantile(gaps, 0.5).toFixed(2)} s, ` +
159 `middle half ${quantile(gaps, 0.25).toFixed(2)}–${quantile(gaps, 0.75).toFixed(2)} s (${gaps.length} gaps)`)
160 if (fades.length) {
161 const fade = quantile(fades, 0.5)
162 console.log(` the ring dies at a median ${fade.toFixed(0)} dB/s ` +
163 `(middle half ${quantile(fades, 0.25).toFixed(0)}–${quantile(fades, 0.75).toFixed(0)} dB/s, ${fades.length} rings): ` +
164 `its amplitude falls to 1/e in ${(8.686 / fade).toFixed(2)} s`)
165 }
166}
1fdf083d 2026-10-10 Pavel Safronov167
168// --- The chirps -----------------------------------------------------------------------
169// A crack some way from the microphone arrives high frequencies first. The model says by
170// how much: a frequency f arrives r / groupSpeed(f) after the crack, r metres away. On
171// deep water under thin ice the group speed goes as f to the power 0.6, so the arrival
172// time should fall as f to the power -0.6, whatever the distance and whatever the ice.
173// That power is what is tested here, since it is the one thing about a chirp that does
174// not depend on knowing where the crack was.
175//
176// The sound is split into bands from just above the note to 2.4 kHz, and for each sudden burst the moment
177// each band comes up out of the quiet is its arrival. A burst counts as a chirp when most bands have a
178// clear arrival, the arrivals are spread over more than a few hundredths of a second,
179// and higher bands arrive earlier.
180{
181 const size = 512
182 const hop = 16
183 const hz = RATE / size
184 const perSecond = RATE / hop
185 const BANDS = 10
186 // From just above the note, where nearly all of a chirp's sound is, up to 2.4 kHz.
187 const lowest = note * 1.08
188 const edge = (b) => lowest * (2400 / lowest) ** (b / BANDS)
189 const middle = Array.from({ length: BANDS }, (_, b) => Math.sqrt(edge(b) * edge(b + 1)))
190 const frames = Math.floor((Math.min(audio.length, SKATING_ENDS * RATE) - size) / hop)
191 const loud = Array.from({ length: BANDS }, () => new Float32Array(frames))
192 const total = new Float32Array(frames)
193 for (let n = 0; n < frames; n++) {
194 const p = spectrum(n * hop, size)
195 for (let b = 0; b < BANDS; b++) {
196 let sum = 0
197 for (let i = Math.round(edge(b) / hz); i < Math.round(edge(b + 1) / hz); i++) sum += p[i]
198 loud[b][n] = sum
199 total[n] += sum
200 }
201 }
202 const before = Math.round(0.08 * perSecond)
203 const after = Math.round(0.4 * perSecond)
204 const lead = Math.round(0.015 * perSecond)
205 const mean = (a, from, to) => { let s = 0; for (let i = from; i < to; i++) s += a[i]; return s / (to - from) }
206 /** How well a straight line fits, 0 to 1, and its slope. */
207 const line = (xs, ts) => {
208 const mx = xs.reduce((a, b) => a + b) / xs.length
209 const mt = ts.reduce((a, b) => a + b) / ts.length
210 let sxx = 0, sxt = 0, stt = 0
211 for (let i = 0; i < xs.length; i++) {
212 sxx += (xs[i] - mx) ** 2
213 sxt += (xs[i] - mx) * (ts[i] - mt)
214 stt += (ts[i] - mt) ** 2
215 }
216 return { fit: sxt ** 2 / (sxx * stt), slope: sxt / sxx }
217 }
218 const chirps = []
219 let bursts = 0
220 let last = -Infinity
221 for (let n = before; n < frames - after; n++) {
222 if (total[n] < 20 * mean(total, n - before, n - 4) || n - last < 0.25 * perSecond) continue
223 last = n
224 bursts++
225 const fs = []
226 const ts = []
227 for (let b = 0; b < BANDS; b++) {
228 const quiet = mean(loud[b], n - before, n - 4)
229 // The band's arrival is when it first comes up out of the quiet and stays up for
230 // a hundredth of a second: its onset, which later sounds in the band cannot move.
231 let onset = -1
232 const stays = Math.round(0.01 * perSecond)
233 for (let i = n - lead; i < n + after && onset < 0; i++) {
234 if (loud[b][i] > 8 * quiet && mean(loud[b], i, i + stays) > 8 * quiet) onset = i
235 }
236 if (onset < 0) continue
237 fs.push(middle[b])
238 ts.push((onset - n) / perSecond)
239 }
240 if (fs.length < 7) continue
241 const spread = Math.max(...ts) - Math.min(...ts)
242 if (spread < 0.03) continue
243 // Higher earlier: of every pair of bands, how many are in that order.
244 let ordered = 0
245 let pairs = 0
246 for (let i = 0; i < fs.length; i++) for (let j = i + 1; j < fs.length; j++) { pairs++; if (ts[j] < ts[i]) ordered++ }
247 if (ordered / pairs < 0.8) continue
248 // The power of f the arrivals fall as: whichever makes them lie straightest.
249 let best = { fit: 0, power: 0 }
250 for (let power = 0.1; power <= 2.001; power += 0.05) {
251 const { fit } = line(fs.map((f) => f ** -power), ts)
252 if (fit > best.fit) best = { fit, power }
253 }
254 const model = line(fs.map((f) => 1 / groupSpeed(wavenumber(f, thickness), thickness)), ts)
255 chirps.push({ t: (n * hop) / RATE, bands: fs.length, spread, power: best.power, powerFit: best.fit, r: model.slope, fit: model.fit })
256 }
257 console.log(`\nTHE CHIRPS (${bursts} bursts; ${chirps.length} are chirps: most bands arriving, over 30 ms or more, higher first)`)
258 for (const c of chirps) {
259 console.log(` ${c.t.toFixed(2).padStart(7)} s ${String(c.bands).padStart(2)} bands over ${(c.spread * 1000).toFixed(0).padStart(3)} ms ` +
260 `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)})`)
261 }
262 const good = chirps.filter((c) => c.fit > 0.8)
263 if (good.length < 5) {
264 // Which is how it comes out for the film. In nearly every burst every band arrives
265 // within 30 ms of every other: the microphones were close to the skater, and a sound
266 // from close by has had no distance to spread out over. The film has the note and
267 // not the chirp, and the chirp's law stands on the paper it is from (README).
268 console.log(` Too few to test anything with: ${bursts - chirps.length} of the ${bursts} bursts show no spread to measure, ` +
269 `and only ${good.length} of the rest fit the model well.`)
270 } else {
271 const powers = good.map((c) => c.power).sort((a, b) => a - b)
272 const far = good.map((c) => c.r).sort((a, b) => a - b)
273 console.log(` ${good.length} fit the model well (over 0.8). Among them the arrivals fall as f to the power ` +
274 `-${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.`)
275 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).`)
276 }
277}