7d910acPhase perception 2: subband phase rotation of arbitrary sounds with filter bank, phase maps, plots, and playbackJeremy Magland 1/**
2 * Gaussian filter bank: the short-time Fourier transform with a hop of one
3 * sample, viewed as a bank of complex bandpass filters.
4 *
5 * Band k has center f_k = k * df and frequency response
6 *
7 * H_k(f) = W(f - f_k) / S(f), W(f) = exp(-f^2 / (2 sigma_f^2)), S(f) = sum_j W(f - f_j)
8 *
9 * where sigma_f = 1 / (2 pi sigma) is the spectral width of a Gaussian window
10 * of time width sigma. The responses sum to one on f >= 0, so the subbands
11 * z_k = x * h_k add back to the signal: x = Re sum_k z_k.
12 *
13 * Rotating every subband by a constant angle phi_k and adding them back is a
14 * linear time-invariant filter with response G(f) = sum_k exp(i phi_k) H_k(f).
15 * We therefore compute G on the FFT grid and apply it with one FFT pair,
16 * which is exact and much cheaper than synthesizing the subbands one by one.
17 */
18import { fftPlan, nextPow2 } from './fft'
20export const sigmaF = (sigmaSec: number): number => 1 / (2 * Math.PI * sigmaSec)
22/** Band centers 0, df, 2 df, ... up to the Nyquist frequency, with df = sigma_f. */
23export function bandCenters(fs: number, sigmaSec: number): Float64Array {
24 const df = sigmaF(sigmaSec)
25 const k = Math.ceil(fs / 2 / df)
26 const c = new Float64Array(k + 1)
27 for (let i = 0; i <= k; i++) c[i] = i * df
28 return c
29}
31export type Response = { re: Float64Array; im: Float64Array; binHz: number }
33/** G(f) on the bins 0..L/2 of an FFT of length L at sample rate fs. */
34export function effectiveFilter(centers: Float64Array, phi: ArrayLike<number>, sigmaSec: number, L: number, fs: number): Response {
35 const sf = sigmaF(sigmaSec)
36 const binHz = fs / L
37 const nb = L / 2 + 1
38 const S = new Float64Array(nb)
39 const re = new Float64Array(nb)
40 const im = new Float64Array(nb)
41 const reach = 5 * sf
42 const inv2s2 = 1 / (2 * sf * sf)
43 for (let k = 0; k < centers.length; k++) {
44 const fc = centers[k]
45 const c = Math.cos(phi[k])
46 const s = Math.sin(phi[k])
47 const i0 = Math.max(0, Math.ceil((fc - reach) / binHz))
48 const i1 = Math.min(nb - 1, Math.floor((fc + reach) / binHz))
49 for (let i = i0; i <= i1; i++) {
50 const d = i * binHz - fc
51 const w = Math.exp(-d * d * inv2s2)
52 S[i] += w
53 re[i] += c * w
54 im[i] += s * w
55 }
56 }
57 for (let i = 0; i < nb; i++) {
58 if (S[i] > 0) {
59 re[i] /= S[i]
60 im[i] /= S[i]
61 }
62 }
63 return { re, im, binHz }
64}
66/** Group delay in seconds, -(d/d omega) arg G, by finite differences between neighbouring bins. */
67export function groupDelay(g: Response): Float64Array {
68 const n = g.re.length
69 const tau = new Float64Array(n)
70 for (let i = 0; i + 1 < n; i++) {
71 // angle of G[i+1] * conj(G[i]), already wrapped to (-pi, pi]
72 const re = g.re[i + 1] * g.re[i] + g.im[i + 1] * g.im[i]
73 const im = g.im[i + 1] * g.re[i] - g.re[i + 1] * g.im[i]
74 tau[i] = -Math.atan2(im, re) / (2 * Math.PI * g.binHz)
75 }
76 tau[n - 1] = tau[n - 2]
77 return tau
78}
80/**
81 * Apply the filter G to x. The signal is placed in the middle of a zero-padded
82 * buffer of length L so that the circular convolution does not wrap; padSec
83 * must cover the window tails plus the largest group delay.
84 */
85export function applyFilter(x: Float32Array, fs: number, sigmaSec: number, phi: ArrayLike<number>, centers: Float64Array, padSec: number): { y: Float32Array; g: Response } {
86 const n = x.length
87 const pad = Math.ceil(padSec * fs)
88 const L = nextPow2(n + 2 * pad)
89 const g = effectiveFilter(centers, phi, sigmaSec, L, fs)
90 const re = new Float64Array(L)
91 const im = new Float64Array(L)
92 re.set(x, pad)
93 const plan = fftPlan(L)
94 plan.transform(re, im)
95 const half = L / 2
96 // Y = X G on the positive half; the negative half is its conjugate mirror so y is real.
97 for (let i = 1; i < half; i++) {
98 const xr = re[i]
99 const xi = im[i]
100 const yr = xr * g.re[i] - xi * g.im[i]
101 const yi = xr * g.im[i] + xi * g.re[i]
102 re[i] = yr
103 im[i] = yi
104 re[L - i] = yr
105 im[L - i] = -yi
106 }
107 re[0] *= g.re[0]
108 im[0] = 0
109 re[half] *= g.re[half]
110 im[half] = 0
111 plan.transform(re, im, true)
112 const y = new Float32Array(n)
113 for (let i = 0; i < n; i++) y[i] = re[pad + i]
114 return { y, g }
115}
117export type Spectrogram = {
118 /** rows x cols, row-major, row 0 at 0 Hz; linear magnitude */
119 data: Float32Array
120 rows: number
121 cols: number
122 durationSec: number
123 fmax: number
124}
126/**
127 * Magnitude spectrogram with a Gaussian window of the given sigma, sampled on
128 * a grid of at most maxCols columns and maxRows rows (bins are max-pooled into
129 * rows). This is |z_k(t)| of the bank above, up to the normalization by S.
130 */
131export function spectrogram(x: Float32Array, fs: number, sigmaSec: number, fmax: number, maxCols: number, maxRows: number): Spectrogram {
132 const n = x.length
133 const ss = sigmaSec * fs
134 const half = Math.ceil(4 * ss)
135 const wlen = 2 * half + 1
136 const N = nextPow2(Math.max(wlen, 2048))
137 const hop = Math.max(Math.ceil(n / maxCols), Math.ceil(ss / 2), 1)
138 const cols = Math.ceil(n / hop)
139 const w = new Float64Array(wlen)
140 for (let j = 0; j < wlen; j++) w[j] = Math.exp(-0.5 * ((j - half) / ss) ** 2)
141 const binHz = fs / N
142 const nb = Math.min(N / 2, Math.floor(fmax / binHz))
143 const rows = Math.min(maxRows, nb + 1)
144 const data = new Float32Array(rows * cols)
145 const plan = fftPlan(N)
146 const re = new Float64Array(N)
147 const im = new Float64Array(N)
148 const mag = new Float64Array(nb + 1)
149 for (let c = 0; c < cols; c++) {
150 const center = c * hop
151 re.fill(0)
152 im.fill(0)
153 for (let j = 0; j < wlen; j++) {
154 const t = center - half + j
155 if (t >= 0 && t < n) re[j] = x[t] * w[j]
156 }
157 plan.transform(re, im)
158 for (let i = 0; i <= nb; i++) mag[i] = Math.hypot(re[i], im[i])
159 for (let r = 0; r < rows; r++) {
160 const b0 = Math.floor((r * (nb + 1)) / rows)
161 const b1 = Math.max(b0 + 1, Math.floor(((r + 1) * (nb + 1)) / rows))
162 let m = 0
163 for (let i = b0; i < b1 && i <= nb; i++) if (mag[i] > m) m = mag[i]
164 data[r * cols + c] = m
165 }
166 }
167 return { data, rows, cols, durationSec: n / fs, fmax: rows * ((nb + 1) / rows) * binHz }
168}
e9f92d1Time-varying random phase maps by decimated subband synthesis; drag-to-zoom, scrollbar pan and zoom buttons for the waveformsJeremy Magland 169
170/** S(f) = sum_k W(f - f_k) on the bins 0..L/2, the normalizer that makes the bank sum to one. */
171export function bandSum(centers: Float64Array, sigmaSec: number, L: number, fs: number): Float64Array {
172 const sf = sigmaF(sigmaSec)
173 const binHz = fs / L
174 const nb = L / 2 + 1
175 const S = new Float64Array(nb)
176 const reach = 5 * sf
177 const inv2s2 = 1 / (2 * sf * sf)
178 for (let k = 0; k < centers.length; k++) {
179 const fc = centers[k]
180 const i0 = Math.max(0, Math.ceil((fc - reach) / binHz))
181 const i1 = Math.min(nb - 1, Math.floor((fc + reach) / binHz))
182 for (let i = i0; i <= i1; i++) {
183 const d = i * binHz - fc
184 S[i] += Math.exp(-d * d * inv2s2)
185 }
186 }
187 return S
188}
190const scratch = new Map<number, { re: Float64Array; im: Float64Array }>()
191const scratchFor = (m: number) => {
192 let s = scratch.get(m)
193 if (!s) {
194 s = { re: new Float64Array(m), im: new Float64Array(m) }
195 scratch.set(m, s)
196 }
197 s.re.fill(0)
198 s.im.fill(0)
199 return s
200}
202/**
203 * y = Re sum_k exp(i phi_k(t)) z_k(t) with a phase that varies in time, by
204 * exact subband synthesis. Each band is narrow, so its subband signal is
205 * computed at a reduced sample rate: the band's slice of the spectrum is
206 * demodulated to baseband with a small inverse FFT, rotated sample by sample,
207 * transformed back, and added into the output spectrum at its original bins.
208 * The small transform is sized to hold the band plus the frequency shift the
209 * rotation can impose (maxRateHz), so nothing aliases.
210 *
211 * phiAt(k, times) must return the phase of band k at the given times, in
212 * seconds from the start of the padded buffer (the signal starts at padSec).
213 */
214export function applyTimeVarying(
215 x: Float32Array,
216 fs: number,
217 sigmaSec: number,
218 centers: Float64Array,
219 padSec: number,
220 phiAt: (k: number, times: Float64Array) => Float64Array,
221 maxRateHz: number,
222): Float32Array {
223 const n = x.length
224 const pad = Math.ceil(padSec * fs)
225 const L = nextPow2(n + 2 * pad)
226 const half = L / 2
227 const binHz = fs / L
228 const sf = sigmaF(sigmaSec)
229 const reach = 5 * sf
230 const inv2s2 = 1 / (2 * sf * sf)
231 const S = bandSum(centers, sigmaSec, L, fs)
233 const re = new Float64Array(L)
234 const im = new Float64Array(L)
235 re.set(x, pad)
236 fftPlan(L).transform(re, im)
237 // Positive half of the spectrum. The rotated subbands are one-sided, and the real part of a
238 // one-sided signal has half its spectrum on each side, so the positive half of X (not the
239 // doubled analytic spectrum) is what the Hermitian mirror below needs.
240 const xr = new Float64Array(half + 1)
241 const xi = new Float64Array(half + 1)
242 for (let i = 0; i <= half; i++) {
243 xr[i] = re[i]
244 xi[i] = im[i]
245 }
247 const yr = new Float64Array(L)
248 const yi = new Float64Array(L)
249 const extra = Math.ceil((1.5 * maxRateHz) / binHz) + 16
250 for (let k = 0; k < centers.length; k++) {
251 const fc = centers[k]
252 const i0 = Math.max(0, Math.ceil((fc - reach) / binHz))
253 const i1 = Math.min(half, Math.floor((fc + reach) / binHz))
254 if (i1 < i0) continue
255 const M = Math.min(L, nextPow2(i1 - i0 + 1 + 2 * extra))
256 const ic = Math.round(fc / binHz)
257 const { re: sr, im: si } = scratchFor(M)
258 for (let i = i0; i <= i1; i++) {
259 const d = i * binHz - fc
260 const w = S[i] > 0 ? Math.exp(-d * d * inv2s2) / S[i] : 0
261 const j = (((i - ic) % M) + M) % M
262 sr[j] = xr[i] * w
263 si[j] = xi[i] * w
264 }
265 const plan = fftPlan(M)
266 plan.transform(sr, si, true) // z_k demodulated, sampled every L/M samples
267 const times = new Float64Array(M)
268 for (let m = 0; m < M; m++) times[m] = (m * L) / (M * fs)
269 const phi = phiAt(k, times)
270 for (let m = 0; m < M; m++) {
271 const c = Math.cos(phi[m])
272 const s = Math.sin(phi[m])
273 const a = sr[m]
274 const b = si[m]
275 sr[m] = a * c - b * s
276 si[m] = a * s + b * c
277 }
278 plan.transform(sr, si)
279 for (let j = 0; j < M; j++) {
280 const off = j < M / 2 ? j : j - M
281 const i = ic + off
282 if (i < 0 || i > half) continue
283 yr[i] += sr[j]
284 yi[i] += si[j]
285 }
286 }
287 // Hermitian mirror so the output is real
288 yi[0] = 0
289 yi[half] = 0
290 for (let i = 1; i < half; i++) {
291 yr[L - i] = yr[i]
292 yi[L - i] = -yi[i]
293 }
294 fftPlan(L).transform(yr, yi, true)
295 const y = new Float32Array(n)
296 for (let i = 0; i < n; i++) y[i] = yr[pad + i]
297 return y
298}