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}