analyse-video.mjs

1	/**
2	* Measure the ice in the reference video: the note it sings, how often a skate strikes it,
3	* how long it rings, and whether its chirps fall the way the model says they should.
4	*
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	*/
11	import { spawnSync } from 'node:child_process'
12	import { readdirSync } from 'node:fs'
13	import { groupSpeed, thicknessFromFrequency, wavenumber } from '../src/ice/acoustics.ts'
14
15	const RATE = 22050
16	/** The film's skating ends here, seconds; after it there is music, which is also tonal. */
17	const SKATING_ENDS = 170
18
19	const video = process.argv[2] ?? readdirSync('.').find((f) => f.endsWith('.mp4'))
20	if (!video) throw new Error('no video: pass its path, or put an .mp4 in this directory')
21	const decoded = spawnSync('ffmpeg', ['-v', 'error', '-i', video, '-ac', '1', '-ar', String(RATE), '-f', 'f32le', '-'],
22	{ maxBuffer: 1 << 28 })
23	if (decoded.status !== 0) throw new Error(`ffmpeg: ${decoded.stderr}`)
24	const 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. */
27	function 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
59	const 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.
64	const 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	}
83	const pitches = tones.map((x) => x.f).sort((a, b) => a - b)
84	const note = quantile(pitches, 0.5)
85	const thickness = thicknessFromFrequency(note)
86	console.log(`THE NOTE (${tones.length} clear frames of ${((tones.length * 1024) / RATE).toFixed(0)} s)`)
87	console.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`)
89	console.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%)`)
91	console.log(' by ten seconds of film:')
92	for (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	}
167
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	}