7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 1/**
2 * Phase maps: the angle phi_k applied to band k, as a function of the band
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 3 * center frequency and, for the random map, optionally of time. Maps must
4 * vary smoothly from one band to the next, otherwise neighbouring bands
5 * partly cancel where they overlap and the magnitude response acquires
6 * notches. A map that varies in time is frequency modulation of each band:
7 * the instantaneous frequency shifts by (d phi / dt) / 2 pi.
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 8 */
10export type PhaseMapSpec =
11 | {
12 kind: 'random'
13 /** Gaussian smoothing width across bands */
14 smoothBands: number
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 15 /** Gaussian smoothing width in time; null means the map does not change in time */
16 timeScaleMs: number | null
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 17 seed: number
18 /** multiplies the map; 1 gives the spread of uniform random phases */
19 strength: number
20 }
21 | {
22 kind: 'delay'
23 /** phi = -2 pi f tau: a pure delay, which must be inaudible */
24 delayMs: number
25 }
26 | {
27 kind: 'chirp'
28 /** group delay grows linearly with frequency, reaching this value at the Nyquist frequency */
29 delayAtNyquistMs: number
30 }
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 32export const defaultPhaseMap: PhaseMapSpec = { kind: 'random', smoothBands: 8, timeScaleMs: null, seed: 1, strength: 1 }
34export const isTimeVarying = (spec: PhaseMapSpec): boolean => spec.kind === 'random' && spec.timeScaleMs !== null
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 35
36/** Small seeded generator (mulberry32) so a map can be reproduced from its seed. */
37export function mulberry32(seed: number): () => number {
38 let a = seed >>> 0
39 return () => {
40 a = (a + 0x6d2b79f5) >>> 0
41 let t = a
42 t = Math.imul(t ^ (t >>> 15), t | 1)
43 t ^= t + Math.imul(t ^ (t >>> 7), t | 61)
44 return ((t ^ (t >>> 14)) >>> 0) / 4294967296
45 }
46}
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 48/** Standard deviation of a phase drawn uniformly from (-pi, pi]. */
49export const UNIFORM_STD = Math.PI / Math.sqrt(3)
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 50
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 51/**
52 * A random phase field on a coarse grid of (band, time) cells, smoothed with
53 * Gaussians in both directions and interpolated bilinearly. Row r holds band
54 * index r * dK and column c holds time c * dT.
55 */
56export type PhaseField = {
57 values: Float32Array
58 nK: number
59 nT: number
60 dK: number
61 dT: number
62 /** largest |d phi / dt| / 2 pi over the grid, in Hz: the largest frequency shift the field imposes */
63 maxRateHz: number
64}
66/** Gaussian smoothing along one axis of a row-major nK x nT grid, with the kernel renormalized at the edges. */
67function smoothAxis(v: Float32Array, nK: number, nT: number, alongTime: boolean, sigmaCells: number): void {
68 if (sigmaCells < 0.3) return
69 const radius = Math.ceil(4 * sigmaCells)
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 70 const kern = new Float64Array(2 * radius + 1)
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 71 for (let j = -radius; j <= radius; j++) kern[j + radius] = Math.exp(-0.5 * (j / sigmaCells) ** 2)
72 const len = alongTime ? nT : nK
73 const lines = alongTime ? nK : nT
74 const stride = alongTime ? 1 : nT
75 const lineStride = alongTime ? nT : 1
76 const tmp = new Float64Array(len)
77 for (let l = 0; l < lines; l++) {
78 const base = l * lineStride
79 for (let i = 0; i < len; i++) {
80 let s = 0
81 let wsum = 0
82 const j0 = Math.max(-radius, -i)
83 const j1 = Math.min(radius, len - 1 - i)
84 for (let j = j0; j <= j1; j++) {
85 const w = kern[j + radius]
86 s += w * v[base + (i + j) * stride]
87 wsum += w
88 }
89 tmp[i] = s / wsum
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 90 }
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 91 for (let i = 0; i < len; i++) v[base + i * stride] = tmp[i]
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 92 }
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 93}
95export function randomField(nBands: number, durationSec: number, smoothBands: number, timeScaleSec: number | null, seed: number, strength: number): PhaseField {
96 const sb = Math.max(0, smoothBands)
97 // grid spacing: about a third of each smoothing scale, bounded so the grid stays small
98 const dK = Math.max(1, Math.floor(sb / 3), Math.ceil(nBands / 1024))
99 const nK = Math.ceil(Math.max(0, nBands - 1) / dK) + 1
100 let nT = 1
101 let dT = durationSec
102 if (timeScaleSec !== null) {
103 dT = Math.max(timeScaleSec / 3, durationSec / 4096)
104 nT = Math.ceil(durationSec / dT) + 1
105 }
106 const v = new Float32Array(nK * nT)
107 const rnd = mulberry32(seed)
108 for (let i = 0; i < v.length; i++) v[i] = (rnd() * 2 - 1) * Math.PI
109 smoothAxis(v, nK, nT, false, sb / dK)
110 if (timeScaleSec !== null) smoothAxis(v, nK, nT, true, timeScaleSec / dT)
111 // rescale so the field keeps the spread of uniform phases, times the strength
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 112 let mean = 0
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 113 for (let i = 0; i < v.length; i++) mean += v[i]
114 mean /= v.length
115 let s2 = 0
116 for (let i = 0; i < v.length; i++) s2 += (v[i] - mean) ** 2
117 const std = Math.sqrt(s2 / Math.max(1, v.length - 1))
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 118 const scale = std > 0 ? (strength * UNIFORM_STD) / std : 0
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 119 for (let i = 0; i < v.length; i++) v[i] = (v[i] - mean) * scale
120 let maxRate = 0
121 if (nT > 1) {
122 for (let r = 0; r < nK; r++) {
123 for (let c = 0; c + 1 < nT; c++) {
124 const d = Math.abs(v[r * nT + c + 1] - v[r * nT + c])
125 if (d > maxRate) maxRate = d
126 }
127 }
128 maxRate /= 2 * Math.PI * dT
129 }
130 return { values: v, nK, nT, dK, dT, maxRateHz: maxRate }
131}
133/** The field for band k as a function of the grid time index (interpolated between grid rows). */
134export function fieldRow(field: PhaseField, k: number): Float64Array {
135 const { values, nK, nT, dK } = field
136 const r = Math.min(nK - 1, k / dK)
137 const r0 = Math.floor(r)
138 const r1 = Math.min(nK - 1, r0 + 1)
139 const w = r - r0
140 const out = new Float64Array(nT)
141 for (let c = 0; c < nT; c++) out[c] = (1 - w) * values[r0 * nT + c] + w * values[r1 * nT + c]
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 142 return out
143}
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 145/** Linear interpolation of a grid row at arbitrary times (seconds). */
146export function sampleRow(row: Float64Array, dT: number, times: ArrayLike<number>, out: Float64Array): Float64Array {
147 const nT = row.length
148 for (let i = 0; i < times.length; i++) {
149 if (nT === 1) {
150 out[i] = row[0]
151 continue
152 }
153 const c = Math.min(nT - 1, Math.max(0, times[i] / dT))
154 const c0 = Math.floor(c)
155 const c1 = Math.min(nT - 1, c0 + 1)
156 const w = c - c0
157 out[i] = (1 - w) * row[c0] + w * row[c1]
158 }
159 return out
160}
162export function fieldAt(field: PhaseField, k: number, tSec: number): number {
163 const row = fieldRow(field, k)
164 return sampleRow(row, field.dT, [tSec], new Float64Array(1))[0]
165}
167/** phi_k for every band center, for maps that do not vary in time. */
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 168export function phaseMap(spec: PhaseMapSpec, centers: Float64Array, fs: number): Float64Array {
169 const n = centers.length
170 switch (spec.kind) {
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 171 case 'random': {
172 const field = randomField(n, 1, spec.smoothBands, null, spec.seed, spec.strength)
173 const out = new Float64Array(n)
174 for (let k = 0; k < n; k++) out[k] = fieldAt(field, k, 0)
175 return out
176 }
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 177 case 'delay': {
178 const tau = spec.delayMs / 1000
179 return centers.map((f) => -2 * Math.PI * f * tau)
180 }
181 case 'chirp': {
182 // tau(f) = tauN f / fN => phi(f) = -2 pi int tau df = -pi tauN f^2 / fN
183 const tauN = spec.delayAtNyquistMs / 1000
184 const fN = fs / 2
185 return centers.map((f) => (-Math.PI * tauN * f * f) / fN)
186 }
187 }
188}
190/** Largest band-to-band slope of the map, as a group delay in seconds. Used to size the FFT padding. */
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 191export function maxGroupDelaySec(phi: ArrayLike<number>, df: number): number {
7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 192 let m = 0
193 for (let i = 1; i < phi.length; i++) {
194 let d = Math.abs(phi[i] - phi[i - 1])
195 d = Math.min(d, 2 * Math.PI) // a jump larger than 2 pi wraps to something smaller anyway
196 if (d > m) m = d
197 }
198 return m / (2 * Math.PI * df)
199}
201export const wrapAngle = (x: number): number => {
202 let y = x % (2 * Math.PI)
203 if (y > Math.PI) y -= 2 * Math.PI
204 if (y <= -Math.PI) y += 2 * Math.PI
205 return y
206}