nordic-skating / commit
feat(ice): why thin ice sings, as a model, and the film measured against it
src/ice/acoustics.ts | 112 ++++++++++++++++++++++++++++++++ src/world/lake.ts | 11 ++++ test/acoustics.test.ts | 57 +++++++++++++++++ test/lake.test.ts | 13 +++- tools/analyse-video.mjs | 166 ++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 358 insertions(+), 1 deletion(-) diff --git a/src/ice/acoustics.ts b/src/ice/acoustics.ts new file mode 100644 index 0000000..fd14673 --- /dev/null +++ b/src/ice/acoustics.ts @@ -0,0 +1,112 @@ +/** + * Why thin ice sings, as numbers. + * + * A sheet of ice floating on water carries bending ("flexural") waves, and they are + * dispersive: the higher the frequency, the faster it travels. Two things follow, and + * they are the two sounds of skating on black ice. + * + * - A crack or the strike of a skate is every frequency at once. Some way off, the high + * ones arrive first and the low ones after: a falling chirp, the "pew" of a laser. + * - Slow bending waves cannot radiate into the air at all. Sound leaves the ice only at + * frequencies whose bending wave travels at least as fast as sound does in air, and + * most strongly right at the frequency where the two speeds match, where the wave in + * the ice and the wave in the air keep step. That frequency is the note the ice sings, + * and it rises as the ice thins: the pitch is a thickness gauge. + * + * The model and its dispersion relation are from Romeyn, Hanssen, Ruud and Johansen, + * "Sea ice thickness from air-coupled flexural waves", The Cryosphere 15, 2939–2955 + * (2021), https://doi.org/10.5194/tc-15-2939-2021 — who measured this very lake from + * the film this app is modelled on: a note of about 725 Hz, from ice about 4.3 cm thick. + */ + +/** Young's modulus of freshwater ice, Pa. */ +export const YOUNGS_MODULUS = 9e9 +export const POISSON_RATIO = 0.33 +/** kg/m³ */ +export const ICE_DENSITY = 917 +export const WATER_DENSITY = 1000 +/** Speed of sound in air at -5 °C, m/s. */ +export const SOUND_SPEED_AIR = 328 +const GRAVITY = 9.81 + +/** How hard a sheet of ice `h` metres thick is to bend, N·m. */ +export function flexuralRigidity(h: number): number { + return (YOUNGS_MODULUS * h ** 3) / (12 * (1 - POISSON_RATIO ** 2)) +} + +/** + * Angular frequency, rad/s, of the bending wave of wavenumber `k` rad/m in ice `h` metres + * thick on deep water. The water matters: it has to be moved out of the way, which + * weighs the ice down and slows the wave. + */ +export function angularFrequency(k: number, h: number): number { + const stiffness = (flexuralRigidity(h) * k ** 5) / WATER_DENSITY + GRAVITY * k + return Math.sqrt(stiffness / (k * h * (ICE_DENSITY / WATER_DENSITY) + 1)) +} + +/** Speed of a single crest of that wave, m/s. */ +export function phaseSpeed(k: number, h: number): number { + return angularFrequency(k, h) / k +} + +/** Speed at which that wave carries a sound across the ice, m/s. */ +export function groupSpeed(k: number, h: number): number { + const dk = k * 1e-4 + return (angularFrequency(k + dk, h) - angularFrequency(k - dk, h)) / (2 * dk) +} + +/** Find where an increasing function crosses zero between `lo` and `hi`. */ +function solve(f: (x: number) => number, lo: number, hi: number): number { + for (let i = 0; i < 80; i++) { + const mid = (lo + hi) / 2 + if (f(mid) < 0) lo = mid + else hi = mid + } + return (lo + hi) / 2 +} + +/** Wavenumber, rad/m, of the bending wave of frequency `f` Hz in ice `h` metres thick. */ +export function wavenumber(f: number, h: number): number { + const omega = 2 * Math.PI * f + return solve((k) => angularFrequency(k, h) - omega, 1e-3, 1e4) +} + +/** + * The note the ice sings, Hz: the frequency whose bending wave moves at the speed of + * sound in air. + */ +export function singingFrequency(h: number): number { + // The wave gets faster as it gets shorter, so the match is where the speed crosses. + const k = solve((k) => phaseSpeed(k, h) - SOUND_SPEED_AIR, 0.5, 1e4) + return (SOUND_SPEED_AIR * k) / (2 * Math.PI) +} + +/** The other way round: how thick, in metres, is ice that sings at `f` Hz. */ +export function thicknessFromFrequency(f: number): number { + // Thicker ice sings lower. + return solve((h) => f - singingFrequency(h), 1e-3, 2) +} + +/** Seconds a sound of `f` Hz takes to cross `distance` metres of ice `h` metres thick. */ +export function travelTime(f: number, h: number, distance: number): number { + return distance / groupSpeed(wavenumber(f, h), h) +} + +/** + * The chirp heard `distance` metres from a crack in ice `h` metres thick: the frequency + * arriving at each moment, at `points` evenly spaced moments from the arrival of `from` + * Hz until the arrival of the ice's own note, which comes last and is where a chirp ends. + */ +export function chirp(h: number, distance: number, from = 6000, points = 64): { duration: number; frequencies: Float32Array } { + // Just above the note itself, where the fall has all but levelled out. + const lowest = singingFrequency(h) * 1.02 + const first = travelTime(from, h, distance) + const duration = travelTime(lowest, h, distance) - first + const frequencies = new Float32Array(points) + for (let i = 0; i < points; i++) { + const t = first + (duration * i) / (points - 1) + // Higher frequencies arrive sooner, so the arrival time falls as the frequency rises. + frequencies[i] = solve((f) => t - travelTime(f, h, distance), lowest, from) + } + return { duration, frequencies } +} diff --git a/src/world/lake.ts b/src/world/lake.ts index 1dcf51f..18e0058 100644 --- a/src/world/lake.ts +++ b/src/world/lake.ts @@ -27,11 +27,19 @@ export interface Lake { halfWidth(s: number): number /** Height of the land `d` metres inland of the left (-1) or right (+1) shore. */ bankHeight(side: -1 | 1, s: number, d: number): number + /** Thickness of the ice, metres, at `x` across (as `centre` is) and `s` along. */ + thickness(x: number, s: number): number } /** The narrowest and widest the channel gets, shore to middle, in metres. */ export const MIN_HALF_WIDTH = 12 export const MAX_HALF_WIDTH = 44 +/** + * The thinnest and thickest the ice gets, metres. The ice in the film this is modelled on + * measures 4.1 cm (tools/analyse-video.mjs), and that is the middle of the range. + */ +export const MIN_THICKNESS = 0.025 +export const MAX_THICKNESS = 0.057 export function createLake(seed: number): Lake { const noise = createNoise2D(mulberry32(seed)) @@ -50,5 +58,8 @@ export function createLake(seed: number): Lake { const rise = 8 * (1 - Math.exp(-d / 25)) + 0.06 * d return rise * (0.55 + 0.45 * noise(s / 140, side * 37 + d / 90)) }, + // Ice thickens and thins in broad patches, longer along the lake than across it. + thickness: (x, s) => + (MIN_THICKNESS + MAX_THICKNESS) / 2 + ((MAX_THICKNESS - MIN_THICKNESS) / 2) * noise(x / 60 + 500, s / 180), } } diff --git a/test/acoustics.test.ts b/test/acoustics.test.ts new file mode 100644 index 0000000..a2f4ddb --- /dev/null +++ b/test/acoustics.test.ts @@ -0,0 +1,57 @@ +import { describe, expect, it } from 'vitest' +import { + SOUND_SPEED_AIR, chirp, groupSpeed, phaseSpeed, singingFrequency, thicknessFromFrequency, + travelTime, wavenumber, +} from '../src/ice/acoustics.js' + +describe('ice acoustics', () => { + it('agrees with the ice measured in the paper the model is from', () => { + // Romeyn et al. 2021: the film's lake, about 725 Hz from about 4.3 cm; and a lake + // whose ice quakes had a 190 Hz component, about 16 cm. + expect(thicknessFromFrequency(725)).toBeCloseTo(0.043, 2) + expect(thicknessFromFrequency(190)).toBeCloseTo(0.16, 2) + }) + + it('sings higher the thinner it is, and very nearly in inverse proportion', () => { + expect(singingFrequency(0.03)).toBeGreaterThan(singingFrequency(0.05)) + for (const h of [0.02, 0.04, 0.08, 0.2]) { + expect(singingFrequency(h) * h).toBeCloseTo(30.3, 0) + expect(thicknessFromFrequency(singingFrequency(h))).toBeCloseTo(h, 6) + } + }) + + it('sings at the frequency whose bending wave keeps step with sound in air', () => { + const h = 0.041 + const note = singingFrequency(h) + expect(phaseSpeed(wavenumber(note, h), h)).toBeCloseTo(SOUND_SPEED_AIR, 3) + // Above the note the wave outruns sound in air and can radiate; below it, it cannot. + expect(phaseSpeed(wavenumber(note * 2, h), h)).toBeGreaterThan(SOUND_SPEED_AIR) + expect(phaseSpeed(wavenumber(note / 2, h), h)).toBeLessThan(SOUND_SPEED_AIR) + }) + + it('carries high frequencies faster than low ones', () => { + const h = 0.041 + let before = 0 + for (const f of [200, 740, 1500, 3000, 6000]) { + const speed = groupSpeed(wavenumber(f, h), h) + expect(speed).toBeGreaterThan(before) + before = speed + } + expect(travelTime(6000, h, 100)).toBeLessThan(travelTime(740, h, 100)) + }) + + it('makes a chirp that falls onto the note, and falls for longer from further off', () => { + const h = 0.041 + const near = chirp(h, 50) + const far = chirp(h, 400) + for (const { frequencies } of [near, far]) { + expect(frequencies[0]).toBeCloseTo(6000, -1) + for (let i = 1; i < frequencies.length; i++) expect(frequencies[i]).toBeLessThan(frequencies[i - 1]!) + expect(frequencies[frequencies.length - 1]! / singingFrequency(h)).toBeCloseTo(1.02, 2) + } + expect(far.duration / near.duration).toBeCloseTo(8, 5) + // A crack 400 m off falls for a third of a second: long enough to hear as a "pew". + expect(far.duration).toBeGreaterThan(0.25) + expect(far.duration).toBeLessThan(0.5) + }) +}) diff --git a/test/lake.test.ts b/test/lake.test.ts index 36ca5c6..dadd621 100644 --- a/test/lake.test.ts +++ b/test/lake.test.ts @@ -1,5 +1,5 @@ import { describe, expect, it } from 'vitest' -import { MAX_HALF_WIDTH, MIN_HALF_WIDTH, createLake } from '../src/world/lake.js' +import { MAX_HALF_WIDTH, MAX_THICKNESS, MIN_HALF_WIDTH, MIN_THICKNESS, createLake } from '../src/world/lake.js' describe('lake', () => { it('is the same lake for the same seed, and another for another', () => { @@ -39,4 +39,15 @@ describe('lake', () => { } } }) + + it('has ice between its thinnest and thickest everywhere', () => { + const lake = createLake(1) + for (let s = 0; s < 50_000; s += 13) { + for (const x of [-80, -20, 0, 35, 90]) { + const h = lake.thickness(x, s) + expect(h).toBeGreaterThanOrEqual(MIN_THICKNESS) + expect(h).toBeLessThanOrEqual(MAX_THICKNESS) + } + } + }) }) diff --git a/tools/analyse-video.mjs b/tools/analyse-video.mjs new file mode 100644 index 0000000..3afae95 --- /dev/null +++ b/tools/analyse-video.mjs @@ -0,0 +1,166 @@ +/** + * Measure the ice in the reference video: the note it sings, how often a skate strikes it, + * and how long it rings. + * + * node tools/analyse-video.mjs [video.mp4] + * + * The video is National Geographic's and is not part of this repository (README, + * Attribution): put a copy beside this project's package.json, or pass its path. Nothing + * from it is kept, only the numbers printed here. Needs ffmpeg. + */ +import { spawnSync } from 'node:child_process' +import { readdirSync } from 'node:fs' +import { thicknessFromFrequency } from '../src/ice/acoustics.ts' + +const RATE = 22050 +/** The film's skating ends here, seconds; after it there is music, which is also tonal. */ +const SKATING_ENDS = 170 + +const video = process.argv[2] ?? readdirSync('.').find((f) => f.endsWith('.mp4')) +if (!video) throw new Error('no video: pass its path, or put an .mp4 in this directory') +const decoded = spawnSync('ffmpeg', ['-v', 'error', '-i', video, '-ac', '1', '-ar', String(RATE), '-f', 'f32le', '-'], + { maxBuffer: 1 << 28 }) +if (decoded.status !== 0) throw new Error(`ffmpeg: ${decoded.stderr}`) +const audio = new Float32Array(decoded.stdout.buffer, decoded.stdout.byteOffset, decoded.stdout.byteLength >> 2) + +/** Power spectrum of `size` samples from `at`, Hann-windowed. `size` is a power of two. */ +function spectrum(at, size) { + const re = new Float64Array(size) + const im = new Float64Array(size) + for (let i = 0; i < size; i++) re[i] = (audio[at + i] ?? 0) * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / size)) + for (let i = 1, j = 0; i < size; i++) { + let bit = size >> 1 + for (; j & bit; bit >>= 1) j ^= bit + j ^= bit + if (i < j) [re[i], re[j]] = [re[j], re[i]] + } + for (let len = 2; len <= size; len <<= 1) { + const angle = (-2 * Math.PI) / len + for (let start = 0; start < size; start += len) { + for (let k = 0; k < len / 2; k++) { + const c = Math.cos(angle * k) + const s = Math.sin(angle * k) + const a = start + k + const b = a + len / 2 + const xr = re[b] * c - im[b] * s + const xi = re[b] * s + im[b] * c + re[b] = re[a] - xr + im[b] = im[a] - xi + re[a] += xr + im[a] += xi + } + } + } + const power = new Float64Array(size / 2) + for (let i = 0; i < size / 2; i++) power[i] = re[i] ** 2 + im[i] ** 2 + return power +} + +const quantile = (sorted, q) => sorted[Math.min(sorted.length - 1, Math.floor(q * sorted.length))] + +// --- The note ------------------------------------------------------------------------- +// Long windows, for a fine reading of pitch. A frame counts when one narrow peak between +// 300 and 1500 Hz stands well clear of the rest of that band. +const tones = [] +{ + const size = 4096 + const hz = RATE / size + const lo = Math.round(300 / hz) + const hi = Math.round(1500 / hz) + for (let at = 0; at + size < Math.min(audio.length, SKATING_ENDS * RATE); at += size / 4) { + const p = spectrum(at, size) + let peak = lo + for (let i = lo; i <= hi; i++) if (p[i] > p[peak]) peak = i + let rest = 0 + let n = 0 + for (let i = lo; i <= hi; i++) if (Math.abs(i - peak) > 6) { rest += p[i]; n++ } + if (p[peak] < 25 * (rest / n)) continue + // The true peak lies between bins; a parabola through the three around it finds it. + const [a, b, c] = [p[peak - 1], p[peak], p[peak + 1]].map(Math.log) + tones.push({ t: at / RATE, f: (peak + (0.5 * (a - c)) / (a - 2 * b + c)) * hz }) + } +} +const pitches = tones.map((x) => x.f).sort((a, b) => a - b) +const note = quantile(pitches, 0.5) +const thickness = thicknessFromFrequency(note) +console.log(`THE NOTE (${tones.length} clear frames of ${((tones.length * 1024) / RATE).toFixed(0)} s)`) +console.log(` median ${note.toFixed(0)} Hz, middle half ${quantile(pitches, 0.25).toFixed(0)}–${quantile(pitches, 0.75).toFixed(0)} Hz, ` + + `middle 90% ${quantile(pitches, 0.05).toFixed(0)}–${quantile(pitches, 0.95).toFixed(0)} Hz`) +console.log(` ice that sings at ${note.toFixed(0)} Hz is ${(thickness * 100).toFixed(1)} cm thick ` + + `(${(thicknessFromFrequency(quantile(pitches, 0.95)) * 100).toFixed(1)}–${(thicknessFromFrequency(quantile(pitches, 0.05)) * 100).toFixed(1)} cm over the middle 90%)`) +console.log(' by ten seconds of film:') +for (let t = 0; t < SKATING_ENDS; t += 10) { + const here = tones.filter((x) => x.t >= t && x.t < t + 10).map((x) => x.f).sort((a, b) => a - b) + if (here.length >= 8) console.log(` ${String(t).padStart(3)} s ${quantile(here, 0.5).toFixed(0)} Hz (${here.length} frames)`) +} + +// --- The strikes ---------------------------------------------------------------------- +// Short windows, for timing. Every push of a skate strikes the ice, and the ice rings at +// its note and dies away. Followed here: the loudness of the band around the note. +{ + const size = 512 + const hop = 128 + const hz = RATE / size + const perSecond = RATE / hop + const lo = Math.round((note * 0.9) / hz) + const hi = Math.round((note * 1.25) / hz) + const frames = Math.floor((Math.min(audio.length, SKATING_ENDS * RATE) - size) / hop) + const loud = new Float64Array(frames) + for (let n = 0; n < frames; n++) { + const p = spectrum(n * hop, size) + for (let i = lo; i <= hi; i++) loud[n] += p[i] + } + // A strike: the band is suddenly ten times louder than over the 150 ms before it. + const before = Math.round(0.15 * perSecond) + const strikes = [] + for (let n = before; n < frames; n++) { + let quiet = 0 + for (let i = n - before; i < n - 3; i++) quiet += loud[i] + quiet /= before - 3 + if (loud[n] > 10 * quiet && (strikes.length === 0 || n - strikes[strikes.length - 1] > 0.35 * perSecond)) strikes.push(n) + } + // Strikes that follow one another within three seconds are one run of skating. + const gaps = [] + for (let i = 1; i < strikes.length; i++) { + const gap = (strikes[i] - strikes[i - 1]) / perSecond + if (gap < 3) gaps.push(gap) + } + gaps.sort((a, b) => a - b) + // How fast the ring dies: a straight line through its loudness in decibels, from its + // loudest moment until the next strike or for a second, whichever is sooner. + const fades = [] + for (let i = 0; i < strikes.length; i++) { + const end = Math.min(strikes[i] + perSecond, (strikes[i + 1] ?? frames) - before, frames) + let top = strikes[i] + for (let n = strikes[i]; n < Math.min(end, strikes[i] + 0.1 * perSecond); n++) if (loud[n] > loud[top]) top = n + const ts = [] + const dbs = [] + for (let n = top; n < end; n++) { + ts.push((n - top) / perSecond) + dbs.push(10 * Math.log10(loud[n])) + } + if (ts.length < 0.3 * perSecond) continue + const mt = ts.reduce((a, b) => a + b) / ts.length + const md = dbs.reduce((a, b) => a + b) / dbs.length + let stt = 0 + let std = 0 + let sdd = 0 + for (let k = 0; k < ts.length; k++) { + stt += (ts[k] - mt) ** 2 + std += (ts[k] - mt) * (dbs[k] - md) + sdd += (dbs[k] - md) ** 2 + } + const slope = std / stt + if (slope < 0 && std ** 2 / (stt * sdd) > 0.5) fades.push(-slope) + } + fades.sort((a, b) => a - b) + console.log(`\nTHE STRIKES (${strikes.length} found)`) + console.log(` time between strikes: median ${quantile(gaps, 0.5).toFixed(2)} s, ` + + `middle half ${quantile(gaps, 0.25).toFixed(2)}–${quantile(gaps, 0.75).toFixed(2)} s (${gaps.length} gaps)`) + if (fades.length) { + const fade = quantile(fades, 0.5) + console.log(` the ring dies at a median ${fade.toFixed(0)} dB/s ` + + `(middle half ${quantile(fades, 0.25).toFixed(0)}–${quantile(fades, 0.75).toFixed(0)} dB/s, ${fades.length} rings): ` + + `its amplitude falls to 1/e in ${(8.686 / fade).toFixed(2)} s`) + } +}