Add documentation figures
11 changed files+571−0
docs/figures/01-stft-basics.pngadded+0−0View file
Binary file not shown.
docs/figures/02-filterbank.pngadded+0−0View file
Binary file not shown.
docs/figures/03-phase-rotation.pngadded+0−0View file
Binary file not shown.
docs/figures/04-spectrograms.pngadded+0−0View file
Binary file not shown.
docs/figures/05-group-delay.pngadded+0−0View file
Binary file not shown.
docs/figures/06-notches.pngadded+0−0View file
Binary file not shown.
docs/figures/07-inconsistency.pngadded+0−0View file
Binary file not shown.
docs/figures/08-time-varying.pngadded+0−0View file
Binary file not shown.
docs/figures/09-pipeline.pngadded+0−0View file
Binary file not shown.
docs/make_figures.pyadded+462−0View file
@@ -0,0 +1,462 @@
1+#!/usr/bin/env python3
2+"""Generate the figures for spectrogram-phase-proposal.md.
3+
4+Run from anywhere: python3 docs/make_figures.py
5+Figures are written to docs/figures/.
6+"""
7+from pathlib import Path
8+
9+import matplotlib
10+
11+matplotlib.use("Agg")
12+import matplotlib.pyplot as plt
13+import numpy as np
14+from scipy.ndimage import gaussian_filter1d
15+from matplotlib.patches import FancyArrowPatch, FancyBboxPatch
16+
17+OUT = Path(__file__).parent / "figures"
18+OUT.mkdir(exist_ok=True)
19+rng = np.random.default_rng(3)
20+
21+FS = 16000
22+T = 0.8
23+N = int(FS * T)
24+t = np.arange(N) / FS
25+F0 = 220
26+CLICK_T = 0.15
27+TONE_T0, TONE_T1 = 0.35, 0.7
28+
29+COL_A = "#1f5fa8"
30+COL_B = "#d9641e"
31+GRAY = "#888888"
32+
33+plt.rcParams.update({"font.size": 9, "axes.titlesize": 10, "figure.dpi": 130})
34+
35+
36+# ----------------------------------------------------------------------------
37+# Test signal: a click followed by a harmonic tone
38+# ----------------------------------------------------------------------------
39+def test_signal():
40+ x = np.zeros(N)
41+ i0 = int(CLICK_T * FS)
42+ x[i0 : i0 + 4] += np.hanning(6)[1:-1]
43+ seg = (t >= TONE_T0) & (t < TONE_T1)
44+ tt = t[seg] - TONE_T0
45+ env = np.ones_like(tt)
46+ fade = int(0.01 * FS)
47+ env[:fade] = 0.5 - 0.5 * np.cos(np.pi * np.arange(fade) / fade)
48+ env[-fade:] = env[:fade][::-1]
49+ tone = sum(np.sin(2 * np.pi * k * F0 * tt) / k for k in range(1, 9))
50+ x[seg] += 0.25 * env * tone
51+ return x
52+
53+
54+x = test_signal()
55+
56+
57+# ----------------------------------------------------------------------------
58+# Gaussian filter bank (the STFT with a hop of one sample)
59+# ----------------------------------------------------------------------------
60+class GaussBank:
61+ """Complex Gaussian bandpass filters whose responses sum to one on f >= 0.
62+
63+ Band k has center f_k = k * df and response W(f - f_k) / S(f), where
64+ W is a Gaussian of width sigma_f = 1 / (2 pi sigma_t) and
65+ S(f) = sum_k W(f - f_k). The subband z_k = x * h_k is the STFT at bin k
66+ evaluated at every sample. Because the responses sum to one,
67+ x = Re sum_k z_k exactly (positive frequencies are doubled first).
68+ """
69+
70+ def __init__(self, n, fs, sigma_t, spacing=1.0):
71+ self.n, self.fs, self.sigma_t = n, fs, sigma_t
72+ pad = int(8 * sigma_t * fs)
73+ self.L = 1 << int(np.ceil(np.log2(n + 2 * pad)))
74+ self.f = np.fft.fftfreq(self.L, 1 / fs)
75+ self.sigma_f = 1 / (2 * np.pi * sigma_t)
76+ self.df = spacing * self.sigma_f
77+ self.centers = np.arange(0, fs / 2 + self.df / 2, self.df)
78+ S = np.zeros(self.L)
79+ for fc in self.centers:
80+ S += self._W(fc)
81+ self.S = np.where(S > 1e-300, S, np.inf) # far outside the covered range the sum underflows; those bins are never used
82+
83+ def _W(self, fc):
84+ return np.exp(-((self.f - fc) ** 2) / (2 * self.sigma_f**2))
85+
86+ def response(self, k):
87+ return self._W(self.centers[k]) / self.S
88+
89+ def analytic_spectrum(self, x):
90+ X = np.fft.fft(x, self.L)
91+ X[self.f < 0] = 0
92+ X[self.f > 0] *= 2
93+ return X
94+
95+ def subband(self, Xa, k):
96+ return np.fft.ifft(Xa * self.response(k))[: self.n]
97+
98+ def apply(self, x, phase=None):
99+ """y = Re sum_k exp(i phi_k(t)) z_k(t). phase(k) returns a scalar or an array over time."""
100+ Xa = self.analytic_spectrum(x)
101+ y = np.zeros(self.n)
102+ for k in range(len(self.centers)):
103+ z = self.subband(Xa, k)
104+ if phase is not None:
105+ z = z * np.exp(1j * phase(k))
106+ y += z.real
107+ return y
108+
109+ def spectrogram(self, x, hop, fmax):
110+ Xa = self.analytic_spectrum(x)
111+ ks = [k for k, fc in enumerate(self.centers) if fc <= fmax]
112+ idx = np.arange(0, self.n, hop)
113+ M = np.zeros((len(ks), len(idx)))
114+ for i, k in enumerate(ks):
115+ M[i] = np.abs(self.subband(Xa, k)[idx])
116+ return M, self.centers[ks], idx / self.fs
117+
118+ def effective_filter(self, phi):
119+ """Frequency response of the LTI filter obtained with constant per-band phases phi_k."""
120+ G = np.zeros(self.L, dtype=complex)
121+ for k, p in enumerate(phi):
122+ G += np.exp(1j * p) * self.response(k)
123+ return G
124+
125+
126+def smooth_random_phases(nbands, smooth_bands, rng):
127+ """Uniform random phases, smoothed across bands with a Gaussian of the given width, then rescaled."""
128+ raw = rng.uniform(-np.pi, np.pi, nbands)
129+ if smooth_bands <= 0:
130+ return raw
131+ sm = gaussian_filter1d(raw, smooth_bands, mode="reflect")
132+ return sm * raw.std() / sm.std()
133+
134+
135+def db(M):
136+ return 20 * np.log10(np.maximum(M, 1e-12) / M.max())
137+
138+
139+def show_spec(ax, M, freqs, times, title, vmin=-60):
140+ ax.imshow(db(M), origin="lower", aspect="auto", extent=[times[0], times[-1], freqs[0], freqs[-1]], vmin=vmin, vmax=0, cmap="magma")
141+ ax.set_title(title)
142+ ax.set_xlabel("time (s)")
143+ ax.set_ylabel("frequency (Hz)")
144+
145+
146+# ----------------------------------------------------------------------------
147+# Figure 1: the STFT as sliding windows, and the spectrogram
148+# ----------------------------------------------------------------------------
149+def fig_stft_basics():
150+ sigma = 0.005
151+ fig, axs = plt.subplots(2, 1, figsize=(8, 5.2), gridspec_kw={"height_ratios": [1, 1.4]})
152+ ax = axs[0]
153+ ax.plot(t, x, color="k", lw=0.6)
154+ hop = 0.02
155+ for i, c in enumerate(np.arange(0.06, 0.56, hop)):
156+ w = np.exp(-((t - c) ** 2) / (2 * sigma**2))
157+ ax.plot(t, 1.05 * w, color=COL_A if i % 2 else COL_B, lw=0.9, alpha=0.8)
158+ ax.annotate("hop", xy=(0.30, 1.12), xytext=(0.30, 1.35), ha="center", arrowprops=dict(arrowstyle="-[, widthB=1.6", lw=0.8))
159+ ax.set_xlim(0, T)
160+ ax.set_ylim(-1.1, 1.5)
161+ ax.set_ylabel("amplitude")
162+ ax.set_title("Test signal (a click, then a harmonic tone) with Gaussian analysis windows, sigma = 5 ms")
163+ bank = GaussBank(N, FS, sigma)
164+ M, freqs, times = bank.spectrogram(x, 32, 4000)
165+ show_spec(axs[1], M, freqs, times, "Spectrogram: magnitude of each windowed Fourier coefficient (dB)")
166+ fig.tight_layout()
167+ fig.savefig(OUT / "01-stft-basics.png")
168+ plt.close(fig)
169+
170+
171+# ----------------------------------------------------------------------------
172+# Figure 2: the filter bank view
173+# ----------------------------------------------------------------------------
174+def fig_filterbank():
175+ sigma = 0.005
176+ bank = GaussBank(N, FS, sigma)
177+ fig = plt.figure(figsize=(8, 7.5))
178+ gs = fig.add_gridspec(4, 1, height_ratios=[1, 1, 1, 1])
179+
180+ ax = fig.add_subplot(gs[0])
181+ pos = bank.f >= 0
182+ f = bank.f[pos]
183+ total = np.zeros(pos.sum())
184+ for k in range(len(bank.centers)):
185+ r = bank.response(k)[pos]
186+ total += r
187+ if bank.centers[k] <= 1500:
188+ ax.plot(f, r, lw=0.8, color=COL_A if k % 2 else COL_B, alpha=0.8)
189+ ax.plot(f, total, "k--", lw=1, label="sum over all bands")
190+ ax.set_xlim(0, 1500)
191+ ax.set_ylim(0, 1.15)
192+ ax.set_xlabel("frequency (Hz)")
193+ ax.set_ylabel("gain")
194+ ax.legend(loc="lower right")
195+ ax.set_title("Each STFT bin is a bandpass filter. The filters overlap and sum to one.")
196+
197+ Xa = bank.analytic_spectrum(x)
198+ picks = [int(round(F0 * m / bank.df)) for m in (1, 3)] + [int(round(1200 / bank.df))]
199+ for row, k in enumerate(picks, start=1):
200+ ax = fig.add_subplot(gs[row])
201+ z = bank.subband(Xa, k)
202+ ax.plot(t, z.real, color=COL_A, lw=0.6, label="Re z_k(t)")
203+ ax.plot(t, np.abs(z), color=COL_B, lw=1.0, label="|z_k(t)| (envelope)")
204+ ax.set_xlim(0, T)
205+ ax.set_ylabel("amplitude")
206+ ax.set_title(f"Subband k = {k}, center {bank.centers[k]:.0f} Hz", loc="left")
207+ if row == 1:
208+ ax.legend(loc="upper right", ncol=2)
209+ ax.set_xlabel("time (s)")
210+ y = bank.apply(x)
211+ fig.text(0.5, 0.005, f"Adding all subbands back together reproduces the signal: max |x - Re sum_k z_k| = {np.max(np.abs(x - y)):.1e}", ha="center", fontsize=9)
212+ fig.tight_layout(rect=(0, 0.02, 1, 1))
213+ fig.savefig(OUT / "02-filterbank.png")
214+ plt.close(fig)
215+
216+
217+# ----------------------------------------------------------------------------
218+# Figure 3: rotating subband phases, for three window lengths
219+# ----------------------------------------------------------------------------
220+SMOOTH = 6 # bands
221+
222+
223+def rotated_versions(sigmas, smooth_bands=SMOOTH, seed=5):
224+ out = {}
225+ for s in sigmas:
226+ bank = GaussBank(N, FS, s)
227+ phi = smooth_random_phases(len(bank.centers), smooth_bands, np.random.default_rng(seed))
228+ y = bank.apply(x, lambda k: phi[k])
229+ out[s] = (bank, phi, y)
230+ return out
231+
232+
233+ROT = rotated_versions([0.0005, 0.005, 0.02])
234+
235+
236+def fig_phase_rotation():
237+ fig, axs = plt.subplots(3, 2, figsize=(8, 7), gridspec_kw={"width_ratios": [1.4, 1]})
238+ for row, (s, (bank, phi, y)) in enumerate(ROT.items()):
239+ ax = axs[row, 0]
240+ sel = (t > CLICK_T - 0.05) & (t < CLICK_T + 0.08)
241+ ax.plot(t[sel] * 1000, x[sel], color=GRAY, lw=0.8, label="original")
242+ ax.plot(t[sel] * 1000, y[sel], color=COL_B, lw=0.8, label="phases rotated")
243+ ax.set_title(f"window sigma = {s*1000:g} ms: the click", loc="left")
244+ ax.set_ylabel("amplitude")
245+ if row == 0:
246+ ax.legend(loc="upper right")
247+ if row == 2:
248+ ax.set_xlabel("time (ms)")
249+ ax = axs[row, 1]
250+ sel = (t > 0.5) & (t < 0.5 + 3 / F0)
251+ ax.plot(t[sel] * 1000, x[sel], color=GRAY, lw=0.8)
252+ ax.plot(t[sel] * 1000, y[sel], color=COL_B, lw=0.8)
253+ ax.set_title("same map, steady tone (3 periods)", loc="left")
254+ if row == 2:
255+ ax.set_xlabel("time (ms)")
256+ fig.suptitle("y = Re sum_k exp(i phi_k) z_k(t) with smooth random phi_k", fontsize=10)
257+ fig.tight_layout(rect=(0, 0, 1, 0.96))
258+ fig.savefig(OUT / "03-phase-rotation.png")
259+ plt.close(fig)
260+
261+
262+# ----------------------------------------------------------------------------
263+# Figure 4: spectrograms before and after, at two analysis resolutions
264+# ----------------------------------------------------------------------------
265+def fig_spectrograms():
266+ _, _, y = ROT[0.02]
267+ fig, axs = plt.subplots(2, 2, figsize=(8, 5.6))
268+ for col, s in enumerate([0.02, 0.002]):
269+ bank = GaussBank(N, FS, s)
270+ for row, (sig, name) in enumerate([(x, "original"), (y, "rotated (20 ms bank)")]):
271+ M, freqs, times = bank.spectrogram(sig, 32, 4000)
272+ show_spec(axs[row, col], M, freqs, times, f"{name}, analysed at sigma = {s*1000:g} ms")
273+ fig.tight_layout()
274+ fig.savefig(OUT / "04-spectrograms.png")
275+ plt.close(fig)
276+
277+
278+# ----------------------------------------------------------------------------
279+# Figure 5: the equivalent all-pass filter and its group delay
280+# ----------------------------------------------------------------------------
281+def fig_group_delay():
282+ fig, axs = plt.subplots(3, 1, figsize=(8, 7), sharex=True)
283+ for s, (bank, phi, _) in ROT.items():
284+ G = bank.effective_filter(phi)
285+ pos = (bank.f > 0) & (bank.f <= 4000)
286+ f = bank.f[pos]
287+ ph = np.unwrap(np.angle(G[pos]))
288+ tau = -np.gradient(ph, f) / (2 * np.pi)
289+ lab = f"sigma = {s*1000:g} ms"
290+ axs[0].plot(f, 20 * np.log10(np.abs(G[pos])), lw=0.8, label=lab)
291+ axs[1].plot(f, ph, lw=0.8, label=lab)
292+ axs[2].plot(f, tau * 1000, lw=0.8, label=lab)
293+ axs[0].set_ylabel("magnitude (dB)")
294+ axs[0].set_ylim(-3, 1)
295+ axs[0].set_title("Constant per-band phases act as a filter G(f): magnitude (top), phase (middle), group delay (bottom)")
296+ axs[0].legend(loc="lower left", ncol=3, fontsize=8)
297+ axs[1].set_ylabel("phase (rad)")
298+ axs[2].axhspan(-1.5, 1.5, color="green", alpha=0.12)
299+ axs[2].text(0.99, 0.95, "green band: roughly the audibility threshold for group delay peaks", ha="right", va="top", transform=axs[2].transAxes, fontsize=8, color="green")
300+ axs[2].set_ylabel("group delay (ms)")
301+ axs[2].set_xlabel("frequency (Hz)")
302+ axs[2].set_xlim(0, 4000)
303+ fig.tight_layout()
304+ fig.savefig(OUT / "05-group-delay.png")
305+ plt.close(fig)
306+
307+
308+# ----------------------------------------------------------------------------
309+# Figure 6: why the phases must vary smoothly across neighbouring bands
310+# ----------------------------------------------------------------------------
311+def fig_notches():
312+ bank = GaussBank(N, FS, 0.005)
313+ fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex="col")
314+ for col, (sb, name) in enumerate([(0, "independent random phase per band"), (SMOOTH, f"random phases smoothed over {SMOOTH} bands")]):
315+ phi = smooth_random_phases(len(bank.centers), sb, np.random.default_rng(5))
316+ G = bank.effective_filter(phi)
317+ pos = (bank.f > 0) & (bank.f <= 2000)
318+ f = bank.f[pos]
319+ ax = axs[0, col]
320+ ax.plot(bank.centers, np.angle(np.exp(1j * phi)), ".", ms=3, color=COL_B)
321+ ax.set_ylabel("phi_k (rad)")
322+ ax.set_title(name)
323+ ax.set_xlim(0, 2000)
324+ ax = axs[1, col]
325+ ax.plot(f, 20 * np.log10(np.abs(G[pos])), color=COL_A, lw=0.8)
326+ ax.set_ylim(-30, 3)
327+ ax.set_ylabel("|G(f)| (dB)")
328+ ax.set_xlabel("frequency (Hz)")
329+ for m in range(1, 9):
330+ ax.axvline(m * F0, color=GRAY, lw=0.5, ls=":")
331+ axs[1, 0].set_title("|G(f)|: notches where neighbours disagree", fontsize=9)
332+ axs[1, 1].set_title("|G(f)|: nearly flat (dotted: harmonics of the tone)", fontsize=9)
333+ fig.tight_layout()
334+ fig.savefig(OUT / "06-notches.png")
335+ plt.close(fig)
336+
337+
338+# ----------------------------------------------------------------------------
339+# Figure 7: hop-based STFT, per-cell random phase, inverse STFT
340+# ----------------------------------------------------------------------------
341+def stft(x, n, hop):
342+ w = np.hanning(n + 1)[:-1]
343+ frames = [x[i : i + n] * w for i in range(0, len(x) - n + 1, hop)]
344+ return np.fft.rfft(np.array(frames), axis=1), w
345+
346+
347+def istft(S, w, hop, length):
348+ n = len(w)
349+ y = np.zeros(length)
350+ norm = np.zeros(length)
351+ for m, spec in enumerate(S):
352+ i = m * hop
353+ y[i : i + n] += np.fft.irfft(spec, n) * w
354+ norm[i : i + n] += w * w
355+ return y / np.maximum(norm, 1e-9)
356+
357+
358+def fig_inconsistency():
359+ n, hop = 512, 128
360+ S, w = stft(x, n, hop)
361+ freqs = np.fft.rfftfreq(n, 1 / FS)
362+ keep = freqs <= 4000
363+ times = np.arange(S.shape[0]) * hop / FS
364+
365+ S_cell = S * np.exp(1j * rng.uniform(-np.pi, np.pi, S.shape))
366+ y_cell = istft(S_cell, w, hop, N)
367+ phi_band = smooth_random_phases(S.shape[1], SMOOTH, np.random.default_rng(5))
368+ S_band = S * np.exp(1j * phi_band)[None, :]
369+ y_band = istft(S_band, w, hop, N)
370+
371+ fig, axs = plt.subplots(2, 3, figsize=(9.5, 5.4), gridspec_kw={"height_ratios": [1, 1.6]})
372+ for col, (sig, title) in enumerate([(x, "original"), (y_cell, "after random phase in every cell"), (y_band, "after one phase per frequency bin")]):
373+ axs[0, col].plot(t, sig, color="k", lw=0.5)
374+ axs[0, col].set_xlim(0, T)
375+ axs[0, col].set_ylim(-1.1, 1.1)
376+ axs[0, col].set_title(title)
377+ M = np.abs(stft(sig, n, hop)[0])[:, keep].T
378+ show_spec(axs[1, col], M, freqs[keep], times, "spectrogram recomputed from the signal")
379+ fig.suptitle("Hop-based STFT (Hann, 32 ms, 75% overlap), modify phases, inverse STFT by overlap-add", fontsize=10)
380+ fig.tight_layout(rect=(0, 0, 1, 0.96))
381+ fig.savefig(OUT / "07-inconsistency.png")
382+ plt.close(fig)
383+
384+
385+# ----------------------------------------------------------------------------
386+# Figure 8: a time-varying phase offset is frequency modulation
387+# ----------------------------------------------------------------------------
388+def fig_time_varying():
389+ fc = 0.5
390+ delta = np.pi
391+ tt = np.linspace(0, 2 / fc, 2000)
392+ s = 0.5 - 0.5 * np.cos(2 * np.pi * fc * tt)
393+ phi = delta * s
394+ dphi = np.gradient(phi, tt)
395+ fig, axs = plt.subplots(1, 3, figsize=(9.5, 2.9))
396+ axs[0].plot(tt, phi, color=COL_A)
397+ axs[0].set_title("phi(t): 0 to pi and back at 0.5 Hz")
398+ axs[0].set_xlabel("time (s)")
399+ axs[0].set_ylabel("rad")
400+ axs[1].plot(tt, dphi / (2 * np.pi), color=COL_B)
401+ axs[1].set_title("frequency shift phi'(t) / 2 pi")
402+ axs[1].set_xlabel("time (s)")
403+ axs[1].set_ylabel("Hz")
404+ fcs = np.logspace(-2, 2, 100)
405+ axs[2].loglog(fcs, delta * fcs / 2, color=COL_B)
406+ axs[2].axhspan(0.5, 2, color="green", alpha=0.12)
407+ axs[2].text(0.011, 2.4, "roughly the pitch JND near 220 Hz", fontsize=8, color="green")
408+ axs[2].set_title("peak shift vs cycle rate (swing pi)")
409+ axs[2].set_xlabel("cycle rate (Hz)")
410+ axs[2].set_ylabel("peak shift (Hz)")
411+ fig.tight_layout()
412+ fig.savefig(OUT / "08-time-varying.png")
413+ plt.close(fig)
414+
415+
416+# ----------------------------------------------------------------------------
417+# Figure 9: proposed processing pipeline
418+# ----------------------------------------------------------------------------
419+def fig_pipeline():
420+ fig, ax = plt.subplots(figsize=(9.5, 3.4))
421+ ax.set_xlim(0, 10)
422+ ax.set_ylim(0, 3.4)
423+ ax.axis("off")
424+
425+ def box(xc, yc, wdt, hgt, text, color="#eef3fa", ec="#333"):
426+ ax.add_patch(FancyBboxPatch((xc - wdt / 2, yc - hgt / 2), wdt, hgt, boxstyle="round,pad=0.05", fc=color, ec=ec, lw=0.9))
427+ ax.text(xc, yc, text, ha="center", va="center", fontsize=8.5)
428+
429+ def arrow(x0, y0, x1, y1):
430+ ax.add_patch(FancyArrowPatch((x0, y0), (x1, y1), arrowstyle="-|>", mutation_scale=12, lw=0.9, color="#333"))
431+
432+ box(0.9, 1.7, 1.5, 1.2, "source sound x(t)\n(click, noise burst,\npluck, file, mic)")
433+ box(3.1, 1.7, 1.7, 1.2, "Gaussian filter bank\nwindow sigma\n\nz_k(t) for each band k")
434+ box(5.55, 2.6, 1.9, 0.75, "phase map A: phi_k = 0", color="#e5edf8", ec=COL_A)
435+ box(5.55, 0.8, 1.9, 0.75, "phase map B: phi_k chosen", color="#fbeadf", ec=COL_B)
436+ box(7.9, 2.6, 1.7, 0.75, "A = Re sum z_k = x", color="#e5edf8", ec=COL_A)
437+ box(7.9, 0.8, 1.7, 0.75, "B = Re sum e^(i phi_k) z_k", color="#fbeadf", ec=COL_B)
438+ box(9.4, 1.7, 1.1, 1.0, "blind test\nA/B, cycling")
439+ arrow(1.65, 1.7, 2.25, 1.7)
440+ arrow(3.95, 1.9, 4.6, 2.6)
441+ arrow(3.95, 1.5, 4.6, 0.8)
442+ arrow(6.5, 2.6, 7.05, 2.6)
443+ arrow(6.5, 0.8, 7.05, 0.8)
444+ arrow(8.75, 2.6, 9.0, 2.0)
445+ arrow(8.75, 0.8, 9.0, 1.4)
446+ ax.text(5.55, 1.7, "same subband envelopes |z_k(t)|,\ndifferent subband phases", ha="center", va="center", fontsize=8, color="#555")
447+ fig.tight_layout()
448+ fig.savefig(OUT / "09-pipeline.png")
449+ plt.close(fig)
450+
451+
452+if __name__ == "__main__":
453+ fig_stft_basics()
454+ fig_filterbank()
455+ fig_phase_rotation()
456+ fig_spectrograms()
457+ fig_group_delay()
458+ fig_notches()
459+ fig_inconsistency()
460+ fig_time_varying()
461+ fig_pipeline()
462+ print("wrote", sorted(p.name for p in OUT.iterdir()))
docs/spectrogram-phase-proposal.mdadded+109−0View file
@@ -0,0 +1,109 @@
1+# Extending the phase test beyond steady tones
2+
3+The current app asks whether two periodic tones with the same harmonic amplitudes but different harmonic phases sound the same. A general sound, a drum hit or a spoken word, has no fixed set of harmonics, so "keep the amplitudes, change the phases" needs a new definition. The natural candidate is the spectrogram: at each time and frequency there is a magnitude and a phase, and we would like to keep every magnitude and change the phases. This document explains why that cannot be done literally, describes a formulation that is exact and easy to implement, shows what it does to a click and to a tone, and proposes how to build it into the app. The figures were generated by `make_figures.py` in this directory and can be regenerated with `python3 docs/make_figures.py`.
4+
5+The document is organized as follows. Section 1 fixes terminology. Section 2 shows what goes wrong when phases are changed cell by cell. Section 3 introduces the filter bank view, which is the STFT with a hop of one sample, and Section 4 defines the proposed operation on it. Sections 5 and 6 work out what the operation means for constant and for time-varying phase offsets, because both have well studied perceptual consequences. Section 7 proposes the app design and Section 8 lists open questions.
6+
7+## 1. Terminology: the STFT, the spectrogram, and the inverse
8+
9+The *short-time Fourier transform* (STFT) multiplies the signal by a window centered at time t, takes the Fourier transform of the product, and repeats for a grid of window positions spaced by the *hop*. The windows overlap when the hop is shorter than the window. Each windowed transform gives one column of complex coefficients X(t, f). The *spectrogram* is the magnitude |X(t, f)|, so strictly speaking the spectrogram has no phase; the object we want to modify is the STFT. With Gaussian windows the STFT is also called the Gabor transform. Figure 1 shows the test signal used throughout, a click followed by a harmonic tone at 220 Hz with eight harmonics of amplitude 1/k, together with the Gaussian windows and the spectrogram.
10+
11+
12+
13+*Figure 1. Top: the test signal with Gaussian analysis windows of width sigma = 5 ms, spaced by the hop. Bottom: the spectrogram, the magnitude of the windowed Fourier coefficients.*
14+
15+The *inverse STFT* (ISTFT) transforms each column back to a windowed piece of signal and adds the pieces, a procedure called *overlap-add*. The loop "STFT, modify, ISTFT" is what the phase vocoder literature calls analysis-modification-synthesis. The key fact about the STFT is that it is redundant: with a hop of one quarter of the window length there are four times as many complex coefficients as there are real samples. Only a small subset of all possible coefficient arrays are actually the STFT of some signal. Such arrays are called *consistent*. When a modified array is not consistent, the ISTFT does not fail; it returns the signal whose STFT is closest in the least-squares sense (Griffin and Lim, 1984), and the STFT of that signal differs from the array we asked for.
16+
17+## 2. Why phases cannot be changed cell by cell
18+
19+Figure 2 makes the consistency problem concrete. We compute a hop-based STFT of the test signal (Hann window, 32 ms, 75% overlap), multiply every coefficient by a random unit phasor, leaving every magnitude untouched, and invert. The recomputed spectrogram of the result is not the original spectrogram. The click has become a burst of noise spread over the window length and the tone has become a noisy band with the right long-term average but a fluctuating envelope. In the right column we instead multiply each frequency bin by one phasor that is the same in every frame. The spectrogram then survives almost unchanged, and the only visible difference is that the click is slightly smeared in time.
20+
21+
22+
23+*Figure 2. Left: original. Middle: after replacing the phase of every STFT cell by a random value and inverting; the magnitudes were untouched before inversion, but the inversion changes them. Right: after applying one phase per frequency bin, constant across frames.*
24+
25+The reason is that overlapping windows see the same samples, so the phases in neighbouring columns are tied together. A steady sinusoid at frequency f advances its phase by exactly 2 pi f times the hop from one column to the next. Random phases break this relation, the overlapping pieces no longer add coherently, and the ISTFT averages them into something else. Indeed, for a sufficiently redundant STFT the magnitudes alone determine the signal up to a global sign, which is why algorithms such as Griffin-Lim can reconstruct audio from a spectrogram. In that strict sense there is no such thing as "the same spectrogram with different phases". What we can do is preserve the magnitudes at the resolution of the analysis, and this turns out to be the useful question anyway.
26+
27+## 3. The filter bank view: a hop of one sample
28+
29+Rather than thinking of the STFT as a grid of windowed transforms, we think of it as a bank of filters. Fix a frequency f_k and consider the STFT coefficient at that frequency as a function of the window position t. It equals the convolution of the signal with h_k(t) = w(t) exp(2 pi i f_k t), a bandpass filter with a Gaussian frequency response centered at f_k. We write z_k(t) for this complex *subband signal*. Evaluating it at every sample, that is, using a hop of one sample, gives a row of the STFT at full time resolution. Its magnitude |z_k(t)| is the subband envelope and its argument is the subband phase.
30+
31+The frequency responses of the bank overlap heavily (Figure 3, top). We normalize them so that they sum to one at every frequency. Then the original signal is recovered by simply adding the subbands, x(t) = Re sum_k z_k(t), with no overlap-add in time at all. (The negative frequencies are folded in by doubling the positive ones first, the usual analytic-signal convention.) In the figures the bank uses a spacing between centers equal to the Gaussian width sigma_f = 1 / (2 pi sigma), so a 5 ms window gives about 32 Hz spacing and about 250 bands up to the 8 kHz Nyquist frequency of the 16 kHz test signal.
32+
33+
34+
35+*Figure 3. Top: the frequency responses of the bands for sigma = 5 ms, and their sum. Below: three subbands of the test signal, one on the fundamental, one on the third harmonic, one between harmonics. The click appears in every band; each harmonic appears in the bands near its frequency.*
36+
37+This is a standard reformulation, sometimes called the filter bank summation view of the phase vocoder, and it dissolves the question of consistency between overlapping windows. There are no separate windows to keep consistent, only subband signals that are added.
38+
39+## 4. The proposed operation: rotate the phase of each subband
40+
41+Here we propose to define the phase-modified sound as
42+
43+ y(t) = Re sum_k exp(i phi_k(t)) z_k(t)
44+
45+where phi_k(t) is the phase offset applied to band k at time t. Every subband envelope |z_k(t)| is exactly unchanged and every subband phase is shifted by exactly the requested amount, at every sample. When phi = 0 we get y = x exactly. This is a direct generalization of what the app does now: for a steady tone whose harmonics are separated by more than the bandwidth, each harmonic lives in its own group of bands, and rotating those bands by a common angle rotates that harmonic's phase.
46+
47+Figure 4 shows the effect on the test signal for three window widths, using a random phase map that varies smoothly across bands (the reason for the smoothness is given in Section 5). The steady tone keeps its period and its amplitude and only changes waveform shape, which is the situation the current app tests. The click spreads out over roughly the window length: a fraction of a millisecond for sigma = 0.5 ms, tens of milliseconds for sigma = 20 ms. This asymmetry between steady and transient sounds is the central new element the extension brings in.
48+
49+
50+
51+*Figure 4. Left: the click before and after phase rotation, for three window widths. Right: three periods of the tone. The steady tone changes shape but not envelope; the click is smeared over the window.*
52+
53+Figure 5 shows what "same spectrogram" means in this setting. Analysed with the same 20 ms window that was used to modify it, the modified signal has essentially the original spectrogram (the click column is a little mottled because different bands were delayed by different amounts within the window). Analysed with a 2 ms window, the original click is a sharp line while the modified click is visibly spread in time. The magnitudes are preserved at the resolution of the analysis and not at finer resolution, and this is unavoidable for the reasons of Section 2.
54+
55+
56+
57+*Figure 5. Spectrograms of the original (top) and of the modified signal from the 20 ms bank (bottom), analysed at sigma = 20 ms (left) and at sigma = 2 ms (right).*
58+
59+## 5. Constant phase offsets are an all-pass filter
60+
61+When the offsets do not depend on time, the operation is linear and time-invariant, so it is a filter, and its frequency response is G(f) = sum_k exp(i phi_k) H_k(f) where H_k is the response of band k. Figure 6 plots the magnitude, phase, and *group delay* of G for the three maps of Figure 4. The group delay is minus the derivative of the phase with respect to angular frequency, tau(f) = -(1 / 2 pi) d phi / d f, and it measures how much the energy near frequency f is delayed by the filter. If the phase varies smoothly across the bands, the magnitude stays within about a decibel of unity and the filter is close to an all-pass filter: it preserves the long-term magnitude spectrum exactly the way the tone app preserves the harmonic amplitudes.
62+
63+
64+
65+*Figure 6. The filter implied by constant per-band phases, for the three window widths of Figure 4. The group delay grows with the window width. The green band marks the approximate audibility threshold for group delay peaks.*
66+
67+The group delay is what makes the click spread. The phase change from one band to the next, divided by the band spacing, sets the delay, and because the spacing scales as 1 / sigma, a random map with a given roughness has a group delay proportional to the window width. Whether a listener can hear this is a well studied question, because all-pass filters occur in loudspeaker crossovers and equalizers. Blauert and Laws (1978) measured thresholds for group delay peaks of about 0.5 to 2 ms depending on frequency, with the lowest thresholds near 2 kHz and with impulsive signals over headphones; steady signals are much less sensitive. The prediction for our experiment is therefore concrete: constant random phase maps should be inaudible on sustained sounds for any window, inaudible on transients for windows below a millisecond or so, and audible as a softening or smearing of transients for windows of several milliseconds and beyond. The tone results already obtained, including the audibility at low fundamentals, are the sustained-sound end of this picture.
68+
69+One requirement follows from the overlap of the bands: the phase map must vary smoothly from band to band. Figure 7 shows why. Where two neighbouring bands have very different phases, their contributions partly cancel at the frequencies they share, and the magnitude response acquires notches of many decibels. A harmonic of the tone that falls on such a notch is attenuated, which changes the amplitude spectrum and defeats the purpose of the test. Smoothing the random map over a handful of bands removes the notches while leaving the group delay large. (An independent random phase per band is the analogue, in this setting, of the inconsistent per-cell modification of Section 2.) The same effect bounds how much delay any map can impose. For a pure delay tau the phase advances by tau / sigma radians per band, and the Gaussian averaging over neighbouring bands reduces the magnitude to exp(-tau^2 / (2 sigma^2)), so the group delay that a map can produce without changing magnitudes is limited to a fraction of the window width sigma. This is the quantitative form of the statement that the window width sets the group delay scale.
70+
71+
72+
73+*Figure 7. Left: independent random phases per band produce a magnitude response with deep notches. Right: the same random values smoothed over six bands give a nearly flat magnitude response. Dotted lines mark the harmonics of the tone.*
74+
75+Since the constant case is a filter, it does not need the filter bank at all: one could compute the FFT of the whole clip, multiply by exp(i phi(f)) for a smooth phase function, and invert. The filter bank version is equivalent, and we keep it because it also handles the time-varying case and because it displays the subband envelopes that the modification preserves.
76+
77+## 6. Time-varying phase offsets are frequency modulation
78+
79+The app's cycling mode glides the phases from A to B and back. In the subband formulation this is phi_k(t) = s(t) delta_k with s a raised cosine. A phase that changes in time is indistinguishable from a change in frequency: the instantaneous frequency of band k is shifted by phi_k'(t) / (2 pi). Figure 8 plots this for a swing of pi radians at a cycle rate of 0.5 Hz, where the peak shift is about 0.8 Hz, and the peak shift as a function of cycle rate. Near 220 Hz the just noticeable difference in pitch is of the order of 1 Hz, so slow cycling is close to the threshold and fast cycling well above it. At 10 Hz the harmonics wobble by about 16 Hz, each by a different amount, so the tone is momentarily inharmonic and what is heard is vibrato and roughness rather than phase.
80+
81+
82+
83+*Figure 8. A raised-cosine phase offset (left), the frequency shift it implies (middle), and the peak shift as a function of cycle rate (right). The green band is roughly the pitch JND near 220 Hz.*
84+
85+This caveat applies to the existing cycling mode as much as to the extension. We therefore propose to keep phi_k(t) as a slow morph between two static maps, exactly as now, and to treat the cycle rate as something to keep low rather than as a free parameter to explore. Slowly varying offsets also keep the filter bank modification close to consistent, so the analysis of Section 5 applies at each instant.
86+
87+## 7. Proposed app design
88+
89+Figure 9 summarizes the processing. Tone A is the original sound and tone B is the same sound with a phase map applied; both come out of the same filter bank, so they have the same subband envelopes at the chosen resolution.
90+
91+
92+
93+*Figure 9. The proposed pipeline. A and B share every subband envelope and differ only in subband phase.*
94+
95+The controls we propose are the following.
96+
97+- **Source sound.** A synthesized click, a noise burst, a plucked string (Karplus-Strong), and the existing harmonic tone cover the interesting cases and need no assets. A file loaded through `decodeAudioData` and a short microphone recording cover speech and music.
98+- **Window width sigma.** A logarithmic slider from about 0.3 ms to 50 ms. This is the parameter that sets the group delay scale and hence the predicted audibility on transients.
99+- **Phase map.** A linear map phi_k = -2 pi f_k tau_0 is a pure delay and serves as a control that must be inaudible. A quadratic map is a chirp all-pass, the classic stimulus in the group delay literature. A random map smoothed over a chosen number of bands is the generic case. A per-harmonic map, for the tone source, reproduces the current app.
100+- **Display.** The waveforms of A and B, and their spectrograms recomputed at the analysis resolution, so the listener can see that the envelopes match while the waveforms differ. A plot of the group delay of the map, with the threshold band, states the prediction before listening.
101+- **Listening.** The existing play A then B, blind same/different test, and slow cycling.
102+
103+Implementation is offline and fits the app's existing render-then-play design. For a clip of n samples, compute one FFT of length about 2n, and for each band multiply by the normalized Gaussian response, inverse transform, rotate by exp(i phi_k(t)), and accumulate the real part. With 256 bands and a two-second clip at 48 kHz this is a few hundred FFTs of length 2^18, on the order of a second in JavaScript, and it can run in a worker. The subband envelopes needed for the display come out of the same loop. The blind test component currently takes ToneParams for A and B; it would instead take two pre-rendered buffers, and the rest of it is unchanged. For the cycling mode, a real-time version would use a hop-based STFT with 75% overlap in an AudioWorklet and per-frame phase rotation, which is equivalent to the filter bank when the maps vary slowly. We suggest doing the offline version first.
104+
105+## 8. Open questions and limitations
106+
107+The claim that A and B "have the same spectrogram" holds only at the analysis resolution, and Section 5 shows that the audible differences are exactly the ones visible at finer resolution. A more principled choice than a free window width would be a bank matched to auditory filters, with bandwidths that grow with frequency (an ERB-spaced gammatone bank). Under the hypothesis that the ear uses only subband envelopes at auditory resolution, phase rotation in such a bank should be inaudible for every sound, and finding an audible case would be informative. We suggest starting with the Gaussian bank because it is simpler and because the window width is then a single interpretable knob, and adding the auditory bank later.
108+
109+The band at zero frequency and the smoothing at the ends of the map need care, and the choice of band spacing relative to bandwidth trades computation against how smooth the map must be. The group delay thresholds quoted above were measured with headphones and impulsive signals; over loudspeakers in a room the thresholds are higher because the room adds its own group delay. Finally, the extension does not change the basic statistics of the blind test, but with transient stimuli it will be important to randomize the interval between sounds so that timing cues do not leak information.