#!/usr/bin/env python3 """Generate the figures for spectrogram-phase-proposal.md. Run from anywhere: python3 docs/make_figures.py Figures are written to docs/figures/. """ from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np from scipy.ndimage import gaussian_filter1d from matplotlib.patches import FancyArrowPatch, FancyBboxPatch OUT = Path(__file__).parent / "figures" OUT.mkdir(exist_ok=True) rng = np.random.default_rng(3) FS = 16000 T = 0.8 N = int(FS * T) t = np.arange(N) / FS F0 = 220 CLICK_T = 0.15 TONE_T0, TONE_T1 = 0.35, 0.7 COL_A = "#1f5fa8" COL_B = "#d9641e" GRAY = "#888888" plt.rcParams.update({"font.size": 9, "axes.titlesize": 10, "figure.dpi": 130}) # ---------------------------------------------------------------------------- # Test signal: a click followed by a harmonic tone # ---------------------------------------------------------------------------- def test_signal(): x = np.zeros(N) i0 = int(CLICK_T * FS) x[i0 : i0 + 4] += np.hanning(6)[1:-1] seg = (t >= TONE_T0) & (t < TONE_T1) tt = t[seg] - TONE_T0 env = np.ones_like(tt) fade = int(0.01 * FS) env[:fade] = 0.5 - 0.5 * np.cos(np.pi * np.arange(fade) / fade) env[-fade:] = env[:fade][::-1] tone = sum(np.sin(2 * np.pi * k * F0 * tt) / k for k in range(1, 9)) x[seg] += 0.25 * env * tone return x x = test_signal() # ---------------------------------------------------------------------------- # Gaussian filter bank (the STFT with a hop of one sample) # ---------------------------------------------------------------------------- class GaussBank: """Complex Gaussian bandpass filters whose responses sum to one on f >= 0. Band k has center f_k = k * df and response W(f - f_k) / S(f), where W is a Gaussian of width sigma_f = 1 / (2 pi sigma_t) and S(f) = sum_k W(f - f_k). The subband z_k = x * h_k is the STFT at bin k evaluated at every sample. Because the responses sum to one, x = Re sum_k z_k exactly (positive frequencies are doubled first). """ def __init__(self, n, fs, sigma_t, spacing=1.0): self.n, self.fs, self.sigma_t = n, fs, sigma_t pad = int(8 * sigma_t * fs) self.L = 1 << int(np.ceil(np.log2(n + 2 * pad))) self.f = np.fft.fftfreq(self.L, 1 / fs) self.sigma_f = 1 / (2 * np.pi * sigma_t) self.df = spacing * self.sigma_f self.centers = np.arange(0, fs / 2 + self.df / 2, self.df) S = np.zeros(self.L) for fc in self.centers: S += self._W(fc) self.S = np.where(S > 1e-300, S, np.inf) # far outside the covered range the sum underflows; those bins are never used def _W(self, fc): return np.exp(-((self.f - fc) ** 2) / (2 * self.sigma_f**2)) def response(self, k): return self._W(self.centers[k]) / self.S def analytic_spectrum(self, x): X = np.fft.fft(x, self.L) X[self.f < 0] = 0 X[self.f > 0] *= 2 return X def subband(self, Xa, k): return np.fft.ifft(Xa * self.response(k))[: self.n] def apply(self, x, phase=None): """y = Re sum_k exp(i phi_k(t)) z_k(t). phase(k) returns a scalar or an array over time.""" Xa = self.analytic_spectrum(x) y = np.zeros(self.n) for k in range(len(self.centers)): z = self.subband(Xa, k) if phase is not None: z = z * np.exp(1j * phase(k)) y += z.real return y def spectrogram(self, x, hop, fmax): Xa = self.analytic_spectrum(x) ks = [k for k, fc in enumerate(self.centers) if fc <= fmax] idx = np.arange(0, self.n, hop) M = np.zeros((len(ks), len(idx))) for i, k in enumerate(ks): M[i] = np.abs(self.subband(Xa, k)[idx]) return M, self.centers[ks], idx / self.fs def effective_filter(self, phi): """Frequency response of the LTI filter obtained with constant per-band phases phi_k.""" G = np.zeros(self.L, dtype=complex) for k, p in enumerate(phi): G += np.exp(1j * p) * self.response(k) return G def smooth_random_phases(nbands, smooth_bands, rng): """Uniform random phases, smoothed across bands with a Gaussian of the given width, then rescaled.""" raw = rng.uniform(-np.pi, np.pi, nbands) if smooth_bands <= 0: return raw sm = gaussian_filter1d(raw, smooth_bands, mode="reflect") return sm * raw.std() / sm.std() def db(M): return 20 * np.log10(np.maximum(M, 1e-12) / M.max()) def show_spec(ax, M, freqs, times, title, vmin=-60): ax.imshow(db(M), origin="lower", aspect="auto", extent=[times[0], times[-1], freqs[0], freqs[-1]], vmin=vmin, vmax=0, cmap="magma") ax.set_title(title) ax.set_xlabel("time (s)") ax.set_ylabel("frequency (Hz)") # ---------------------------------------------------------------------------- # Figure 1: the STFT as sliding windows, and the spectrogram # ---------------------------------------------------------------------------- def fig_stft_basics(): sigma = 0.005 fig, axs = plt.subplots(2, 1, figsize=(8, 5.2), gridspec_kw={"height_ratios": [1, 1.4]}) ax = axs[0] ax.plot(t, x, color="k", lw=0.6) hop = 0.02 for i, c in enumerate(np.arange(0.06, 0.56, hop)): w = np.exp(-((t - c) ** 2) / (2 * sigma**2)) ax.plot(t, 1.05 * w, color=COL_A if i % 2 else COL_B, lw=0.9, alpha=0.8) ax.annotate("hop", xy=(0.30, 1.12), xytext=(0.30, 1.35), ha="center", arrowprops=dict(arrowstyle="-[, widthB=1.6", lw=0.8)) ax.set_xlim(0, T) ax.set_ylim(-1.1, 1.5) ax.set_ylabel("amplitude") ax.set_title("Test signal (a click, then a harmonic tone) with Gaussian analysis windows, sigma = 5 ms") bank = GaussBank(N, FS, sigma) M, freqs, times = bank.spectrogram(x, 32, 4000) show_spec(axs[1], M, freqs, times, "Spectrogram: magnitude of each windowed Fourier coefficient (dB)") fig.tight_layout() fig.savefig(OUT / "01-stft-basics.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 2: the filter bank view # ---------------------------------------------------------------------------- def fig_filterbank(): sigma = 0.005 bank = GaussBank(N, FS, sigma) fig = plt.figure(figsize=(8, 7.5)) gs = fig.add_gridspec(4, 1, height_ratios=[1, 1, 1, 1]) ax = fig.add_subplot(gs[0]) pos = bank.f >= 0 f = bank.f[pos] total = np.zeros(pos.sum()) for k in range(len(bank.centers)): r = bank.response(k)[pos] total += r if bank.centers[k] <= 1500: ax.plot(f, r, lw=0.8, color=COL_A if k % 2 else COL_B, alpha=0.8) ax.plot(f, total, "k--", lw=1, label="sum over all bands") ax.set_xlim(0, 1500) ax.set_ylim(0, 1.15) ax.set_xlabel("frequency (Hz)") ax.set_ylabel("gain") ax.legend(loc="lower right") ax.set_title("Each STFT bin is a bandpass filter. The filters overlap and sum to one.") Xa = bank.analytic_spectrum(x) picks = [int(round(F0 * m / bank.df)) for m in (1, 3)] + [int(round(1200 / bank.df))] for row, k in enumerate(picks, start=1): ax = fig.add_subplot(gs[row]) z = bank.subband(Xa, k) ax.plot(t, z.real, color=COL_A, lw=0.6, label="Re z_k(t)") ax.plot(t, np.abs(z), color=COL_B, lw=1.0, label="|z_k(t)| (envelope)") ax.set_xlim(0, T) ax.set_ylabel("amplitude") ax.set_title(f"Subband k = {k}, center {bank.centers[k]:.0f} Hz", loc="left") if row == 1: ax.legend(loc="upper right", ncol=2) ax.set_xlabel("time (s)") y = bank.apply(x) 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) fig.tight_layout(rect=(0, 0.02, 1, 1)) fig.savefig(OUT / "02-filterbank.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 3: rotating subband phases, for three window lengths # ---------------------------------------------------------------------------- SMOOTH = 6 # bands def rotated_versions(sigmas, smooth_bands=SMOOTH, seed=5): out = {} for s in sigmas: bank = GaussBank(N, FS, s) phi = smooth_random_phases(len(bank.centers), smooth_bands, np.random.default_rng(seed)) y = bank.apply(x, lambda k: phi[k]) out[s] = (bank, phi, y) return out ROT = rotated_versions([0.0005, 0.005, 0.02]) def fig_phase_rotation(): fig, axs = plt.subplots(3, 2, figsize=(8, 7), gridspec_kw={"width_ratios": [1.4, 1]}) for row, (s, (bank, phi, y)) in enumerate(ROT.items()): ax = axs[row, 0] sel = (t > CLICK_T - 0.05) & (t < CLICK_T + 0.08) ax.plot(t[sel] * 1000, x[sel], color=GRAY, lw=0.8, label="original") ax.plot(t[sel] * 1000, y[sel], color=COL_B, lw=0.8, label="phases rotated") ax.set_title(f"window sigma = {s*1000:g} ms: the click", loc="left") ax.set_ylabel("amplitude") if row == 0: ax.legend(loc="upper right") if row == 2: ax.set_xlabel("time (ms)") ax = axs[row, 1] sel = (t > 0.5) & (t < 0.5 + 3 / F0) ax.plot(t[sel] * 1000, x[sel], color=GRAY, lw=0.8) ax.plot(t[sel] * 1000, y[sel], color=COL_B, lw=0.8) ax.set_title("same map, steady tone (3 periods)", loc="left") if row == 2: ax.set_xlabel("time (ms)") fig.suptitle("y = Re sum_k exp(i phi_k) z_k(t) with smooth random phi_k", fontsize=10) fig.tight_layout(rect=(0, 0, 1, 0.96)) fig.savefig(OUT / "03-phase-rotation.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 4: spectrograms before and after, at two analysis resolutions # ---------------------------------------------------------------------------- def fig_spectrograms(): _, _, y = ROT[0.02] fig, axs = plt.subplots(2, 2, figsize=(8, 5.6)) for col, s in enumerate([0.02, 0.002]): bank = GaussBank(N, FS, s) for row, (sig, name) in enumerate([(x, "original"), (y, "rotated (20 ms bank)")]): M, freqs, times = bank.spectrogram(sig, 32, 4000) show_spec(axs[row, col], M, freqs, times, f"{name}, analysed at sigma = {s*1000:g} ms") fig.tight_layout() fig.savefig(OUT / "04-spectrograms.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 5: the equivalent all-pass filter and its group delay # ---------------------------------------------------------------------------- def fig_group_delay(): fig, axs = plt.subplots(3, 1, figsize=(8, 7), sharex=True) for s, (bank, phi, _) in ROT.items(): G = bank.effective_filter(phi) pos = (bank.f > 0) & (bank.f <= 4000) f = bank.f[pos] ph = np.unwrap(np.angle(G[pos])) tau = -np.gradient(ph, f) / (2 * np.pi) lab = f"sigma = {s*1000:g} ms" axs[0].plot(f, 20 * np.log10(np.abs(G[pos])), lw=0.8, label=lab) axs[1].plot(f, ph, lw=0.8, label=lab) axs[2].plot(f, tau * 1000, lw=0.8, label=lab) axs[0].set_ylabel("magnitude (dB)") axs[0].set_ylim(-3, 1) axs[0].set_title("Constant per-band phases act as a filter G(f): magnitude (top), phase (middle), group delay (bottom)") axs[0].legend(loc="lower left", ncol=3, fontsize=8) axs[1].set_ylabel("phase (rad)") axs[2].axhspan(-1.5, 1.5, color="green", alpha=0.12) 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") axs[2].set_ylabel("group delay (ms)") axs[2].set_xlabel("frequency (Hz)") axs[2].set_xlim(0, 4000) fig.tight_layout() fig.savefig(OUT / "05-group-delay.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 6: why the phases must vary smoothly across neighbouring bands # ---------------------------------------------------------------------------- def fig_notches(): bank = GaussBank(N, FS, 0.005) fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex="col") for col, (sb, name) in enumerate([(0, "independent random phase per band"), (SMOOTH, f"random phases smoothed over {SMOOTH} bands")]): phi = smooth_random_phases(len(bank.centers), sb, np.random.default_rng(5)) G = bank.effective_filter(phi) pos = (bank.f > 0) & (bank.f <= 2000) f = bank.f[pos] ax = axs[0, col] ax.plot(bank.centers, np.angle(np.exp(1j * phi)), ".", ms=3, color=COL_B) ax.set_ylabel("phi_k (rad)") ax.set_title(name) ax.set_xlim(0, 2000) ax = axs[1, col] ax.plot(f, 20 * np.log10(np.abs(G[pos])), color=COL_A, lw=0.8) ax.set_ylim(-30, 3) ax.set_ylabel("|G(f)| (dB)") ax.set_xlabel("frequency (Hz)") for m in range(1, 9): ax.axvline(m * F0, color=GRAY, lw=0.5, ls=":") axs[1, 0].set_title("|G(f)|: notches where neighbours disagree", fontsize=9) axs[1, 1].set_title("|G(f)|: nearly flat (dotted: harmonics of the tone)", fontsize=9) fig.tight_layout() fig.savefig(OUT / "06-notches.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 7: hop-based STFT, per-cell random phase, inverse STFT # ---------------------------------------------------------------------------- def stft(x, n, hop): w = np.hanning(n + 1)[:-1] frames = [x[i : i + n] * w for i in range(0, len(x) - n + 1, hop)] return np.fft.rfft(np.array(frames), axis=1), w def istft(S, w, hop, length): n = len(w) y = np.zeros(length) norm = np.zeros(length) for m, spec in enumerate(S): i = m * hop y[i : i + n] += np.fft.irfft(spec, n) * w norm[i : i + n] += w * w return y / np.maximum(norm, 1e-9) def fig_inconsistency(): n, hop = 512, 128 S, w = stft(x, n, hop) freqs = np.fft.rfftfreq(n, 1 / FS) keep = freqs <= 4000 times = np.arange(S.shape[0]) * hop / FS S_cell = S * np.exp(1j * rng.uniform(-np.pi, np.pi, S.shape)) y_cell = istft(S_cell, w, hop, N) phi_band = smooth_random_phases(S.shape[1], SMOOTH, np.random.default_rng(5)) S_band = S * np.exp(1j * phi_band)[None, :] y_band = istft(S_band, w, hop, N) fig, axs = plt.subplots(2, 3, figsize=(9.5, 5.4), gridspec_kw={"height_ratios": [1, 1.6]}) for col, (sig, title) in enumerate([(x, "original"), (y_cell, "after random phase in every cell"), (y_band, "after one phase per frequency bin")]): axs[0, col].plot(t, sig, color="k", lw=0.5) axs[0, col].set_xlim(0, T) axs[0, col].set_ylim(-1.1, 1.1) axs[0, col].set_title(title) M = np.abs(stft(sig, n, hop)[0])[:, keep].T show_spec(axs[1, col], M, freqs[keep], times, "spectrogram recomputed from the signal") fig.suptitle("Hop-based STFT (Hann, 32 ms, 75% overlap), modify phases, inverse STFT by overlap-add", fontsize=10) fig.tight_layout(rect=(0, 0, 1, 0.96)) fig.savefig(OUT / "07-inconsistency.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 8: a time-varying phase offset is frequency modulation # ---------------------------------------------------------------------------- def fig_time_varying(): fc = 0.5 delta = np.pi tt = np.linspace(0, 2 / fc, 2000) s = 0.5 - 0.5 * np.cos(2 * np.pi * fc * tt) phi = delta * s dphi = np.gradient(phi, tt) fig, axs = plt.subplots(1, 3, figsize=(9.5, 2.9)) axs[0].plot(tt, phi, color=COL_A) axs[0].set_title("phi(t): 0 to pi and back at 0.5 Hz") axs[0].set_xlabel("time (s)") axs[0].set_ylabel("rad") axs[1].plot(tt, dphi / (2 * np.pi), color=COL_B) axs[1].set_title("frequency shift phi'(t) / 2 pi") axs[1].set_xlabel("time (s)") axs[1].set_ylabel("Hz") fcs = np.logspace(-2, 2, 100) axs[2].loglog(fcs, delta * fcs / 2, color=COL_B) axs[2].axhspan(0.5, 2, color="green", alpha=0.12) axs[2].text(0.011, 2.4, "roughly the pitch JND near 220 Hz", fontsize=8, color="green") axs[2].set_title("peak shift vs cycle rate (swing pi)") axs[2].set_xlabel("cycle rate (Hz)") axs[2].set_ylabel("peak shift (Hz)") fig.tight_layout() fig.savefig(OUT / "08-time-varying.png") plt.close(fig) # ---------------------------------------------------------------------------- # Figure 9: proposed processing pipeline # ---------------------------------------------------------------------------- def fig_pipeline(): fig, ax = plt.subplots(figsize=(9.5, 3.4)) ax.set_xlim(0, 10) ax.set_ylim(0, 3.4) ax.axis("off") def box(xc, yc, wdt, hgt, text, color="#eef3fa", ec="#333"): ax.add_patch(FancyBboxPatch((xc - wdt / 2, yc - hgt / 2), wdt, hgt, boxstyle="round,pad=0.05", fc=color, ec=ec, lw=0.9)) ax.text(xc, yc, text, ha="center", va="center", fontsize=8.5) def arrow(x0, y0, x1, y1): ax.add_patch(FancyArrowPatch((x0, y0), (x1, y1), arrowstyle="-|>", mutation_scale=12, lw=0.9, color="#333")) box(0.9, 1.7, 1.5, 1.2, "source sound x(t)\n(click, noise burst,\npluck, file, mic)") box(3.1, 1.7, 1.7, 1.2, "Gaussian filter bank\nwindow sigma\n\nz_k(t) for each band k") box(5.55, 2.6, 1.9, 0.75, "phase map A: phi_k = 0", color="#e5edf8", ec=COL_A) box(5.55, 0.8, 1.9, 0.75, "phase map B: phi_k chosen", color="#fbeadf", ec=COL_B) box(7.9, 2.6, 1.7, 0.75, "A = Re sum z_k = x", color="#e5edf8", ec=COL_A) box(7.9, 0.8, 1.7, 0.75, "B = Re sum e^(i phi_k) z_k", color="#fbeadf", ec=COL_B) box(9.4, 1.7, 1.1, 1.0, "blind test\nA/B, cycling") arrow(1.65, 1.7, 2.25, 1.7) arrow(3.95, 1.9, 4.6, 2.6) arrow(3.95, 1.5, 4.6, 0.8) arrow(6.5, 2.6, 7.05, 2.6) arrow(6.5, 0.8, 7.05, 0.8) arrow(8.75, 2.6, 9.0, 2.0) arrow(8.75, 0.8, 9.0, 1.4) ax.text(5.55, 1.7, "same subband envelopes |z_k(t)|,\ndifferent subband phases", ha="center", va="center", fontsize=8, color="#555") fig.tight_layout() fig.savefig(OUT / "09-pipeline.png") plt.close(fig) if __name__ == "__main__": fig_stft_basics() fig_filterbank() fig_phase_rotation() fig_spectrograms() fig_group_delay() fig_notches() fig_inconsistency() fig_time_varying() fig_pipeline() print("wrote", sorted(p.name for p in OUT.iterdir()))