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
3 * center frequency. The maps must vary smoothly from one band to the next,
4 * otherwise neighbouring bands partly cancel where they overlap and the
5 * magnitude response acquires notches.
6 */
8export type PhaseMapSpec =
9 | {
10 kind: 'random'
11 /** Gaussian smoothing width across bands */
12 smoothBands: number
13 seed: number
14 /** multiplies the map; 1 gives the spread of uniform random phases */
15 strength: number
16 }
17 | {
18 kind: 'delay'
19 /** phi = -2 pi f tau: a pure delay, which must be inaudible */
20 delayMs: number
21 }
22 | {
23 kind: 'chirp'
24 /** group delay grows linearly with frequency, reaching this value at the Nyquist frequency */
25 delayAtNyquistMs: number
26 }
28export const defaultPhaseMap: PhaseMapSpec = { kind: 'random', smoothBands: 8, seed: 1, strength: 1 }
30/** Small seeded generator (mulberry32) so a map can be reproduced from its seed. */
31export function mulberry32(seed: number): () => number {
32 let a = seed >>> 0
33 return () => {
34 a = (a + 0x6d2b79f5) >>> 0
35 let t = a
36 t = Math.imul(t ^ (t >>> 15), t | 1)
37 t ^= t + Math.imul(t ^ (t >>> 7), t | 61)
38 return ((t ^ (t >>> 14)) >>> 0) / 4294967296
39 }
40}
42const UNIFORM_STD = Math.PI / Math.sqrt(3)
44function smoothedRandom(n: number, smoothBands: number, seed: number, strength: number): Float64Array {
45 const rnd = mulberry32(seed)
46 const raw = new Float64Array(n)
47 for (let i = 0; i < n; i++) raw[i] = (rnd() * 2 - 1) * Math.PI
48 if (smoothBands <= 0) return raw.map((v) => v * strength)
49 const radius = Math.ceil(4 * smoothBands)
50 const kern = new Float64Array(2 * radius + 1)
51 for (let j = -radius; j <= radius; j++) kern[j + radius] = Math.exp(-0.5 * (j / smoothBands) ** 2)
52 const out = new Float64Array(n)
53 for (let i = 0; i < n; i++) {
54 let s = 0
55 let wsum = 0
56 for (let j = -radius; j <= radius; j++) {
57 const m = i + j
58 if (m < 0 || m >= n) continue
59 const w = kern[j + radius]
60 s += w * raw[m]
61 wsum += w
62 }
63 out[i] = s / wsum
64 }
65 // Smoothing shrinks the spread; rescale so the map keeps the std of uniform phases.
66 let mean = 0
67 for (let i = 0; i < n; i++) mean += out[i]
68 mean /= n
69 let v = 0
70 for (let i = 0; i < n; i++) v += (out[i] - mean) ** 2
71 const std = Math.sqrt(v / Math.max(1, n - 1))
72 const scale = std > 0 ? (strength * UNIFORM_STD) / std : 0
73 for (let i = 0; i < n; i++) out[i] = (out[i] - mean) * scale
74 return out
75}
77/** phi_k for every band center. */
78export function phaseMap(spec: PhaseMapSpec, centers: Float64Array, fs: number): Float64Array {
79 const n = centers.length
80 switch (spec.kind) {
81 case 'random':
82 return smoothedRandom(n, spec.smoothBands, spec.seed, spec.strength)
83 case 'delay': {
84 const tau = spec.delayMs / 1000
85 return centers.map((f) => -2 * Math.PI * f * tau)
86 }
87 case 'chirp': {
88 // tau(f) = tauN f / fN => phi(f) = -2 pi int tau df = -pi tauN f^2 / fN
89 const tauN = spec.delayAtNyquistMs / 1000
90 const fN = fs / 2
91 return centers.map((f) => (-Math.PI * tauN * f * f) / fN)
92 }
93 }
94}
96/** Largest band-to-band slope of the map, as a group delay in seconds. Used to size the FFT padding. */
97export function maxGroupDelaySec(phi: Float64Array, df: number): number {
98 let m = 0
99 for (let i = 1; i < phi.length; i++) {
100 let d = Math.abs(phi[i] - phi[i - 1])
101 d = Math.min(d, 2 * Math.PI) // a jump larger than 2 pi wraps to something smaller anyway
102 if (d > m) m = d
103 }
104 return m / (2 * Math.PI * df)
105}
107export const wrapAngle = (x: number): number => {
108 let y = x % (2 * Math.PI)
109 if (y > Math.PI) y -= 2 * Math.PI
110 if (y <= -Math.PI) y += 2 * Math.PI
111 return y
112}