# SPDX-License-Identifier: MIT # © 2026 DTMFRAME """Synthetic-signal demonstration for the article mid-side-processing. Claims it tests with numbers (no product, no licence, no recorded material): C1 Is M/S encode -> decode lossless? (null residual in float64, float32 and 24-bit storage, for the 1/2 and the 1/sqrt(2) conventions) C2 Does going through M/S change the level? (the four combinations of the two conventions) C3 What does raising S do to the mono fold-down and to the stereo/mono difference? (synthetic mix: centre noise, hard-left and hard-right noise; S gain -12..+6 dB) C4 What does an EQ applied to S alone do to a hard-left source? (peaking +6 dB at 3 kHz, Q 1, as an IIR biquad and as a linear-phase FIR with the same magnitude) C5 What does a high-pass on S (bass made mono) do above the corner? (100 Hz, Q 0.707, IIR biquad vs linear-phase) C6 What happens when the S path is late by a few samples and M is not? (S delayed by 32 samples, no EQ) It also writes the article's concept diagram, out/concept_ms_flow.svg (figures.py; the M/S round trip and the mono fold-down as a signal flow, no data). 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 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 = 20261019 OUT = HERE / "out" EVIDENCE_ID = "ev_msr_synth_mid_side_1" SQRT2 = np.sqrt(2.0) # ---------------------------------------------------------------- M/S matrices def encode(L, R, conv: str): if conv == "half": # M = (L+R)/2, S = (L-R)/2 return (L + R) / 2, (L - R) / 2 if conv == "sqrt2": # M = (L+R)/sqrt2, S = (L-R)/sqrt2 return (L + R) / SQRT2, (L - R) / SQRT2 raise ValueError(conv) def decode(M, S, conv: str): if conv == "half": # inverse of "half": L = M+S, R = M-S return M + S, M - S if conv == "sqrt2": # inverse of "sqrt2": L = (M+S)/sqrt2, R = (M-S)/sqrt2 return (M + S) / SQRT2, (M - S) / SQRT2 raise ValueError(conv) def dbfs(x) -> float: v = float(np.max(np.abs(x))) return float("-inf") if v == 0.0 else 20 * np.log10(v) def rms_dbfs(x) -> float: v = float(np.sqrt(np.mean(np.square(x)))) return float("-inf") if v == 0.0 else 20 * np.log10(v) def q24(x): """Store as 24-bit PCM (round to the nearest step of 2^-23, clip to the range).""" step = 2.0 ** -23 return np.clip(np.round(x / step), -2 ** 23, 2 ** 23 - 1) * step # ---------------------------------------------------------------- filters def rbj_peaking(f0: float, gain_db: float, q: float, fs: int = FS): """Audio EQ Cookbook peaking EQ (b, a), normalised by a0.""" 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 b / a[0], a / a[0] def rbj_highpass(f0: float, q: float, fs: int = FS): """Audio EQ Cookbook high-pass (b, a), normalised by a0.""" w0 = 2 * np.pi * f0 / fs alpha = np.sin(w0) / (2 * q) c = np.cos(w0) b = np.array([(1 + c) / 2, -(1 + c), (1 + c) / 2]) a = np.array([1 + alpha, -2 * c, 1 - alpha]) return b / a[0], a / a[0] def linear_phase_twin(b, a, n_taps: int = 16383) -> np.ndarray: """FIR with the biquad's magnitude and zero phase, delayed by (n_taps-1)/2 samples (frequency sampling).""" nfft = 1 << 16 w = np.linspace(0, np.pi, nfft // 2 + 1) _, H = signal.freqz(b, a, worN=w) mag = np.abs(H) h = np.fft.irfft(mag, nfft) # zero-phase, circular h = np.roll(h, (n_taps - 1) // 2)[:n_taps] h *= signal.windows.blackmanharris(n_taps) return h def process_side(L, R, side_filter, conv="half", m_delay: int = 0): """Encode, filter S only, delay M by m_delay samples, decode. side_filter maps an array to an array.""" M, S = encode(L, R, conv) S2 = side_filter(S) n = len(S2) M2 = np.zeros(n) M2[m_delay:m_delay + len(M)] = M[: max(0, n - m_delay)] return decode(M2, S2, conv) def transfer(y: np.ndarray, delay: int, nfft: int = 1 << 17): """Complex transfer function of an impulse response, with a known delay removed.""" Y = np.fft.rfft(y, nfft) f = np.fft.rfftfreq(nfft, 1 / FS) return f, Y * np.exp(2j * np.pi * f * delay / FS) # ---------------------------------------------------------------- test material def pinkish(n: int, rng: np.random.Generator) -> np.ndarray: """Gaussian noise through a -3 dB/octave approximation (Kellet's 3-pole filter).""" b = np.array([0.049922035, -0.095993537, 0.050612699, -0.004408786]) a = np.array([1.0, -2.494956002, 2.017265875, -0.522189400]) return signal.lfilter(b, a, rng.standard_normal(n)) def test_mix(rng: np.random.Generator, seconds: float = 3.0): """Centre element (identical in L and R), a hard-left and a hard-right element, independent noises.""" n = int(seconds * FS) centre = pinkish(n, rng) left = pinkish(n, rng) right = pinkish(n, rng) centre /= np.sqrt(np.mean(centre ** 2)) left /= np.sqrt(np.mean(left ** 2)) right /= np.sqrt(np.mean(right ** 2)) g_side = 10 ** (-6 / 20) # side elements 6 dB below the centre element L = centre + g_side * left R = centre + g_side * right fade = int(0.01 * FS) ramp = np.ones(n) ramp[:fade] = np.linspace(0, 1, fade) ramp[-fade:] = np.linspace(1, 0, fade) L, R = L * ramp, R * ramp peak = max(np.max(np.abs(L)), np.max(np.abs(R))) k = 10 ** (-6 / 20) / peak # -6 dBFS peak, room for +6 dB on S return L * k, R * k def stereo_level_db(L, R) -> float: """Mean power of the two channels, dB re full-scale sine-less RMS (10 log10 of mean square).""" return 10 * np.log10((np.mean(L ** 2) + np.mean(R ** 2)) / 2) def mono_level_db(L, R) -> float: return 10 * np.log10(np.mean(((L + R) / 2) ** 2)) def correlation(L, R) -> float: return float(np.sum(L * R) / np.sqrt(np.sum(L ** 2) * np.sum(R ** 2))) # ---------------------------------------------------------------- main def main() -> None: dc.clear_outputs(HERE, OUT) dc.setup_style() rng = np.random.default_rng(SEED) results: dict = {} L, R = test_mix(rng) # C1: lossless round trip c1 = [] for conv in ("half", "sqrt2"): M, S = encode(L, R, conv) L2, R2 = decode(M, S, conv) res64 = np.concatenate([L2 - L, R2 - R]) L32, R32 = L.astype(np.float32), R.astype(np.float32) M32, S32 = encode(L32, R32, conv) L3, R3 = decode(M32.astype(np.float32), S32.astype(np.float32), conv) res32 = np.concatenate([L3.astype(np.float64) - L32, R3.astype(np.float64) - R32]) L24, R24 = q24(L), q24(R) M24, S24 = encode(L24, R24, conv) L4, R4 = decode(q24(M24), q24(S24), conv) # M and S stored as 24-bit files in between res24 = np.concatenate([L4 - L24, R4 - R24]) c1.append({"convention": conv, "float64_residual_peak_dbfs": dbfs(res64), "float64_residual_rms_dbfs": rms_dbfs(res64), "float32_residual_peak_dbfs": dbfs(res32), "float32_residual_rms_dbfs": rms_dbfs(res32), "pcm24_between_residual_peak_dbfs": dbfs(res24), "pcm24_between_residual_rms_dbfs": rms_dbfs(res24)}) by_conv = {r["convention"]: r for r in c1} dc.write_curve_csv(OUT / "c1_round_trip.csv", "curve_csv", ["ビット数 [bit]", "2で割る形 [dBFS]", "√2で割る形 [dBFS]"], [np.array([24.0, 32.0, 64.0])] + [np.array([by_conv[c]["pcm24_between_residual_peak_dbfs"], by_conv[c]["float32_residual_peak_dbfs"], by_conv[c]["float64_residual_peak_dbfs"]]) for c in ("half", "sqrt2")], [0, 3, 3]) results["C1_round_trip"] = {"signal_peak_dbfs": dbfs(np.concatenate([L, R])), "signal_rms_dbfs": rms_dbfs(np.concatenate([L, R])), "rows": c1} # C2: level after mixed conventions; also how a centre-only signal reads on an M meter c2 = [] for enc in ("half", "sqrt2"): for dec in ("half", "sqrt2"): M, S = encode(L, R, enc) L2, R2 = decode(M, S, dec) c2.append({"encode": enc, "decode": dec, "level_change_db": stereo_level_db(L2, R2) - stereo_level_db(L, R)}) centre = np.sin(2 * np.pi * 1000 * np.arange(FS) / FS) * 0.5 meter = {} for conv in ("half", "sqrt2"): M, S = encode(centre, centre, conv) meter[conv] = {"M_peak_minus_L_peak_db": dbfs(M) - dbfs(centre), "S_peak_dbfs": dbfs(S)} results["C2_conventions"] = {"rows": c2, "centre_sine_on_M": meter} # C3: S gain vs stereo level, mono level, correlation gains = np.arange(-12.0, 6.01, 1.0) st, mo, corr, mono_dev = [], [], [], [] M, S = encode(L, R, "half") mono_ref = (L + R) / 2 for gdb in gains: L2, R2 = decode(M, S * 10 ** (gdb / 20), "half") st.append(stereo_level_db(L2, R2)) mo.append(mono_level_db(L2, R2)) corr.append(correlation(L2, R2)) mono_dev.append(dbfs((L2 + R2) / 2 - mono_ref)) st, mo, corr = np.array(st), np.array(mo), np.array(corr) i0 = int(np.where(gains == 0)[0][0]) i6 = int(np.where(gains == 6)[0][0]) im6 = int(np.where(gains == -6)[0][0]) # a hard-left element alone, S +6 dB Lh = np.ones(4) Rh = np.zeros(4) Mh, Sh = encode(Lh, Rh, "half") Lh2, Rh2 = decode(Mh, Sh * 2.0, "half") results["C3_side_gain"] = { "mix": "centre pinkish noise in both channels; hard-left and hard-right pinkish noises 6 dB lower; -6 dBFS peak", "rows": [{"side_gain_db": float(g), "stereo_level_db": float(a), "mono_level_db": float(b_), "stereo_minus_mono_db": float(a - b_), "correlation": float(c), "mono_change_peak_dbfs": float(d)} for g, a, b_, c, d in zip(gains, st, mo, corr, mono_dev)], "stereo_minus_mono_at_0db": float(st[i0] - mo[i0]), "stereo_minus_mono_at_plus6db": float(st[i6] - mo[i6]), "stereo_minus_mono_at_minus6db": float(st[im6] - mo[im6]), "max_mono_change_peak_dbfs": float(np.max(mono_dev)), "hard_left_element_side_plus6db": {"L_gain_db": float(20 * np.log10(abs(Lh2[0]))), "R_gain_db": float(20 * np.log10(abs(Rh2[0]))), "R_polarity": "inverted" if Rh2[0] < 0 else "same", "mono_gain_db": float(20 * np.log10(abs((Lh2[0] + Rh2[0]) / 2) / 0.5))}, } fig, axes = plt.subplots(1, 2, figsize=(8.4, 3.4)) ax = axes[0] ax.plot(gains, st - st[i0], color=dc.SERIES[0], label="stereo (mean of L and R)") ax.plot(gains, mo - mo[i0], color=dc.SERIES[1], label="mono fold-down (L+R)/2") ax.set_xlabel("gain applied to S (dB)") ax.set_ylabel("level change (dB)") ax.set_title("C3 Raising S: stereo moves, mono does not") ax.legend(loc="upper left") ax = axes[1] ax.plot(gains, corr, color=dc.SERIES[0]) ax.axhline(0, color=dc.TEXT_2, lw=0.8, ls="--") ax.set_ylim(-1, 1) ax.set_xlabel("gain applied to S (dB)") ax.set_ylabel("L/R correlation") ax.set_title("Correlation of the same mix") fig.tight_layout() dc.save_figure(fig, OUT / "c3_side_gain.png") dc.write_curve_csv(OUT / "c3_side_gain_levels.csv", "curve_csv", ["Sのゲイン [dB]", "ステレオの音量の変化 [dB]", "モノラルの音量の変化 [dB]"], [gains, st - st[i0], mo - mo[i0]], [0, 4, 4]) dc.write_curve_csv(OUT / "c3_side_gain_correlation.csv", "curve_csv", ["Sのゲイン [dB]", "左右の相関係数 [-]"], [gains, corr], [0, 4]) # C4-C6: filters on S, impulse in L only (a hard-left source) n_imp = 1 << 15 imp = np.zeros(n_imp) imp[0] = 1.0 zero = np.zeros(n_imp) b_pk, a_pk = rbj_peaking(3000.0, 6.0, 1.0) b_hp, a_hp = rbj_highpass(100.0, 1 / SQRT2) taps = 16383 D = (taps - 1) // 2 fir_pk = linear_phase_twin(b_pk, a_pk, taps) fir_hp = linear_phase_twin(b_hp, a_hp, taps) def iir(b, a): return lambda s: signal.lfilter(b, a, s) def fir(h): return lambda s: signal.fftconvolve(s, h) cases = { "peak_iir": (iir(b_pk, a_pk), 0, 0), "peak_lin": (fir(fir_pk), D, D), "side_late_32": (lambda s_: np.concatenate([np.zeros(32), s_])[: len(s_)], 0, 0), "hpf_iir": (iir(b_hp, a_hp), 0, 0), "hpf_lin": (fir(fir_hp), D, D), } tf = {} for name, (flt, mdel, ref_delay) in cases.items(): Lo, Ro = process_side(np.pad(imp, (0, taps)), np.pad(zero, (0, taps)), flt, "half", m_delay=mdel) f, HL = transfer(Lo, ref_delay) _, HR = transfer(Ro, ref_delay) tf[name] = (f, HL, HR) # centre source: L = R = impulse -> output must equal the (delayed) input Lc, Rc = process_side(np.pad(imp, (0, taps)), np.pad(imp, (0, taps)), flt, "half", m_delay=mdel) ref = np.zeros(len(Lc)) ref[mdel] = 1.0 tf[name + "_centre_dev"] = max(dbfs(Lc - ref), dbfs(Rc - ref)) # mono fold-down of the hard-left source must equal the (delayed) input / 2 whatever S went through mono_ref = np.zeros(len(Lo)) mono_ref[mdel] = 0.5 tf[name + "_mono_dev"] = dbfs((Lo + Ro) / 2 - mono_ref) def at(name, which, hz): f, HL, HR = tf[name] i = int(np.argmin(np.abs(f - hz))) H = HL if which == "L" else HR return float(20 * np.log10(max(abs(H[i]), 1e-300))) def max_in(name, which, lo, hi): f, HL, HR = tf[name] H = HL if which == "L" else HR sel = (f >= lo) & (f <= hi) v = 20 * np.log10(np.maximum(np.abs(H[sel]), 1e-300)) j = int(np.argmax(v)) return {"db": float(v[j]), "at_hz": float(f[sel][j])} def min_in(name, which, lo, hi): f, HL, HR = tf[name] H = HL if which == "L" else HR sel = (f >= lo) & (f <= hi) v = 20 * np.log10(np.maximum(np.abs(H[sel]), 1e-300)) j = int(np.argmin(v)) return {"db": float(v[j]), "at_hz": float(f[sel][j])} def l_out_max_abs_diff(a, b, lo, hi): f, HLa, _ = tf[a] _, HLb, _ = tf[b] sel = (f >= lo) & (f <= hi) d = np.abs(20 * np.log10(np.abs(HLa[sel])) - 20 * np.log10(np.abs(HLb[sel]))) j = int(np.argmax(d)) return {"db": float(d[j]), "at_hz": float(f[sel][j])} probe = [500, 1000, 1500, 2000, 3000, 4500, 6000, 9000, 15000] results["C4_side_peaking_eq"] = { "filter": "peaking +6 dB at 3 kHz, Q 1 (Audio EQ Cookbook biquad) on S; linear-phase twin: FIR 16383 taps, same magnitude", "source": "impulse in L only (hard-left element)", "centre_source_max_deviation_dbfs": {k: tf[k + "_centre_dev"] for k in ("peak_iir", "peak_lin")}, "mono_fold_down_deviation_peak_dbfs": {k: tf[k + "_mono_dev"] for k in ("peak_iir", "peak_lin")}, "L_out_db": {k: {str(h): at(k, "L", h) for h in probe} for k in ("peak_iir", "peak_lin")}, "R_out_db": {k: {str(h): at(k, "R", h) for h in probe} for k in ("peak_iir", "peak_lin")}, "R_out_max_20_20k": {k: max_in(k, "R", 20, 20000) for k in ("peak_iir", "peak_lin")}, "L_out_max_abs_diff_iir_vs_lin_100_20k": l_out_max_abs_diff("peak_iir", "peak_lin", 100, 20000), } results["C5_side_highpass"] = { "filter": "high-pass 100 Hz, Q 1/sqrt2 (Audio EQ Cookbook biquad) on S; linear-phase twin as above", "source": "impulse in L only (hard-left element)", "centre_source_max_deviation_dbfs": {k: tf[k + "_centre_dev"] for k in ("hpf_iir", "hpf_lin")}, "mono_fold_down_deviation_peak_dbfs": {k: tf[k + "_mono_dev"] for k in ("hpf_iir", "hpf_lin")}, "L_out_db": {k: {str(h): at(k, "L", h) for h in (25, 50, 100, 200, 400, 800, 1600)} for k in ("hpf_iir", "hpf_lin")}, "R_out_db": {k: {str(h): at(k, "R", h) for h in (25, 50, 100, 200, 400, 800, 1600)} for k in ("hpf_iir", "hpf_lin")}, } f6, HL6, HR6 = tf["side_late_32"] sel6 = (f6 >= 20) & (f6 <= 20000) def notches(H, count=4): v = 20 * np.log10(np.maximum(np.abs(H[sel6]), 1e-15)) # an exact zero reads as -300 dB idx = signal.argrelmin(v, order=50)[0] idx = idx[v[idx] < -40][:count] return [{"at_hz": float(f6[sel6][i]), "db": float(v[i])} for i in idx] results["C6_side_late"] = { "case": "S delayed by 32 samples (0.667 ms at 48 kHz), no EQ, M not delayed; impulse in L only", "delay_samples": 32, "L_out_notches": notches(HL6), "R_out_notches": notches(HR6), "L_out_max_db": float(np.max(20 * np.log10(np.maximum(np.abs(HL6[sel6]), 1e-15)))), "centre_source_max_deviation_dbfs": tf["side_late_32_centre_dev"], "mono_fold_down_deviation_peak_dbfs": tf["side_late_32_mono_dev"], "expected_L_notches_hz": "fs (2k+1) / (2*32) = 750, 2250, 3750 ... Hz", "expected_R_notches_hz": "fs k / 32 = 0, 1500, 3000 ... Hz", "L_out_db_at": {str(h): at("side_late_32", "L", h) for h in (750, 1500, 2250, 3000, 3750)}, "R_out_db_at": {str(h): at("side_late_32", "R", h) for h in (750, 1500, 2250, 3000, 3750)}, } def plot_pair(ax, names, labels, which, title, xlim, ylim): for i, (nm, lb) in enumerate(zip(names, labels)): f, HL, HR = tf[nm] H = HL if which == "L" else HR sel = (f >= xlim[0]) & (f <= xlim[1]) ax.plot(f[sel], 20 * np.log10(np.maximum(np.abs(H[sel]), 1e-12)), color=dc.SERIES[i], label=lb, lw=1.8) ax.set_xscale("log") ax.set_xlim(*xlim) dc.hz_axis(ax) ax.set_ylim(*ylim) ax.set_title(title) ax.set_ylabel("dB re input") fig, axes = plt.subplots(2, 1, figsize=(7.4, 5.6), sharex=True) plot_pair(axes[0], ["peak_iir", "peak_lin"], ["IIR biquad EQ on S", "linear-phase EQ on S (M delayed to match)"], "L", "C4 Hard-left source, +6 dB at 3 kHz on S: left output", (100, 20000), (-3, 5)) plot_pair(axes[1], ["peak_iir", "peak_lin"], ["IIR biquad EQ on S", "linear-phase EQ on S (M delayed to match)"], "R", "Right output (what leaks to the other side)", (100, 20000), (-60, 0)) axes[1].set_xlabel("frequency (Hz)") axes[0].legend(loc="upper left") fig.tight_layout() dc.save_figure(fig, OUT / "c4_side_eq_hard_left.png") fig, axes = plt.subplots(2, 1, figsize=(7.4, 5.6), sharex=True) plot_pair(axes[0], ["hpf_iir", "hpf_lin"], ["IIR biquad high-pass on S", "linear-phase high-pass on S (M delayed)"], "L", "C5 Hard-left source, S high-passed at 100 Hz: left output", (20, 5000), (-9, 1)) plot_pair(axes[1], ["hpf_iir", "hpf_lin"], ["IIR biquad high-pass on S", "linear-phase high-pass on S (M delayed)"], "R", "Right output", (20, 5000), (-60, 0)) axes[1].set_xlabel("frequency (Hz)") axes[0].legend(loc="lower right") fig.tight_layout() dc.save_figure(fig, OUT / "c5_side_highpass.png") fig, ax = plt.subplots(figsize=(7.4, 3.4)) for i, (H, lab) in enumerate(((HL6, "left output"), (HR6, "right output"))): sel = (f6 >= 100) & (f6 <= 8000) ax.plot(f6[sel], 20 * np.log10(np.maximum(np.abs(H[sel]), 1e-12)), color=dc.SERIES[i], label=lab, lw=1.5) ax.set_xscale("log") ax.set_xlim(100, 8000) dc.hz_axis(ax) ax.set_ylim(-40, 3) ax.set_xlabel("frequency (Hz)") ax.set_ylabel("dB re input") ax.set_title("C6 Hard-left source, S path 32 samples late, no EQ: comb filtering") ax.legend(loc="lower left") fig.tight_layout() dc.save_figure(fig, OUT / "c6_side_late.png") def response_grid(lo, hi, per_octave, extra): 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)) def response_csv(name, cases_, labels, which, grid): cols = [grid] for nm, ch in zip(cases_, which if isinstance(which, list) else [which] * len(cases_)): f, HL, HR = tf[nm] H = HL if ch == "L" else HR mag = np.abs(H) v = np.interp(grid, f, 20 * np.log10(np.where(mag > 0, mag, np.nan))) cols.append(v) dc.write_curve_csv(OUT / name, "frequency_response_csv", ["周波数 [Hz]"] + [f"{lb} [dB]" for lb in labels], cols, [3] + [4] * len(labels)) g4 = response_grid(20, 20000, 48, probe) response_csv("c4_side_eq_left.csv", ["peak_iir", "peak_lin"], ["IIR", "リニアフェイズ"], "L", g4) response_csv("c4_side_eq_right.csv", ["peak_iir", "peak_lin"], ["IIR", "リニアフェイズ"], "R", g4) g5 = response_grid(20, 20000, 48, [25, 50, 100, 200, 400, 800, 1600]) response_csv("c5_side_highpass_left.csv", ["hpf_iir", "hpf_lin"], ["IIR", "リニアフェイズ"], "L", g5) response_csv("c5_side_highpass_right.csv", ["hpf_iir", "hpf_lin"], ["IIR", "リニアフェイズ"], "R", g5) g6 = response_grid(100, 8000, 96, []) response_csv("c6_side_late.csv", ["side_late_32", "side_late_32"], ["左の出力", "右の出力"], ["L", "R"], g6) # Sound examples (stereo, 24-bit): original, S +6 dB, and the mono fold-down of each M, S = encode(L, R, "half") Lp, Rp = decode(M, S * 2.0, "half") clips = { "s1_mix_original.wav": np.column_stack([L, R]), "s2_mix_side_plus6db.wav": np.column_stack([Lp, Rp]), "s3_mono_of_original.wav": np.column_stack([(L + R) / 2, (L + R) / 2]), "s4_mono_of_side_plus6db.wav": np.column_stack([(Lp + Rp) / 2, (Lp + Rp) / 2]), } for name, x in clips.items(): dc.write_wav(OUT / name, x) results["sound_examples"] = {"files": list(clips), "committed": False, "note": "same gain for all files; s3 and s4 are the same signal (see C3)", "s3_minus_s4_peak_dbfs": dbfs(clips["s3_mono_of_original.wav"] - clips["s4_mono_of_side_plus6db.wav"])} detail = { "evidence_id": EVIDENCE_ID, "article_slug": "mid-side-processing", "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": { "conventions": "half: M=(L+R)/2, S=(L-R)/2, L=M+S, R=M-S; sqrt2: M=(L+R)/sqrt2, S=(L-R)/sqrt2, L=(M+S)/sqrt2, R=(M-S)/sqrt2", "levels": "stereo level = 10 log10 of the mean of the two channels' mean squares; mono = (L+R)/2", "filters": "Audio EQ Cookbook biquads (peaking, high-pass) via scipy.signal.lfilter; linear-phase twins by frequency sampling of the biquad magnitude, 16383 taps, Blackman-Harris window", "transfer_functions": "impulse in L only through encode -> filter on S -> decode, rfft with 131072 points, known latency removed; CSV values interpolated in dB on a log grid", }, "results": results, } half = by_conv["half"] ledger = { "summary_editorial": "合成した信号で、M/Sの変換と逆変換、Sの増減、Sだけへのフィルターと遅延の影響を確かめた実証の記録です。", "subject_ja": "ステレオ信号をM/Sに変換して戻したときの残差と、Sの増減やSだけへのEQ・ハイパス・遅延が左右とモノラルに与える影響を確かめる実証です。", "material_description_ja": "スクリプトが乱数から生成した3秒のステレオの雑音（中央・左だけ・右だけの成分）と、インパルス、1 kHzのサイン波です。", "settings": [ {"name": "サンプリングレート", "value": "48000 Hz"}, {"name": "M/Sの係数（結果の残差）", "value": "M=(L+R)/2、S=(L-R)/2、64ビット浮動小数点"}, {"name": "比べた係数", "value": "2で割る形と√2で割る形"}, {"name": "Sのゲイン", "value": "-12 dBから+6 dBまで1 dB刻み"}, {"name": "SのEQ", "value": "3 kHz、+6 dB、Q 1の双2次フィルターと、同じ振幅特性のリニアフェイズFIR"}, {"name": "Sのハイパス", "value": "100 Hz、Q 0.707の双2次フィルターと、同じ振幅特性のリニアフェイズFIR"}, {"name": "リニアフェイズFIR", "value": "16383タップ、Mを8191サンプル遅らせて合わせる"}, {"name": "Sの経路の遅延", "value": "32サンプル（EQなし、Mは遅らせない）"}, ], "comparison_conditions_ja": ["同じ seed で生成した同じ信号を、変換の前後と、1つの条件だけ変えた処理の間で比べました。"], "procedure_ja": [ "乱数の seed を固定し、中央・左だけ・右だけの成分を持つステレオの雑音を生成しました。", "MとSに変換してすぐに戻し、元の信号との差のピークとRMSを求めました。", "変換と逆変換の係数の組み合わせを変え、左右の音量の変化を求めました。", "Sのゲインを変え、ステレオとモノラルの音量、左右の相関を求めました。", "左だけのインパルスを、Sだけにフィルターや遅延をかけた経路に通し、左右の出力の周波数特性を求めました。", ], "results": { "null_residual_peak_dbfs": {"value": round(half["float64_residual_peak_dbfs"], 2), "decimals": 2}, "null_residual_rms_dbfs": {"value": round(half["float64_residual_rms_dbfs"], 2), "decimals": 2}, }, "signal_seed": SEED, } csv_kinds = {"c1_round_trip.csv": "curve_csv", "c3_side_gain_levels.csv": "curve_csv", "c3_side_gain_correlation.csv": "curve_csv", "c4_side_eq_left.csv": "frequency_response_csv", "c4_side_eq_right.csv": "frequency_response_csv", "c5_side_highpass_left.csv": "frequency_response_csv", "c5_side_highpass_right.csv": "frequency_response_csv", "c6_side_late.csv": "frequency_response_csv"} figures.write_all(OUT) rec = dc.write_outputs(HERE, OUT, EVIDENCE_ID, csv_kinds, detail, ledger) print(rec) if __name__ == "__main__": main()