# SPDX-License-Identifier: MIT # © 2026 DTMFRAME """Synthetic-signal demonstration for the article linear-phase-eq. Questions it answers with numbers (no product, no licence, no recorded material): D1 Same magnitude, different time: a peaking EQ (+6 dB, 100 Hz, Q 1) as an Audio EQ Cookbook biquad and as a linear-phase FIR with the same magnitude. Magnitude difference, group delay at 50, 100, 200 and 1000 Hz, latency. D2 Pre-ringing: a low cut at 80 Hz (Butterworth, 4th and 2nd order) as IIR and as linear-phase FIR, fed a synthetic kick (55 Hz decaying sine burst with a click) and an impulse. Energy before the sound starts and how early the output passes -60 dB re its peak. D3 A processed copy summed with an unprocessed copy (parallel processing, two mics on one source): y = x + HPF(x), 2nd-order high-pass at 120 Hz on the copy, IIR versus linear-phase with and without the dry path delayed by the FIR latency. D4 A 2-band split at 200 Hz summed back: Butterworth 2nd-order pair (both polarities), Linkwitz-Riley LR4 pair, complementary linear-phase FIRs. D5 Latency against low-frequency precision: the shortest FIR that keeps the D1 peaking EQ within 0.5 dB at 50, 100, 400 and 1000 Hz. It also writes the article's concept diagram, out/concept_phase_time.svg (figures.py; D1's group delays at the four probe frequencies, IIR against linear-phase, on one time axis). Run: ../.venv/bin/python demo.py (writes out/ and out/record.json) """ from __future__ import annotations import sys from pathlib import Path import numpy as np from scipy import optimize, signal HERE = Path(__file__).resolve().parent sys.path.insert(0, str(HERE.parent / "common")) import demo_common as dc # noqa: E402 import figures # noqa: E402 import matplotlib.pyplot as plt # noqa: E402 FS = dc.FS SEED = 20261012 OUT = HERE / "out" EVIDENCE_ID = "ev_msr_synth_linear_phase_eq_1" N_TAPS = 16383 # linear-phase FIRs in D1-D4 LATENCY = (N_TAPS - 1) // 2 # 8191 samples N_DESIGN = 1 << 16 # frequency-sampling grid of the FIR design (0.73 Hz) N_EVAL = 1 << 18 # magnitude comparison grid (0.18 Hz) F_EVAL = np.fft.rfftfreq(N_EVAL, 1 / FS) SEL = (F_EVAL >= 20) & (F_EVAL <= 20000) GD_STEP_HZ = 0.01 # group delay: central difference of the phase over +-0.01 Hz D1_PROBES = (50.0, 100.0, 200.0, 1000.0) D3_PROBES = (60.0, 120.0, 240.0, 1000.0) D5_LIMIT_DB = 0.5 # ---------------------------------------------------------------- filters def rbj_peaking(f0: float, gain_db: float, q: float) -> np.ndarray: """Audio EQ Cookbook peaking EQ as one second-order section (b0 b1 b2 1 a1 a2).""" A = 10 ** (gain_db / 40) w0 = 2 * np.pi * f0 / FS alpha = np.sin(w0) / (2 * q) b = np.array([1 + alpha * A, -2 * np.cos(w0), 1 - alpha * A]) a = np.array([1 + alpha / A, -2 * np.cos(w0), 1 - alpha / A]) return np.concatenate([b / a[0], a / a[0]])[None, :] def butter(order: int, fc: float, btype: str) -> np.ndarray: return signal.butter(order, fc, btype=btype, fs=FS, output="sos") def h_sos(sos: np.ndarray, f) -> np.ndarray: return signal.sosfreqz(sos, worN=np.atleast_1d(np.asarray(f, dtype=float)), fs=FS)[1] def h_fir(h: np.ndarray, f) -> np.ndarray: return signal.freqz(h, [1.0], worN=np.atleast_1d(np.asarray(f, dtype=float)), fs=FS)[1] def linear_phase_twin(mag, n_taps: int, window: str = "boxcar") -> np.ndarray: """Linear-phase FIR with a given magnitude (frequency sampling). mag(f) is sampled on N_DESIGN/2+1 points from 0 Hz to Nyquist, the inverse real FFT gives the zero-phase response, which is moved by (n_taps-1)/2 samples, cut to n_taps (rectangular window unless another window is named) and made exactly symmetric, so the group delay is (n_taps-1)/2 samples at every frequency. """ f = np.linspace(0.0, FS / 2, N_DESIGN // 2 + 1) h0 = np.fft.irfft(mag(f), N_DESIGN) return cut_twin(h0, n_taps, window) def cut_twin(h0: np.ndarray, n_taps: int, window: str = "boxcar") -> np.ndarray: h = np.roll(h0, (n_taps - 1) // 2)[:n_taps] if window != "boxcar": h = h * signal.windows.get_window(window, n_taps, fftbins=False) return 0.5 * (h + h[::-1]) def mag_of(sos: np.ndarray): return lambda f: np.abs(h_sos(sos, f)) def group_delay_ms(H, f) -> np.ndarray: """-d(phase)/d(omega) by a central difference over f +- 0.01 Hz (ratio of the two responses, no unwrapping).""" f = np.atleast_1d(np.asarray(f, dtype=float)) r = H(f + GD_STEP_HZ) / H(f - GD_STEP_HZ) return -np.angle(r) / (2 * np.pi * 2 * GD_STEP_HZ) * 1000.0 def db(x) -> np.ndarray: return 20 * np.log10(np.abs(x)) def eval_db_fir(h: np.ndarray) -> np.ndarray: return db(np.fft.rfft(h, N_EVAL)[SEL]) def eval_db_sos(sos: np.ndarray) -> np.ndarray: return db(h_sos(sos, F_EVAL[SEL])) def max_abs(d: np.ndarray, f: np.ndarray) -> dict: j = int(np.argmax(np.abs(d))) return {"db": float(np.abs(d[j])), "at_hz": float(f[j])} def log_grid(lo: float, hi: float, per_octave: int, extra=()) -> np.ndarray: n = int(round(np.log2(hi / lo) * per_octave)) + 1 g = np.concatenate([lo * 2.0 ** (np.arange(n) / per_octave), np.array(extra, dtype=float)]) return np.unique(np.round(g[(g >= lo) & (g <= hi)], 3)) # ---------------------------------------------------------------- test material def kick(rng: np.random.Generator, hits_s=(0.25,), total_s: float = 1.0) -> tuple[np.ndarray, list[int]]: """Synthetic kick: 55 Hz sine from phase 0 decaying with a 80 ms time constant (400 ms, 10 ms fade-out), plus a 3 ms click of Gaussian noise decaying with 0.5 ms; peak normalised to 1.""" nb = int(0.4 * FS) t = np.arange(nb) / FS body = np.sin(2 * np.pi * 55.0 * t) * np.exp(-t / 0.080) fade = int(0.01 * FS) body[-fade:] *= np.linspace(1.0, 0.0, fade) nc = int(0.003 * FS) click = rng.standard_normal(nc) * np.exp(-np.arange(nc) / FS / 0.0005) * 0.5 one = body.copy() one[:nc] += click one /= np.max(np.abs(one)) x = np.zeros(int(round(total_s * FS))) onsets = [] for s in hits_s: n0 = int(round(s * FS)) x[n0:n0 + nb] += one onsets.append(n0) return x, onsets def pre_ringing(y: np.ndarray, n_ref: int) -> dict: """Energy before a reference sample (dB re total) and how early |y| first reaches -60 dB re its peak. Reported twice: against the output's main peak (max |y|) and against the onset of the input (the sample where the input starts, moved by the FIR latency for the linear-phase output). Lead times are positive when the output moves before the reference. """ e = y ** 2 total = float(np.sum(e)) pk = int(np.argmax(np.abs(y))) first = int(np.argmax(np.abs(y) >= np.abs(y[pk]) * 1e-3)) def ratio_db(n): v = float(np.sum(e[:n])) / total return float("-inf") if v == 0.0 else 10 * np.log10(v) return {"peak_after_onset_ms": (pk - n_ref) / FS * 1000, "energy_before_onset_db": ratio_db(n_ref), "minus60db_lead_before_onset_ms": (n_ref - first) / FS * 1000, "energy_before_main_peak_db": ratio_db(pk), "minus60db_lead_before_main_peak_ms": (pk - first) / FS * 1000} def envelope_db(y: np.ndarray, n_ref: int, lo_ms: float, hi_ms: float, bin_ms: float = 0.5, floor_db: float = -120.0): """Peak |y| in bins of bin_ms from lo_ms to hi_ms around n_ref, dB re max |y|, floored.""" nb = int(round(bin_ms * FS / 1000)) a = int(round(lo_ms * FS / 1000)) + n_ref n_bins = int(round((hi_ms - lo_ms) / bin_ms)) seg = np.abs(y[a:a + n_bins * nb]).reshape(n_bins, nb).max(axis=1) v = dc.db(seg / np.max(np.abs(y)), floor=10 ** (floor_db / 20)) return lo_ms + np.arange(n_bins) * bin_ms, np.maximum(v, floor_db) # ---------------------------------------------------------------- main def main() -> None: dc.clear_outputs(HERE, OUT) dc.setup_style() rng = np.random.default_rng(SEED) results: dict = {} f_sel = F_EVAL[SEL] # ------------------------------------------------------------ D1 same magnitude, different time sos_pk = rbj_peaking(100.0, 6.0, 1.0) h_pk = linear_phase_twin(mag_of(sos_pk), N_TAPS) H_iir = lambda f: h_sos(sos_pk, f) # noqa: E731 H_lin = lambda f: h_fir(h_pk, f) # noqa: E731 diff = eval_db_fir(h_pk) - eval_db_sos(sos_pk) gd_iir = group_delay_ms(H_iir, D1_PROBES) gd_lin = group_delay_ms(H_lin, D1_PROBES) gd_scipy = signal.group_delay((sos_pk[0, :3], sos_pk[0, 3:]), w=np.array(D1_PROBES), fs=FS)[1] / FS * 1000 g1 = log_grid(20, 20000, 48, D1_PROBES) gd_iir_curve = group_delay_ms(H_iir, g1) gd_lin_curve = group_delay_ms(H_lin, g1) i_max, i_min = int(np.argmax(gd_iir_curve)), int(np.argmin(gd_iir_curve)) results["D1_same_magnitude_different_time"] = { "iir": "Audio EQ Cookbook peaking EQ, 100 Hz, +6 dB, Q 1 (one biquad, scipy.signal.sosfilt)", "linear_phase": {"n_taps": N_TAPS, "design": "frequency sampling of the biquad magnitude, rectangular window, exactly symmetric", "symmetry_max_abs_difference": float(np.max(np.abs(h_pk - h_pk[::-1])))}, "max_magnitude_difference_db_20_20k": max_abs(diff, f_sel), "group_delay_ms": {"iir": {f"{f:g}": float(v) for f, v in zip(D1_PROBES, gd_iir)}, "linear_phase": {f"{f:g}": float(v) for f, v in zip(D1_PROBES, gd_lin)}}, "group_delay_samples_linear_phase": {f"{f:g}": float(v * FS / 1000) for f, v in zip(D1_PROBES, gd_lin)}, "iir_group_delay_extremes_20_20k": {"max_ms": float(gd_iir_curve[i_max]), "max_at_hz": float(g1[i_max]), "min_ms": float(gd_iir_curve[i_min]), "min_at_hz": float(g1[i_min])}, "linear_phase_group_delay_spread_20_20k_samples": float((gd_lin_curve.max() - gd_lin_curve.min()) * FS / 1000), "latency": {"iir_samples": 0, "iir_ms": 0.0, "linear_phase_samples": LATENCY, "linear_phase_ms": LATENCY / FS * 1000}, "check_iir_group_delay_vs_scipy_group_delay_max_abs_ms": float(np.max(np.abs(gd_iir - gd_scipy))), "note": "negative group delay: the envelope of a narrow-band component comes out earlier than without the filter (no output before the input)", } fig_gd = figures.iir_group_delay_ms(np.array(D1_PROBES)) assert np.max(np.abs(fig_gd - gd_iir)) < 1e-6 and figures.N_TAPS == N_TAPS, "figures.py disagrees with D1" g1m = log_grid(20, 20000, 48, D1_PROBES) dc.write_curve_csv(OUT / "d1_magnitude.csv", "frequency_response_csv", ["周波数 [Hz]", "IIR [dB]", "線形位相 [dB]"], [g1m, db(H_iir(g1m)), db(H_lin(g1m))], [3, 4, 4]) dc.write_curve_csv(OUT / "d1_group_delay.csv", "curve_csv", ["周波数 [Hz]", "IIR [ms]", "線形位相 [ms]"], [g1, gd_iir_curve, gd_lin_curve], [3, 4, 4]) fig, axes = plt.subplots(2, 1, figsize=(7.4, 5.6), sharex=True) ax = axes[0] ax.plot(g1m, db(H_iir(g1m)), color=dc.SERIES[0], label="IIR biquad") ax.plot(g1m, db(H_lin(g1m)), color=dc.SERIES[1], ls="--", label=f"linear-phase FIR, {N_TAPS} taps") ax.set_xscale("log") ax.set_xlim(20, 20000) dc.hz_axis(ax) ax.set_ylim(-1, 7) ax.set_ylabel("dB") ax.set_title("D1 Peaking +6 dB at 100 Hz, Q 1: the same magnitude") ax.legend(loc="upper right") ax = axes[1] ax.plot(g1, gd_iir_curve, color=dc.SERIES[0], label="IIR biquad") ax.axhline(0, color=dc.TEXT_2, lw=0.8) ax.set_ylim(-1.5, 2.8) ax.set_ylabel("group delay (ms)") ax.set_xlabel("frequency (Hz)") ax.set_title(f"Group delay: IIR varies with frequency; linear phase is {LATENCY / FS * 1000:.2f} ms everywhere") ax.legend(loc="upper right") fig.tight_layout() dc.save_figure(fig, OUT / "d1_group_delay.png") # ------------------------------------------------------------ D2 pre-ringing x_kick, (n_on,) = kick(rng) n_imp = 1 << 15 imp = np.zeros(n_imp) imp_at = 1 << 13 imp[imp_at] = 1.0 d2_rows, d2_filters, d2_out = [], {}, {} for order in (4, 2): sos = butter(order, 80.0, "highpass") h = linear_phase_twin(mag_of(sos), N_TAPS) dd = eval_db_fir(h) - eval_db_sos(sos) above = eval_db_sos(sos) >= -40 d2_filters[f"order_{order}"] = {"iir": f"Butterworth high-pass 80 Hz, order {order} (scipy.signal.butter, sos)", "max_magnitude_difference_db_20_20k": max_abs(dd, f_sel), "max_magnitude_difference_db_where_iir_above_minus40db": max_abs(dd[above], f_sel[above])} pad = np.zeros(FS) for inp, x, n0 in (("impulse", np.concatenate([imp, pad]), imp_at), ("kick", np.concatenate([x_kick, pad]), n_on)): y_lin = signal.fftconvolve(x, h) y_iir = signal.sosfilt(sos, x) for kind, y, ref in (("linear_phase", y_lin, n0 + LATENCY), ("iir", y_iir, n0)): d2_rows.append({"order": order, "filter": kind, "input": inp, **pre_ringing(y, ref)}) d2_out[(order, kind, inp)] = (y, ref) n_long = 2 * N_TAPS + 1 # 32767 taps, latency 16383 samples h_long = linear_phase_twin(mag_of(butter(4, 80.0, "highpass")), n_long) y_long = signal.fftconvolve(np.concatenate([x_kick, np.zeros(FS)]), h_long) check_long = pre_ringing(y_long, n_on + (n_long - 1) // 2) def row(order, kind, inp): return next(r for r in d2_rows if r["order"] == order and r["filter"] == kind and r["input"] == inp) results["D2_pre_ringing"] = { "filters": d2_filters, "linear_phase_n_taps": N_TAPS, "kick": "55 Hz sine from phase 0, decay time constant 80 ms, 400 ms with a 10 ms fade-out, plus a 3 ms click " "(seeded Gaussian noise, 0.5 ms decay, half amplitude); peak normalised to 1; onset at 0.25 s", "impulse": "unit impulse", "rows": d2_rows, "headline": {f"order_{o}": {"linear_phase_kick_energy_before_onset_db": row(o, "linear_phase", "kick")["energy_before_onset_db"], "linear_phase_kick_minus60db_lead_ms": row(o, "linear_phase", "kick")["minus60db_lead_before_onset_ms"], "linear_phase_impulse_energy_before_peak_db": row(o, "linear_phase", "impulse")["energy_before_main_peak_db"], "linear_phase_impulse_minus60db_lead_ms": row(o, "linear_phase", "impulse")["minus60db_lead_before_main_peak_ms"], "iir_kick_energy_before_onset_db": row(o, "iir", "kick")["energy_before_onset_db"], "iir_kick_minus60db_lead_ms": row(o, "iir", "kick")["minus60db_lead_before_onset_ms"]} for o in (4, 2)}, "check_order_4_kick_with_n_taps_32767": check_long, "note": "For the impulse the main peak is the onset. For the kick the main peak comes 0.1-0.4 ms after the onset, " "inside the click, so the energy before the main peak of the IIR output is the click's own start, not pre-ringing; " "the onset is the reference for pre-ringing.", } t_lo, t_hi, step = -20.0, 10.0, 2 y4l, r4l = d2_out[(4, "linear_phase", "kick")] y4i, r4i = d2_out[(4, "iir", "kick")] k = np.arange(int(round(t_lo * FS / 1000)), int(round(t_hi * FS / 1000)) + 1, step) dc.write_curve_csv(OUT / "d2_preringing.csv", "curve_csv", ["時間 [ms]", "線形位相 [-]", "IIR [-]"], [k / FS * 1000, y4l[r4l + k], y4i[r4i + k]], [4, 6, 6]) env = {} for o in (4, 2): for kind in ("linear_phase", "iir"): y, r = d2_out[(o, kind, "kick")] env[(o, kind)] = envelope_db(y, r, -60.0, 60.0) t_env = env[(4, "linear_phase")][0] dc.write_curve_csv(OUT / "d2_preringing_level.csv", "decay_curve_csv", ["時間 [ms]", "線形位相 4次 [dB]", "IIR 4次 [dB]", "線形位相 2次 [dB]", "IIR 2次 [dB]"], [t_env] + [env[(o, kd)][1] for o in (4, 2) for kd in ("linear_phase", "iir")], [1, 2, 2, 2, 2]) fig, axes = plt.subplots(2, 1, figsize=(7.4, 5.8)) ax = axes[0] tt = k / FS * 1000 ax.plot(tt, y4i[r4i + k], color=dc.SERIES[0], lw=1.2, label="IIR, 4th order") ax.plot(tt, y4l[r4l + k], color=dc.SERIES[1], lw=1.2, ls="--", label="linear phase, 4th order (latency removed)") ax.axvline(0, color=dc.TEXT_2, lw=0.8, ls=":") ax.set_xlim(t_lo, t_hi) ax.set_xlabel("time from the kick's onset (ms)") ax.set_ylabel("amplitude (input peak = 1)") ax.set_title("D2 Synthetic kick through an 80 Hz low cut") ax.legend(loc="upper left") ax = axes[1] for o, ls in ((4, "-"), (2, "--")): ax.plot(t_env, env[(o, "linear_phase")][1], color=dc.SERIES[1], ls=ls, lw=1.4, label=f"linear phase, order {o}") ax.plot(t_env, env[(o, "iir")][1], color=dc.SERIES[0], ls=ls, lw=1.4, label=f"IIR, order {o}") ax.axhline(-60, color=dc.TEXT_2, lw=0.8, ls=":") ax.set_xlim(-60, 60) ax.set_ylim(-122, 3) ax.set_xlabel("time from the kick's onset (ms)") ax.set_ylabel("peak level in 0.5 ms (dB re peak)") ax.legend(loc="lower right", ncol=2) fig.tight_layout() dc.save_figure(fig, OUT / "d2_preringing.png") # ------------------------------------------------------------ D3 parallel sum sos_hp = butter(2, 120.0, "highpass") h_hp = linear_phase_twin(mag_of(sos_hp), N_TAPS) d3_hp_diff = eval_db_fir(h_hp) - eval_db_sos(sos_hp) sums = { "iir": lambda f: 1.0 + h_sos(sos_hp, f), "linear_phase_compensated": lambda f: np.exp(-2j * np.pi * np.atleast_1d(np.asarray(f, dtype=float)) * LATENCY / FS) + h_fir(h_hp, f), "linear_phase_uncompensated": lambda f: 1.0 + h_fir(h_hp, f), } def lowest(H, lo, hi, n=20000): g = np.geomspace(lo, hi, n) v = db(H(g)) i = int(np.argmin(v)) a, b = g[max(i - 1, 0)], g[min(i + 1, n - 1)] r = optimize.minimize_scalar(lambda f: float(db(H(f))[0]), bounds=(a, b), method="bounded", options={"xatol": 1e-6}) return {"db": float(min(r.fun, v[i])), "at_hz": float(r.x if r.fun <= v[i] else g[i])} notch_spacing = FS / LATENCY fk = (np.arange(0, int(2000 / notch_spacing) + 2) + 0.5) * notch_spacing fk = fk[(fk >= 20) & (fk <= 2000)] vk = db(sums["linear_phase_uncompensated"](fk)) ik = int(np.argmin(vk)) rk = optimize.minimize_scalar(lambda f: float(db(sums["linear_phase_uncompensated"](f))[0]), bounds=(fk[ik] - notch_spacing / 4, fk[ik] + notch_spacing / 4), method="bounded", options={"xatol": 1e-9}) results["D3_parallel_sum"] = { "eq_on_copy": "Butterworth high-pass 120 Hz, order 2 (scipy.signal.butter, sos); linear-phase twin as in D1", "linear_phase_max_magnitude_difference_db_20_20k": max_abs(d3_hp_diff, f_sel), "cases": {"iir": "x + IIR(x)", "linear_phase_compensated": f"x delayed by {LATENCY} samples + FIR(x)", "linear_phase_uncompensated": "x + FIR(x), the dry path not delayed"}, "sum_db_re_one_copy": {k_: {f"{f:g}": float(v) for f, v in zip(D3_PROBES, db(H(np.array(D3_PROBES))))} for k_, H in sums.items()}, "lowest_point_20_2000": {"iir": lowest(sums["iir"], 20, 2000), "linear_phase_compensated": {**lowest(sums["linear_phase_compensated"], 20, 2000), "note": "the band edge: the compensated sum has no dip (1 + |H| >= 1)"}, "linear_phase_uncompensated": {"db": float(min(rk.fun, vk[ik])), "at_hz": float(rk.x if rk.fun <= vk[ik] else fk[ik])}}, "uncompensated_comb": {"delay_samples": LATENCY, "delay_ms": LATENCY / FS * 1000, "notch_spacing_hz": notch_spacing, "first_notch_hz": float(fk[0]), "notches_20_2000": int(len(fk)), "notches_below_minus20db_20_2000": int(np.sum(vk < -20)), "notch_frequencies": "(k + 1/2) fs / 8191"}, } g3 = np.round(np.arange(20.0, 500.0 + 1e-9, 0.5), 1) dc.write_curve_csv(OUT / "d3_parallel_sum.csv", "frequency_response_csv", ["周波数 [Hz]", "IIR [dB]", "線形位相 遅れを補正 [dB]", "線形位相 補正なし [dB]"], [g3] + [db(sums[k_](g3)) for k_ in ("iir", "linear_phase_compensated", "linear_phase_uncompensated")], [1, 4, 4, 4]) # The article's figure: the two compensated cases only, so the IIR dip near the corner is readable # (the uncompensated comb stays in d3_parallel_sum.csv and in results_detail.json). dc.write_curve_csv(OUT / "d3_parallel_sum_iir_vs_linear.csv", "frequency_response_csv", ["周波数 [Hz]", "IIR [dB]", "線形位相 遅れを補正 [dB]"], [g3] + [db(sums[k_](g3)) for k_ in ("iir", "linear_phase_compensated")], [1, 4, 4]) fig, ax = plt.subplots(figsize=(7.4, 3.8)) ax.plot(g3, db(sums["linear_phase_uncompensated"](g3)), color=dc.SERIES[2], lw=0.8, label="linear phase, dry not delayed") ax.plot(g3, db(sums["iir"](g3)), color=dc.SERIES[0], lw=2.0, label="IIR") ax.plot(g3, db(sums["linear_phase_compensated"](g3)), color=dc.SERIES[1], lw=2.0, ls="--", label="linear phase, dry delayed 8191 samples") ax.set_xscale("log") ax.set_xlim(20, 500) dc.hz_axis(ax) ax.set_ylim(-30, 8) ax.set_xlabel("frequency (Hz)") ax.set_ylabel("dB re one copy") ax.set_title("D3 Dry + high-passed copy (120 Hz, 2nd order)") ax.legend(loc="lower right") fig.tight_layout() dc.save_figure(fig, OUT / "d3_parallel_sum.png") # ------------------------------------------------------------ D4 crossover summed back lp2, hp2 = butter(2, 200.0, "lowpass"), butter(2, 200.0, "highpass") lr_lp, lr_hp = np.vstack([lp2, lp2]), np.vstack([hp2, hp2]) h_low = linear_phase_twin(mag_of(lr_lp), N_TAPS) h_high = -h_low.copy() h_high[LATENCY] += 1.0 pairs = { "butterworth2_in_phase": lambda f: h_sos(lp2, f) + h_sos(hp2, f), "butterworth2_high_inverted": lambda f: h_sos(lp2, f) - h_sos(hp2, f), "linkwitz_riley_lr4": lambda f: h_sos(lr_lp, f) + h_sos(lr_hp, f), "linear_phase_complementary": lambda f: h_fir(h_low, f) + h_fir(h_high, f), } g4d = np.unique(np.concatenate([np.geomspace(20, 20000, 20000), [200.0]])) d4_rows = [] unwrap_grid = np.linspace(1.0, 200.0, 20000) for name, H in pairs.items(): v = db(H(g4d)) j = int(np.argmax(np.abs(v))) rowd = {"pair": name, "max_abs_deviation_db_20_20k": float(abs(v[j])), "at_hz": float(g4d[j]), "sum_db_at_200hz": float(db(H([200.0]))[0])} if name == "butterworth2_in_phase": below = g4d[v < -3] rowd.update({"phase_deg_at_200hz": None, "group_delay_ms_at_200hz": None, "below_minus3db_hz": [float(below.min()), float(below.max())], "note": "the sum has a zero on the unit circle at 200 Hz (an exact null in exact arithmetic): " "the phase jumps by 180 degrees there and the group delay is undefined"}) else: lat = LATENCY if name == "linear_phase_complementary" else 0 ph = np.unwrap(np.angle(H(unwrap_grid) * np.exp(2j * np.pi * unwrap_grid * lat / FS))) rowd.update({"phase_deg_at_200hz": float(np.degrees(ph[-1])), "phase_latency_removed_samples": lat, "group_delay_ms_at_200hz": float(group_delay_ms(H, [200.0])[0])}) d4_rows.append(rowd) above60 = eval_db_sos(lr_hp) >= -60 hp_diff = eval_db_fir(h_high) - eval_db_sos(lr_hp) x_split = 0.5 * x_kick # -6 dBFS peak y_sum = signal.fftconvolve(x_split, h_low) + signal.fftconvolve(x_split, h_high) ref = np.zeros(len(y_sum)) ref[LATENCY:LATENCY + len(x_split)] = x_split null_peak = float(dc.db(np.max(np.abs(y_sum - ref)))) results["D4_crossover_sum"] = { "split_hz": 200.0, "pairs": {"butterworth2_in_phase": "Butterworth low-pass + high-pass, order 2 each, summed", "butterworth2_high_inverted": "the same with the high band's polarity inverted", "linkwitz_riley_lr4": "each band two Butterworth order-2 sections in cascade (LR4), summed", "linear_phase_complementary": f"low band: linear-phase FIR ({N_TAPS} taps) with the LR4 low-pass magnitude; " f"high band: a {LATENCY}-sample delay minus the low band"}, "rows": d4_rows, "fir_band_magnitude_vs_lr4": {"low_max_abs_db_20_20k": max_abs(eval_db_fir(h_low) - eval_db_sos(lr_lp), f_sel), "high_max_abs_db_where_lr4_above_minus60db": max_abs(hp_diff[above60], f_sel[above60])}, "linear_phase_split_and_sum_null_residual_peak_dbfs": null_peak, "null_test": "the D2 kick at -6 dBFS peak through both FIR bands separately (fftconvolve), summed, minus the input delayed by 8191 samples", } g4 = log_grid(20, 20000, 48) dc.write_curve_csv(OUT / "d4_crossover_sum.csv", "frequency_response_csv", ["周波数 [Hz]", "バターワース2次 [dB]", "バターワース2次 高域を逆相 [dB]", "LR4 [dB]", "線形位相FIR [dB]"], [g4] + [db(pairs[p](g4)) for p in pairs], [3, 4, 4, 4, 4]) g4p = np.unique(np.concatenate([np.linspace(1.0, 20.0, 200), g4])) phase_cols = [] for p, lat in (("butterworth2_high_inverted", 0), ("linkwitz_riley_lr4", 0), ("linear_phase_complementary", LATENCY)): ph = np.degrees(np.unwrap(np.angle(pairs[p](g4p) * np.exp(2j * np.pi * g4p * lat / FS)))) phase_cols.append(ph[g4p >= 20]) dc.write_curve_csv(OUT / "d4_crossover_phase.csv", "phase_response_csv", ["周波数 [Hz]", "バターワース2次 高域を逆相 [度]", "LR4 [度]", "線形位相FIR 遅れを除く [度]"], [g4p[g4p >= 20]] + phase_cols, [3, 3, 3, 3]) fig, axes = plt.subplots(2, 1, figsize=(7.4, 5.8), sharex=True) ax = axes[0] styles = (("Butterworth 2 + 2", dc.SERIES[2], "-"), ("Butterworth 2 + 2, high inverted", dc.SERIES[3], ":"), ("LR4", dc.SERIES[0], "-"), ("linear-phase FIR pair", dc.SERIES[1], "--")) for (lab, col, ls), p in zip(styles, pairs): ax.plot(g4, db(pairs[p](g4)), color=col, ls=ls, label=lab) ax.set_xscale("log") ax.set_xlim(20, 20000) dc.hz_axis(ax) ax.set_ylim(-30, 6) ax.set_ylabel("sum (dB)") ax.set_title("D4 Two bands split at 200 Hz, summed back") ax.legend(loc="lower right") ax = axes[1] for (lab, col, ls), ph in zip((styles[1], styles[2], styles[3]), phase_cols): ax.plot(g4p[g4p >= 20], ph, color=col, ls=ls, label=lab) ax.set_ylabel("phase (deg), FIR latency removed") ax.set_xlabel("frequency (Hz)") ax.set_ylim(-380, 20) fig.tight_layout() dc.save_figure(fig, OUT / "d4_crossover_sum.png") # ------------------------------------------------------------ D5 latency against precision d5_rows = [] windows = ("boxcar", "hann", "blackmanharris") for f0 in (50.0, 100.0, 400.0, 1000.0): sos = rbj_peaking(f0, 6.0, 1.0) target = eval_db_sos(sos) f = np.linspace(0.0, FS / 2, N_DESIGN // 2 + 1) h0 = np.fft.irfft(mag_of(sos)(f), N_DESIGN) for win in windows: errs = {} first = None n = 3 while first is None or n <= 2 * first: errs[n] = float(np.max(np.abs(eval_db_fir(cut_twin(h0, n, win)) - target))) if first is None and errs[n] <= D5_LIMIT_DB: first = n n += 2 longer_fail = [m for m, e in errs.items() if m > first and e > D5_LIMIT_DB] d5_rows.append({"f0_hz": f0, "window": win, "n_taps": first, "latency_samples": (first - 1) // 2, "latency_ms": (first - 1) / 2 / FS * 1000, "max_error_db": errs[first], "max_error_db_at_n_minus_2": errs[first - 2], "every_longer_n_passes_up_to": max(errs) if not longer_fail else None, "longer_n_that_fail": longer_fail}) by = {(r["f0_hz"], r["window"]): r for r in d5_rows} results["D5_latency_vs_precision"] = { "eq": "Audio EQ Cookbook peaking, +6 dB, Q 1, at each f0", "criterion": "max |magnitude difference| <= 0.5 dB over 20 Hz-20 kHz (0.18 Hz grid), shortest odd N (type I linear phase)", "search": "every odd N from 3 upward; the first N that passes is checked against every odd N up to twice its length", "main_design": "boxcar (rectangular window), the same design as D1-D4; hann and blackmanharris show how much the window choice moves N", "rows": d5_rows, "d1_n_taps_for_comparison": N_TAPS, } f5 = np.array([50.0, 100.0, 400.0, 1000.0]) dc.write_curve_csv(OUT / "d5_latency.csv", "curve_csv", ["持ち上げる周波数 [Hz]", "矩形窓 [ms]", "ハン窓 [ms]", "ブラックマン・ハリス窓 [ms]"], [f5] + [np.array([by[(f, w)]["latency_ms"] for f in f5]) for w in windows], [0, 4, 4, 4]) fig, ax = plt.subplots(figsize=(6.4, 3.8)) for (w, lab, col, ls) in (("boxcar", "rectangular window (D1-D4 design)", dc.SERIES[0], "-"), ("hann", "Hann window", dc.SERIES[1], "--"), ("blackmanharris", "Blackman-Harris window", dc.SERIES[2], ":")): ys = [by[(f, w)]["latency_ms"] for f in f5] ax.plot(f5, ys, color=col, ls=ls, marker="o", ms=6, label=lab) for f in f5: r = by[(f, "boxcar")] ax.annotate(f"N {r['n_taps']}", (f, r["latency_ms"]), textcoords="offset points", xytext=(0, -14), ha="center", fontsize=8, color=dc.TEXT_2) ax.set_xscale("log") ax.set_yscale("log") ax.set_xticks(f5, ["50", "100", "400", "1k"]) ax.minorticks_off() ax.set_xlabel("peaking EQ frequency (Hz), +6 dB, Q 1") ax.set_ylabel("latency (ms)") ax.set_title("D5 Shortest linear-phase FIR within 0.5 dB") ax.legend(loc="upper right") fig.tight_layout() dc.save_figure(fig, OUT / "d5_latency.png") # ------------------------------------------------------------ sound examples (not committed) rng_s = np.random.default_rng(SEED) x4, _ = kick(rng_s, hits_s=(0.25, 0.75, 1.25, 1.75), total_s=2.5) sos4 = butter(4, 80.0, "highpass") h4 = linear_phase_twin(mag_of(sos4), N_TAPS) clips = {"s1_kick_dry.wav": x4, "s2_kick_lowcut80_iir.wav": signal.sosfilt(sos4, x4), "s3_kick_lowcut80_linear_phase.wav": signal.fftconvolve(x4, h4)[LATENCY:LATENCY + len(x4)]} peak = max(np.max(np.abs(c)) for c in clips.values()) common_gain = 10 ** (-1.0 / 20) / peak for name, c in clips.items(): dc.write_wav(OUT / name, c * common_gain) results["sound_examples"] = {"files": list(clips), "committed": False, "common_gain_db": float(20 * np.log10(common_gain)), "note": "four D2 kicks 0.5 s apart; 4th-order 80 Hz low cut; the linear-phase file has its 8191-sample latency removed"} detail = { "evidence_id": EVIDENCE_ID, "article_slug": "linear-phase-eq", "signal_seed": SEED, "sample_rate_hz": FS, "common_module_sha256": dc.sha256_file(Path(dc.__file__)), "figures_module_sha256": dc.sha256_file(HERE / "figures.py"), "svg_module_sha256": dc.sha256_file(HERE.parent / "common" / "svgkit.py"), "method": { "arithmetic": "float64 throughout", "iir": "Audio EQ Cookbook peaking biquad; Butterworth and LR4 sections from scipy.signal.butter (bilinear, cutoff prewarped); applied with scipy.signal.sosfilt", "linear_phase_fir": f"frequency sampling: the IIR magnitude on {N_DESIGN // 2 + 1} points from 0 Hz to Nyquist, inverse real FFT ({N_DESIGN} points) " f"to the zero-phase response, moved by (N-1)/2 samples and cut to N taps (rectangular window), made exactly symmetric; " f"N = {N_TAPS} in D1-D4 (latency {LATENCY} samples); applied with scipy.signal.fftconvolve", "magnitude_difference": f"20 log10 of the FIR's rfft ({N_EVAL} points, 0.18 Hz) against the IIR's response at the same frequencies, 20 Hz-20 kHz", "group_delay": "-d(phase)/d(omega) by a central difference of the complex response over f +- 0.01 Hz; checked against scipy.signal.group_delay for the D1 biquad", "pre_ringing": "energy before the reference sample divided by the total energy of the whole output (dB); lead = reference minus the first sample " "whose |y| reaches 1e-3 of max |y| (-60 dB), in ms; reference = the input's onset (plus the FIR latency for the linear-phase output) " "and, separately, the output's main peak", "pre_ringing_level_csv": "peak |y| in 0.5 ms bins from -60 to +60 ms around the onset, dB re max |y|; values below -120 dB are written as -120 dB", "sums": "transfer functions evaluated at exact frequencies (scipy.signal.sosfreqz, freqz); lowest points refined with scipy.optimize.minimize_scalar (bounded)", "phase": "unwrapped from 1 Hz; the linear-phase FIR's pure delay removed", }, "results": results, } r1 = results["D1_same_magnitude_different_time"] ledger = { "summary_editorial": "合成した信号で、IIRのEQと同じ振幅の線形位相FIRの、時間の振る舞いと遅れの違いを確かめた実証の記録です。", "subject_ja": "同じ振幅のIIRのEQと線形位相FIRで、群遅延、プリリンギング、並列の足し算、帯域分割の和、遅れと精度の関係を比べる実証です。", "material_description_ja": "スクリプトが式と乱数から生成したインパルスと、55 Hzの減衰する正弦波にクリックを足した合成のキックです。", "settings": [ {"name": "標本化周波数", "value": "48000 Hz、64ビット浮動小数"}, {"name": "線形位相FIRの作り方", "value": "IIRの振幅を周波数で標本化し、逆FFTで零位相の応答にして切り出す（矩形窓）"}, {"name": "線形位相FIRの長さ", "value": f"{N_TAPS}タップ、遅れ{LATENCY}サンプル（{LATENCY / FS * 1000:.2f} ms）"}, {"name": "群遅延を比べたEQ", "value": "100 Hz、+6 dB、Q 1のピーキング（Audio EQ Cookbookの双2次フィルター）"}, {"name": "プリリンギングを見たローカット", "value": "80 Hzのバターワースのハイパス、4次と2次"}, {"name": "プリリンギングを見た素材", "value": "55 Hzの減衰する正弦波とクリックの合成のキック、インパルス"}, {"name": "並列の足し算", "value": "元の音に、120 Hzの2次のハイパスをかけた複製を足す"}, {"name": "帯域分割", "value": "200 Hz、バターワース2次の組、LR4の組、相補の線形位相FIR"}, {"name": "最短のタップ数の基準", "value": "20 Hz〜20 kHzで振幅の差が0.5 dB以内（ピーキング、+6 dB、Q 1）"}, ], "comparison_conditions_ja": ["同じ振幅の特性を持つIIRのフィルターと線形位相FIRを、同じ入力と同じ計算の方法で比べました。"], "procedure_ja": [ "IIRのフィルターの振幅を周波数で標本化し、逆FFTで零位相の応答にして、同じ振幅の線形位相FIRを作りました。", "100 Hzを+6 dB持ち上げるEQで、両者の振幅の差と、50 Hz、100 Hz、200 Hz、1000 Hzの群遅延を求めました。", "80 Hzのローカットに合成のキックとインパルスを通し、音の始まりより前に出るエネルギーと、-60 dBに達する時刻を求めました。", "元の音と、120 Hzのハイパスをかけた複製を足した特性を、遅れの補正の有無で比べました。", "200 Hzで2つの帯域に分けて足し戻した特性と位相を、フィルターの組ごとに求めました。", "持ち上げる周波数ごとに、振幅の差が0.5 dB以内になる最短のFIRの長さと、その遅れを求めました。", ], "results": { "measured_latency_samples": {"value": round(r1["group_delay_samples_linear_phase"]["100"], 0), "decimals": 0}, "null_residual_peak_dbfs": {"value": round(null_peak, 2), "decimals": 2}, }, "signal_seed": SEED, } csv_kinds = {"d1_magnitude.csv": "frequency_response_csv", "d1_group_delay.csv": "curve_csv", "d2_preringing.csv": "curve_csv", "d2_preringing_level.csv": "decay_curve_csv", "d3_parallel_sum.csv": "frequency_response_csv", "d3_parallel_sum_iir_vs_linear.csv": "frequency_response_csv", "d4_crossover_sum.csv": "frequency_response_csv", "d4_crossover_phase.csv": "phase_response_csv", "d5_latency.csv": "curve_csv"} figures.write_all(OUT) rec = dc.write_outputs(HERE, OUT, EVIDENCE_ID, csv_kinds, detail, ledger) print(rec) if __name__ == "__main__": main()