# SPDX-License-Identifier: MIT # © 2026 DTMFRAME """Synthetic-signal demonstration for the article amp-sim-hardware-software (アンプシミュレーターはハードかソフトか). Questions it answers with numbers (no product, no licence, no recorded material): D1 Latency arithmetic: the length of one audio buffer (16-1024 samples at 44.1, 48 and 96 kHz) in ms and as the distance sound travels in that time, a round-trip estimate (input buffer + output buffer + a swept converter-and-driver overhead), and a few round trips derived from the components two makers state (MOTU M2, RME Babyface Pro FS). D2 Loading of a passive pickup by the input, with the pickup model of Paiva, Pakarinen and Välimäki (2012): an EMF with series inductance L and resistance R, a shunt capacitance C with a magnetic-loss resistance Rl across it, then the cable's capacitance and the input impedance. Resonance, peak height re 100 Hz and the 100 Hz level, for 10 kOhm to 1 MOhm inputs. D3 Aliasing of a static distortion stage y = tanh(g x)/tanh(g) fed a sine at a guitar note, run at 48 kHz directly and with 2x-16x oversampling: alias power, the strongest alias, and the numerical floor (coherent sampling). Every parameter is in PARAMS with its source: a citation for a value from literature or a maker's document, None for a condition chosen for the demonstration. All of PARAMS goes into results_detail.json. It also writes the article's concept diagrams (figures.py): out/concept_signal_path.svg, and from the results computed here concept_latency_ruler.svg (D1), concept_pickup_circuit.svg (D2), concept_aliasing_fold.svg (D3) and concept_choice_flow.svg (the conclusion table). 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 = 20261016 OUT = HERE / "out" EVIDENCE_ID = "ev_msr_synth_amp_sim_1" SLUG = "amp-sim-hardware-software" def _p(value, unit: str, source: str | None, note: str) -> dict: return {"value": value, "unit": unit, "source": source, "note": note} SRC_SPEED_OF_SOUND = ('J. O. Smith III, Physical Audio Signal Processing, "Speed of Sound in Air", Eq. (B.33): ' '"At 20 degrees Celsius, we get 1124.1 feet per second" ' '(https://www.dsprelated.com/freebooks/pasp/Speed_Sound_Air.html); 1124.1 ft/s x 0.3048 = 342.6 m/s') SRC_MOTU_M2 = "MOTU M2 specifications, https://motu.com/en-us/products/m-series/m2/specs/ : round-trip latency 2.5 ms at 96 kHz with a 32-sample buffer" SRC_RME_BFPFS = ("RME Babyface Pro FS User's Guide v1.5, section 28.2, pp. 81-82: total = 2 x buffer + AD 5 + DA 7 samples " "+ safety buffer 32 samples (output) + safety offset 16 samples per direction (Core Audio, macOS); " "without the offsets under ASIO") SRC_PICKUP_MODEL = ('Paiva, Pakarinen and Välimäki, "Acoustics and Modeling of Pickups", ' 'JAES 60(10), 2012, Fig. 12(a) and Eq. (30), p. 775') SRC_CABLE = 'Paiva and Penttinen, "Cable matters", AES 131st Convention, 2011, paper 8466' SRC_OVERSAMPLING = ("Bilbao, Esqueda, Parker and Välimäki, IEEE Signal Processing Letters 24(7), 2017, " "DOI 10.1109/LSP.2017.2675541: 8x or 16x oversampling is typical in audio") SRC_ALIASING = "Parker, Zavalishin and Le Bivic, DAFx-16, 2016, pp. 137-144" PARAMS = { # ---- D1 latency arithmetic "d1_sample_rates_hz": _p([44100, 48000, 96000], "Hz", None, "sample rates offered in DAW settings; chosen"), "d1_buffer_sizes_samples": _p([16, 32, 64, 128, 256, 512, 1024], "samples", None, "buffer sizes offered in DAW and driver settings; chosen"), "d1_speed_of_sound": _p(342.6, "m/s", SRC_SPEED_OF_SOUND, "air at 20 degrees C"), "d1_overhead_sweep": _p([0.0, 1.0, 2.0, 3.0, 4.0], "ms", None, "converter filters, USB transfer and driver safety buffers together; a sweep, not a product value " "(the published examples below have 0.92-1.83 ms outside the two buffers)"), "d1_overhead_for_figure": _p(2.0, "ms", None, "the middle of the sweep, used for the round-trip figure CSV"), "d1_figure_case": _p({"buffer_samples": 64, "sample_rate_hz": 48000}, "-", None, "the case quoted in the concept diagram and the headline"), "d1_example_motu_m2": _p({"round_trip_ms": 2.5, "sample_rate_hz": 96000, "buffer_samples": 32}, "-", SRC_MOTU_M2, "a stated total; the part outside the two buffers is derived here"), "d1_example_rme_babyface_pro_fs": _p({"ad_samples": 5, "da_samples": 7, "output_safety_buffer_samples": 32, "core_audio_safety_offset_samples_per_direction": 16, "cases_sample_rate_hz": 48000, "cases_buffer_samples": [32, 64, 128]}, "samples", SRC_RME_BFPFS, "stated components; totals derived here"), # ---- D2 pickup loading "d2_pickup_baseline": _p({"inductance_h": 2.0, "resistance_ohm": 10e3, "capacitance_f": 50e-12, "loss_resistance_ohm": 1000e3}, "-", SRC_PICKUP_MODEL + "; values from the caption of Fig. 13", "series L and R, shunt C with the magnetic-loss resistance Rl in parallel"), "d2_pickup_L5H_R12k_C100pF": _p({"inductance_h": 5.0, "resistance_ohm": 12e3, "capacitance_f": 100e-12, "loss_resistance_ohm": 1000e3}, "-", None, "a second pickup chosen inside the stated ranges (d2_pickup_ranges); Rl as the baseline; " "labelled by its values only (no source for a pickup type)"), "d2_pickup_ranges": _p({"inductance_h": [1.0, 10.0], "resistance_ohm": [5e3, 15e3], "capacitance_pf": [15.0, 200.0]}, "-", SRC_PICKUP_MODEL, "ranges stated for passive pickups"), "d2_cable_capacitance_total": _p({"main": 400e-12, "bounds": [225e-12, 843e-12]}, "F", SRC_CABLE + ': "typical cable" 400 pF in total (p. 8); measured 38-141 pF/m, ' '225-843 pF in total (Table 2)', "total cable capacitance, main case and bounds"), "d2_measured_resonance_range": _p([2760.0, 4140.0], "Hz", SRC_CABLE, "resonances measured with real 2.5-9 m cables; a sanity check only, the model is not tuned to it"), "d2_input_impedances": _p([ {"ohm": 10e3, "source": "Line 6 Support FAQ: the Helix Aux In is 10 kOhm (a line-level input)"}, {"ohm": 20e3, "source": "Zoom AMS-22 specifications: LINE IN 20 kOhm"}, {"ohm": 330e3, "source": "Zoom AMS-24 specifications: GUITAR input 330 kOhm"}, {"ohm": 500e3, "source": 'Roland Rubix22/24 specifications: Hi-Z input "500 k ohms or greater"'}, {"ohm": 1e6, "source": "1 MOhm stated by Kemper, Fractal, BOSS (GT-1000, GX-100), Strymon (Iridium), " "RME (Babyface Pro FS) and Zoom (AMS-22 GUITAR)"}, ], "ohm", "makers' specifications (one per value)", "input impedance R_in"), "d2_probe_frequencies": _p([1000.0, 2000.0, 3000.0, 5000.0], "Hz", None, "levels re 100 Hz reported at these frequencies"), "d2_band": _p([20.0, 20000.0], "Hz", None, "frequency range searched and drawn"), "d2_peak_min_db": _p(0.5, "dB", None, "a maximum inside the band counts as a resonant peak when it is at least this far above the 100 Hz level"), # ---- D3 aliasing "d3_a4_reference": _p(440.0, "Hz", None, "stated convention: A4 = 440 Hz, twelve-tone equal temperament"), "d3_notes": _p({"E4": {"semitones_from_a4": -5, "on_guitar": "open high E string"}, "E5": {"semitones_from_a4": 7, "on_guitar": "high E string, 12th fret"}, "E6": {"semitones_from_a4": 19, "on_guitar": "high E string, 24th fret"}}, "-", None, "chosen"), "d3_gains": _p([3.0, 10.0, 30.0, 100.0], "-", None, "g in y = tanh(g x)/tanh(g); chosen"), "d3_input_peak": _p(1.0, "-", None, "peak of the input sine (full scale of the shaper's input)"), "d3_oversampling_factors": _p([1, 2, 4, 8, 16], "-", SRC_OVERSAMPLING, "1 = the shaper runs at 48 kHz directly"), "d3_resampling_filter": _p({"passband_edge_hz": 20000.0, "stopband_edge_hz": 24000.0, "stopband_attenuation_db": 150.0}, "-", None, "Kaiser-window FIR (scipy.signal.kaiserord, firwin) at the oversampled rate, the same for up and down"), "d3_fft_length": _p(1 << 20, "samples", None, "0.0458 Hz bins at 48 kHz; each note is rounded to the nearest bin (coherent sampling)"), "d3_window": _p("rectangular", "-", None, "coherent sampling: an integer number of cycles in the block, so no window is needed"), "d3_near_linear_gain": _p(0.01, "-", None, "the floor run: nearly linear shaper"), "d3_spectrum_case": _p({"gain": 30.0, "note": "E6", "factors": [1, 8], "bin_hz": 10.0, "display_floor_db": -200.0}, "-", None, "the spectrum CSV: the largest bin in each 10 Hz band, dB re the fundamental's bin; values below the display floor " "are written as the display floor"), "d3_figure_note": _p("E6", "-", None, "the note of the alias-versus-oversampling figure CSV"), "d3_pad_samples": _p(4096, "samples", None, "samples dropped at each end (the resampling filters' start-up; checked to be longer)"), } P = {k: v["value"] for k, v in PARAMS.items()} 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)) def ohm_label(r: float) -> str: return f"{r / 1e6:g} MΩ" if r >= 1e6 else f"{r / 1e3:g} kΩ" # ---------------------------------------------------------------- D1 EXAMPLE_LABEL = "derived from the maker's stated components" def d1_examples(c: float) -> list[dict]: """Round trips derived from what two makers state (totals in samples, ms and cm; the part outside the two buffers).""" def row(device, driver, fs, buf, total, source, how): return {"label": EXAMPLE_LABEL, "device": device, "driver": driver, "sample_rate_hz": fs, "buffer_samples": buf, "total_samples": total, "total_ms": total / fs * 1000.0, "total_cm": total / fs * c * 100.0, "two_buffers_samples": 2 * buf, "outside_buffers_samples": total - 2 * buf, "outside_buffers_ms": (total - 2 * buf) / fs * 1000.0, "how": how, "source": source} m = P["d1_example_motu_m2"] rows = [row("MOTU M2", "as stated", m["sample_rate_hz"], m["buffer_samples"], m["round_trip_ms"] * m["sample_rate_hz"] / 1000.0, PARAMS["d1_example_motu_m2"]["source"], "stated round trip x sample rate; outside the buffers = total - 2 x buffer")] r = P["d1_example_rme_babyface_pro_fs"] fixed = r["ad_samples"] + r["da_samples"] + r["output_safety_buffer_samples"] for buf in r["cases_buffer_samples"]: for driver, extra, how in (("Core Audio (macOS)", 2 * r["core_audio_safety_offset_samples_per_direction"], "2 x buffer + AD + DA + output safety buffer + 2 x safety offset"), ("ASIO (no safety offsets)", 0, "2 x buffer + AD + DA + output safety buffer")): rows.append(row("RME Babyface Pro FS", driver, r["cases_sample_rate_hz"], buf, 2 * buf + fixed + extra, PARAMS["d1_example_rme_babyface_pro_fs"]["source"], how)) return rows def d1_latency() -> tuple[dict, dict]: c = P["d1_speed_of_sound"] rates, sizes, sweep = P["d1_sample_rates_hz"], P["d1_buffer_sizes_samples"], P["d1_overhead_sweep"] rows = [] for fs in rates: for n in sizes: one_ms = n / fs * 1000.0 rt = {f"{o:g}": 2 * one_ms + o for o in sweep} rows.append({"sample_rate_hz": fs, "buffer_samples": n, "one_buffer_ms": one_ms, "one_buffer_cm": one_ms / 1000.0 * c * 100.0, "round_trip_ms_by_overhead_ms": rt, "round_trip_cm_by_overhead_ms": {k: v / 1000.0 * c * 100.0 for k, v in rt.items()}}) case = P["d1_figure_case"] head = next(r for r in rows if r["sample_rate_hz"] == case["sample_rate_hz"] and r["buffer_samples"] == case["buffer_samples"]) o_fig = f"{P['d1_overhead_for_figure']:g}" ex = d1_examples(c) res = { "formula": "one buffer (ms) = buffer samples / sample rate * 1000; distance (cm) = ms / 1000 * c * 100; " "round trip (ms) = 2 * one buffer + overhead", "speed_of_sound_m_per_s": c, "speed_of_sound_source": PARAMS["d1_speed_of_sound"]["source"], "cm_per_ms": c / 1000.0 * 100.0, "overhead_ms_sweep": sweep, "overhead_note": "the overhead is a sweep standing for converter filters, USB transfer and driver safety buffers; " "it is not a value of any product. The published examples below show what two makers' own figures give.", "rows": rows, "headline": {"buffer_samples": head["buffer_samples"], "sample_rate_hz": head["sample_rate_hz"], "one_buffer_ms": head["one_buffer_ms"], "one_buffer_cm": head["one_buffer_cm"], "two_buffers_ms": 2 * head["one_buffer_ms"], "two_buffers_cm": 2 * head["one_buffer_cm"], "round_trip_ms_with_figure_overhead": head["round_trip_ms_by_overhead_ms"][o_fig], "round_trip_cm_with_figure_overhead": head["round_trip_cm_by_overhead_ms"][o_fig], "figure_overhead_ms": P["d1_overhead_for_figure"]}, "published_examples": ex, "published_examples_outside_buffers_ms_range": [min(e["outside_buffers_ms"] for e in ex), max(e["outside_buffers_ms"] for e in ex)], } x = np.array(sizes, dtype=float) one = {fs: np.array([n / fs * 1000.0 for n in sizes]) for fs in rates} rtf = {fs: 2 * one[fs] + P["d1_overhead_for_figure"] for fs in rates} return res, {"x": x, "one": one, "rt": rtf} # ---------------------------------------------------------------- D2 PICKUPS = {"baseline": "d2_pickup_baseline", "L5H_R12k_C100pF": "d2_pickup_L5H_R12k_C100pF"} def pickup_h(f, pu: dict, c_cable: float, r_in: float) -> np.ndarray: """Paiva, Pakarinen and Välimäki (2012), Fig. 12(a): H = Z / (R + jwL + Z), Z = Rl || R_in || 1/(jw (C + C_cable)).""" s = 2j * np.pi * np.atleast_1d(np.asarray(f, dtype=float)) z = 1.0 / (1.0 / pu["loss_resistance_ohm"] + 1.0 / r_in + s * (pu["capacitance_f"] + c_cable)) return z / (pu["resistance_ohm"] + s * pu["inductance_h"] + z) def d2_case(pu: dict, c_cable: float, r_in: float) -> dict: lo, hi = P["d2_band"] H = lambda f: pickup_h(f, pu, c_cable, r_in) # noqa: E731 lev = lambda f: float(dc.db(H(f))[0]) # noqa: E731 g = np.geomspace(lo, hi, 20001) v = dc.db(H(g)) i = int(np.argmax(v)) l100 = lev(100.0) interior = 0 < i < len(g) - 1 if interior: r = optimize.minimize_scalar(lambda f: -lev(f), bounds=(g[i - 1], g[i + 1]), method="bounded", options={"xatol": 1e-6}) f_pk, l_pk = float(r.x), -float(r.fun) else: f_pk, l_pk = float(g[i]), float(v[i]) has_peak = interior and l_pk - l100 >= P["d2_peak_min_db"] def cross_above(f_from: float, level: float): j0 = int(np.searchsorted(g, f_from)) below = np.nonzero(v[j0:] <= level)[0] if len(below) == 0: return None j = j0 + int(below[0]) if j == 0: return float(g[0]) return float(optimize.brentq(lambda f: lev(f) - level, g[j - 1], g[j], xtol=1e-6)) start = max(f_pk, 100.0) if has_peak else 100.0 r_shunt = 1.0 / (1.0 / pu["loss_resistance_ohm"] + 1.0 / r_in) c_total = pu["capacitance_f"] + c_cable return { "r_in_ohm": r_in, "cable_capacitance_pf": c_cable * 1e12, "total_capacitance_pf": c_total * 1e12, "shunt_resistance_ohm": r_shunt, "undamped_lc_resonance_hz": 1.0 / (2 * np.pi * np.sqrt(pu["inductance_h"] * c_total)), "has_peak_in_band": has_peak, "resonance_hz": f_pk if has_peak else None, "peak_db_re_100hz": (l_pk - l100) if has_peak else None, "level_db_re_emf_at_100hz": l100, "max_in_band_hz": f_pk, "max_in_band_db_re_100hz": l_pk - l100, "minus3db_re_100hz_hz": cross_above(start, l100 - 3.0), "minus3db_re_peak_above_peak_hz": cross_above(f_pk, l_pk - 3.0) if has_peak else None, "level_db_re_100hz": {f"{f:g}": lev(f) - l100 for f in P["d2_probe_frequencies"]}, "dc_gain_db": float(20 * np.log10(r_shunt / (pu["resistance_ohm"] + r_shunt))), } def d2_pickups() -> tuple[dict, dict]: cab = P["d2_cable_capacitance_total"] cables = [cab["main"]] + list(cab["bounds"]) r_ins = [e["ohm"] for e in P["d2_input_impedances"]] cases = [] for name, key in PICKUPS.items(): for c_cable in cables: for r_in in r_ins: cases.append({"pickup": name, **d2_case(P[key], c_cable, r_in)}) def pick(pu, c_pf, r): return next(c for c in cases if c["pickup"] == pu and round(c["cable_capacitance_pf"]) == c_pf and c["r_in_ohm"] == r) main_pf = round(cab["main"] * 1e12) sanity = pick("baseline", main_pf, 1e6) lo_m, hi_m = P["d2_measured_resonance_range"] res = { "model": "Paiva, Pakarinen and Välimäki (2012), Fig. 12(a) and Eq. (30): H(f) = Z / (R + j 2 pi f L + Z), " "Z = Rl || R_in || 1/(j 2 pi f (C + C_cable)); series L and R, shunt C with the magnetic-loss resistance Rl across it", "model_source": SRC_PICKUP_MODEL, "not_modelled": "the guitar's own controls (volume and tone pots, tone capacitor) are not modelled: no citable value was found; " "the load is the pickup's own C and Rl, the cable's capacitance and R_in only", "pickups": {name: {**P[key], "source": PARAMS[key]["source"]} for name, key in PICKUPS.items()}, "cable_capacitance_pf": {"main": main_pf, "bounds": [round(b * 1e12) for b in cab["bounds"]], "source": PARAMS["d2_cable_capacitance_total"]["source"]}, "input_impedances": P["d2_input_impedances"], "level_reference": "dB re the pickup's open-circuit EMF", "resonance": "the maximum of |H| on a 20001-point log grid from 20 Hz to 20 kHz, refined with scipy.optimize.minimize_scalar; " "a resonant peak (has_peak_in_band) is a maximum inside the band at least d2_peak_min_db above the 100 Hz level; " "otherwise resonance_hz is null and max_in_band_* still give the maximum (the low edge for a falling response)", "minus3db_points": "minus3db_re_100hz_hz: the first frequency above the resonant peak (above 100 Hz when there is none) " "where the level is 3 dB below the 100 Hz level; minus3db_re_peak_above_peak_hz: 3 dB below the peak level", "cases": cases, f"headline_baseline_{main_pf}pf": {ohm_label(r): {k: pick("baseline", main_pf, r)[k] for k in ("resonance_hz", "peak_db_re_100hz", "level_db_re_emf_at_100hz", "minus3db_re_100hz_hz", "minus3db_re_peak_above_peak_hz", "level_db_re_100hz")} for r in r_ins}, "sanity_check": {"case": f"baseline pickup, {main_pf} pF cable, R_in 1 MOhm", "resonance_hz": sanity["resonance_hz"], "measured_range_hz": [lo_m, hi_m], "measured_range_source": PARAMS["d2_measured_resonance_range"]["source"] + " (real 2.5-9 m cables)", "inside_measured_range": sanity["resonance_hz"] is not None and lo_m <= sanity["resonance_hz"] <= hi_m, "note": "reported as found; the model is not tuned to the measured range"}, } f = log_grid(P["d2_band"][0], P["d2_band"][1], 48) curves = {name: {r: dc.db(pickup_h(f, P[key], cab["main"], r)) for r in r_ins} for name, key in PICKUPS.items()} return res, {"f": f, "curves": curves, "r_ins": r_ins} # ---------------------------------------------------------------- D3 N_FFT = P["d3_fft_length"] DF = FS / N_FFT # 0.0457763671875 Hz PAD = P["d3_pad_samples"] def note_et_hz(name: str) -> float: """Equal-tempered frequency of a note from A4.""" return P["d3_a4_reference"] * 2.0 ** (P["d3_notes"][name]["semitones_from_a4"] / 12.0) def note_bin(name: str) -> int: """The note rounded to the nearest bin of the N-point FFT (coherent sampling).""" return int(round(note_et_hz(name) / DF)) def note_hz(name: str) -> float: return note_bin(name) * DF def shaper(x: np.ndarray, g: float) -> np.ndarray: return np.tanh(g * x) / np.tanh(g) def os_filter(factor: int) -> np.ndarray: spec = P["d3_resampling_filter"] rate = factor * FS width = (spec["stopband_edge_hz"] - spec["passband_edge_hz"]) / (rate / 2) n, beta = signal.kaiserord(spec["stopband_attenuation_db"], width) n |= 1 return signal.firwin(n, (spec["passband_edge_hz"] + spec["stopband_edge_hz"]) / 2, window=("kaiser", beta), fs=rate) def filter_report(h: np.ndarray, factor: int) -> dict: spec = P["d3_resampling_filter"] rate = factor * FS w, H = signal.freqz(h, worN=1 << 16, fs=rate) pb, sb = w <= spec["passband_edge_hz"], w >= spec["stopband_edge_hz"] return {"taps": len(h), "start_up_samples_at_48khz": int(np.ceil((len(h) - 1) / 2 / factor)), "passband_ripple_db": float(np.max(np.abs(dc.db(H[pb])))), "stopband_max_db": float(np.max(dc.db(H[sb])))} def render(g: float, k0: int, factor: int, filters: dict) -> np.ndarray: """N steady-state output samples: the input sine has exactly k0 cycles in N samples (phase index (k0 n) mod N, exact in integers), so after the PAD samples that hold the filters' start-up are dropped, the block is exactly periodic.""" n = np.arange(N_FFT + 2 * PAD, dtype=np.int64) x = P["d3_input_peak"] * np.sin(2 * np.pi * ((k0 * n) % N_FFT) / N_FFT) if factor == 1: y = shaper(x, g) else: h = filters[factor] y = signal.resample_poly(shaper(signal.resample_poly(x, factor, 1, window=h), g), 1, factor, window=h) return y[PAD:PAD + N_FFT] def power_spectrum(y: np.ndarray) -> np.ndarray: """One-sided power, rectangular window (coherent sampling: every component sits on one bin).""" p = np.abs(np.fft.rfft(y)) ** 2 p[1:-1] *= 2.0 return p def harmonic_bins(k0: int) -> np.ndarray: """Bins of the harmonics k f0 below fs/2 (k = 1, 2, ...).""" k = np.arange(1, (N_FFT // 2) // k0 + 1) b = k * k0 return b[b < N_FFT // 2] def amplitude_dbfs(bin_power: float) -> float: """Peak amplitude of a sinusoid on one bin (rectangular window: |X| = A N / 2), dB re 1.""" return float(10 * np.log10(2.0 * bin_power) - 20 * np.log10(N_FFT)) def db10(a: float, b: float) -> float: return float("-inf") if a == 0 else float(10 * np.log10(a / b)) def measure(p: np.ndarray, k0: int) -> dict: mask = np.ones(len(p), dtype=bool) mask[harmonic_bins(k0)] = False total, fund = float(np.sum(p)), float(p[k0]) pa = np.where(mask, p, 0.0) i = int(np.argmax(pa)) return {"alias_power_db_re_total": db10(float(np.sum(pa)), total), "strongest_alias_db_re_fundamental": db10(float(pa[i]), fund), "strongest_alias_hz": i * DF, "strongest_alias_dbfs": amplitude_dbfs(float(pa[i])) if pa[i] > 0 else float("-inf"), "fundamental_dbfs": amplitude_dbfs(fund)} def source_harmonic(alias_bin: int, k0: int, factor: int) -> dict | None: """The lowest harmonic above fs/2 whose frequency folds (mod fs) onto alias_bin, found in exact integer bin arithmetic. Above L fs/2 it folded at the shaper's rate and passed the decimation filter; between fs/2 and L fs/2 (L > 1) it leaked through the decimation filter's stopband and folded when the rate went back to fs.""" k = np.arange(1, 400001, dtype=np.int64) b = (k * k0) % N_FFT folded = np.minimum(b, N_FFT - b) hit = k[(folded == alias_bin) & (k * k0 > N_FFT // 2)] if len(hit) == 0: return None kk = int(hit[0]) fk = kk * k0 * DF return {"harmonic": kk, "harmonic_hz": fk, "folds_to_hz": alias_bin * DF, "shaper_rate_hz": factor * FS, "path": "folded at the shaper's rate" if fk > factor * FS / 2 else "leaked through the decimation filter's stopband"} def d3_aliasing() -> tuple[dict, dict]: notes = list(P["d3_notes"]) factors = P["d3_oversampling_factors"] filters = {L: os_filter(L) for L in factors if L > 1} assert all(int(np.ceil((len(h) - 1) / 2 / L)) < PAD for L, h in filters.items()), "PAD shorter than a filter's start-up" k0s = {nm: note_bin(nm) for nm in notes} floor_lin = {(nm, L): measure(power_spectrum(render(P["d3_near_linear_gain"], k0s[nm], L, filters)), k0s[nm]) for nm in notes for L in factors} rows, spectra = [], {} sc = P["d3_spectrum_case"] for g in P["d3_gains"]: for nm in notes: for L in factors: p = power_spectrum(render(g, k0s[nm], L, filters)) m = measure(p, k0s[nm]) fl_p = floor_lin[(nm, L)]["alias_power_db_re_total"] fl_s = floor_lin[(nm, L)]["strongest_alias_db_re_fundamental"] ok_s = m["strongest_alias_db_re_fundamental"] > fl_s + 3.0 rows.append({"gain": g, "note": nm, "note_hz": note_hz(nm), "oversampling": L, **m, "floor_alias_power_db": fl_p, "floor_strongest_alias_db": fl_s, "alias_power_resolved": m["alias_power_db_re_total"] > fl_p + 3.0, "strongest_alias_resolved": ok_s, "strongest_alias_source": source_harmonic(int(round(m["strongest_alias_hz"] / DF)), k0s[nm], L) if ok_s else None}) if g == sc["gain"] and nm == sc["note"] and L in sc["factors"]: spectra[L] = p def row(g, nm, L): return next(r for r in rows if r["gain"] == g and r["note"] == nm and r["oversampling"] == L) # spectrum case: the aliases by size, and the strongest one below the note itself k0 = k0s[sc["note"]] mask_s = np.ones(N_FFT // 2 + 1, dtype=bool) mask_s[harmonic_bins(k0)] = False peaks = {} for L, p in spectra.items(): rel = np.where(mask_s, p, 0.0) / float(p[k0]) idx = np.nonzero(rel >= 10 ** (-60 / 10))[0] top = idx[np.argsort(rel[idx])[::-1]][:8] j = int(np.argmax(rel[:k0])) fl_s = row(sc["gain"], sc["note"], L)["floor_strongest_alias_db"] lo_db = 10 * np.log10(rel[j]) if rel[j] > 0 else float("-inf") peaks[f"{L}x"] = {"strongest_below_fundamental": {"hz": float(j * DF), "db_re_fundamental": float(lo_db), "resolved": bool(lo_db > fl_s + 3.0), "source": source_harmonic(j, k0, L) if lo_db > fl_s + 3.0 else None}, "aliases_above_minus60db_re_fundamental": int(len(idx)), "top": [{"hz": float(i * DF), "db_re_fundamental": float(10 * np.log10(rel[i])), "source": source_harmonic(int(i), k0, L)} for i in sorted(top)]} floor_all = [floor_lin[k_]["alias_power_db_re_total"] for k_ in floor_lin] floor_one = [floor_lin[k_]["strongest_alias_db_re_fundamental"] for k_ in floor_lin] res = { "shaper": "y = tanh(g x) / tanh(g), x = sine of peak 1 (so a full-scale input reaches output 1)", "notes": {nm: {"equal_tempered_hz": note_et_hz(nm), "fft_bin": k0s[nm], "f0_hz": note_hz(nm), "deviation_hz": note_hz(nm) - note_et_hz(nm), "deviation_cents": float(1200 * np.log2(note_hz(nm) / note_et_hz(nm)))} for nm in notes}, "oversampling": "x upsampled by L with scipy.signal.resample_poly (the Kaiser FIR below), shaped at L * 48 kHz, " "then filtered and downsampled to 48 kHz with the same FIR; L = 1 shapes at 48 kHz with no filter", "resampling_filters": {f"{L}x": filter_report(filters[L], L) for L in filters}, "measurement": f"coherent sampling: each note rounded to the nearest bin of a {N_FFT}-point FFT ({DF} Hz bins), so the " f"input has an integer number of cycles in the block and every harmonic k f0 and every alias |m fs +- k f0| " f"falls exactly on a bin; the input's phase is computed from (k0 n) mod N in integers; {PAD} samples are dropped " f"at each end (the resampling filters' start-up is at most " f"{max(filter_report(filters[L], L)['start_up_samples_at_48khz'] for L in filters)} samples), leaving N exactly " f"periodic samples; rectangular window; harmonic bins = k k0 below fs/2, every other bin is alias; " f"alias power = their sum relative to the total output power; strongest alias = the largest such bin " f"relative to the fundamental's bin", "filter_limit_note": "where the strongest alias's path is 'leaked through the decimation filter's stopband', the value is set by the " "resampling filter (stopband about 146-149 dB, a design choice in PARAMS), not by the floor; this is why " "the low-gain curves flatten from 2x or 4x on", "classification_check": "an alias can only land on a harmonic bin when (k +- j) k0 is a multiple of N; with these k0 " "that needs harmonic numbers above 2^18, whose level is far below the float64 floor", "floor": { "method": f"the same measurement with the shaper nearly linear (g = {P['d3_near_linear_gain']:g}) at every note and factor; " "its harmonics above fs/2 are far below float64 resolution, so what remains is rounding; " "a value within 3 dB of its row's floor is marked not resolved", "near_linear_alias_power_db": {f"{nm}_{L}x": floor_lin[(nm, L)]["alias_power_db_re_total"] for nm in notes for L in factors}, "near_linear_strongest_alias_db": {f"{nm}_{L}x": floor_lin[(nm, L)]["strongest_alias_db_re_fundamental"] for nm in notes for L in factors}, "near_linear_fundamental_dbfs": {f"{nm}_{L}x": floor_lin[(nm, L)]["fundamental_dbfs"] for nm in notes for L in factors}, "alias_power_db_range": [min(floor_all), max(floor_all)], "strongest_alias_db_range": [min(floor_one), max(floor_one)], }, "rows": rows, "headline": {f"g{sc['gain']:g}_{sc['note']}": {f"{L}x": {k_: row(sc["gain"], sc["note"], L)[k_] for k_ in ("alias_power_db_re_total", "alias_power_resolved", "floor_alias_power_db", "strongest_alias_db_re_fundamental", "strongest_alias_hz", "strongest_alias_dbfs", "strongest_alias_resolved", "strongest_alias_source")} for L in factors}}, "spectrum_case_peaks": peaks, } return res, {"rows": rows, "spectra": spectra, "row": row} # ---------------------------------------------------------------- main def main() -> None: dc.clear_outputs(HERE, OUT) dc.setup_style() results: dict = {} # ------------------------------------------------------------ D1 r1, c1 = d1_latency() results["D1_latency_arithmetic"] = r1 case = P["d1_figure_case"] assert figures.FS_HZ == case["sample_rate_hz"] and figures.BUFFER_SAMPLES == case["buffer_samples"], "figures.py case differs from D1" assert figures.SPEED_OF_SOUND == P["d1_speed_of_sound"], "figures.py c differs from D1" assert abs(figures.buffer_ms() - r1["headline"]["one_buffer_ms"]) < 1e-12, "figures.py disagrees with D1" rates = P["d1_sample_rates_hz"] heads = [f"{fs / 1000:g} kHz [ms]" for fs in rates] dc.write_curve_csv(OUT / "d1_buffer_ms.csv", "curve_csv", ["バッファーサイズ [サンプル]"] + heads, [c1["x"]] + [c1["one"][fs] for fs in rates], [0] + [4] * len(rates)) dc.write_curve_csv(OUT / "d1_round_trip_ms.csv", "curve_csv", ["バッファーサイズ [サンプル]"] + heads, [c1["x"]] + [c1["rt"][fs] for fs in rates], [0] + [4] * len(rates)) fig, axes = plt.subplots(1, 2, figsize=(9.6, 3.8)) for ax, key, title in ((axes[0], "one", "one buffer"), (axes[1], "rt", f"round trip: 2 buffers + {P['d1_overhead_for_figure']:g} ms overhead")): for i, fs in enumerate(rates): ax.plot(c1["x"], c1[key][fs], color=dc.SERIES[i], marker="o", ms=4, label=f"{fs / 1000:g} kHz") ax.set_xscale("log", base=2) ax.set_yscale("log") ax.set_xticks(c1["x"], [f"{int(v)}" for v in c1["x"]]) lo_y = min(float(np.min(c1[key][fs])) for fs in rates) hi_y = max(float(np.max(c1[key][fs])) for fs in rates) yt = [t for t in (0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100) if lo_y / 1.5 <= t <= hi_y * 1.5] ax.set_yticks(yt, [f"{t:g}" for t in yt]) ax.minorticks_off() ax.set_xlabel("buffer size (samples)") ax.set_ylabel("ms") ax.set_title(f"D1 {title}") ax.legend(loc="upper left") fig.tight_layout() dc.save_figure(fig, OUT / "d1_latency.png") # ------------------------------------------------------------ D2 r2, c2 = d2_pickups() results["D2_pickup_loading"] = r2 r_ins = c2["r_ins"] main_pf = round(P["d2_cable_capacitance_total"]["main"] * 1e12) # The site's chart draws at most 4 series (apps/web article-figure FIGURE_SERIES_MAX), so the figure # leaves out 500 kohm, which lies between 330 kohm and 1 Mohm; the table and results keep all five. fig_r_ins = [r for r in r_ins if r != 500e3] dc.write_curve_csv(OUT / "d2_pickup_loading.csv", "frequency_response_csv", ["周波数 [Hz]"] + [f"入力 {ohm_label(r)} [dB]" for r in fig_r_ins], [c2["f"]] + [c2["curves"]["baseline"][r] for r in fig_r_ins], [3] + [3] * len(fig_r_ins)) r2["figure_csv"] = {"file": "d2_pickup_loading.csv", "case": f"baseline pickup, {main_pf} pF cable", "r_in_ohm": fig_r_ins, "values": "dB re the pickup's open-circuit EMF", "note": "500 kohm is left out of the figure (the site's chart draws at most 4 series); it is in the table"} fig, axes = plt.subplots(2, 1, figsize=(7.4, 6.4), sharex=True) for ax, key, title in ((axes[0], "baseline", "baseline pickup 2 H, 10 kohm, 50 pF"), (axes[1], "L5H_R12k_C100pF", "pickup 5 H, 12 kohm, 100 pF")): for i, r in enumerate(r_ins): ax.plot(c2["f"], c2["curves"][key][r], color=dc.SERIES[i], lw=1.6, label=ohm_label(r).replace("Ω", "ohm")) ax.set_xscale("log") ax.set_xlim(P["d2_band"][0], P["d2_band"][1]) dc.hz_axis(ax) ax.set_ylim(-40, 15) ax.set_ylabel("dB re pickup EMF") ax.set_title(f"D2 {title}, {main_pf} pF cable: input impedance") ax.legend(loc="lower left", ncol=3) axes[1].set_xlabel("frequency (Hz)") fig.tight_layout() dc.save_figure(fig, OUT / "d2_pickup_loading.png") # ------------------------------------------------------------ D3 r3, c3 = d3_aliasing() results["D3_aliasing"] = r3 rows = c3["rows"] dc.write_curve_csv(OUT / "d3_alias_power_table.csv", "curve_csv", ["行 [番号]", "ゲイン [倍]", "音の周波数 [Hz]", "オーバーサンプリング [倍]", "エイリアシングノイズの電力 [dB]", "最大のエイリアシングノイズ [dB]", "エイリアシングノイズの電力の下限 [dB]", "最大のエイリアシングノイズの下限 [dB]"], [np.arange(1, len(rows) + 1)] + [np.array([r[k] for r in rows], dtype=float) for k in ("gain", "note_hz", "oversampling", "alias_power_db_re_total", "strongest_alias_db_re_fundamental", "floor_alias_power_db", "floor_strongest_alias_db")], [0, 0, 2, 0, 2, 2, 2, 2]) nm_fig = P["d3_figure_note"] factors = np.array(P["d3_oversampling_factors"], dtype=float) fig_cols, clamped = [], [] for g in P["d3_gains"]: col = [] for L in P["d3_oversampling_factors"]: r = c3["row"](g, nm_fig, L) v = max(r["alias_power_db_re_total"], r["floor_alias_power_db"]) if not r["alias_power_resolved"]: clamped.append(f"g{g:g}_{L}x") col.append(v) fig_cols.append(np.array(col)) r3["figure_csv"] = {"file": "d3_alias_vs_oversampling.csv", "note": nm_fig, "values": "alias power in dB re total" + ("; a value within 3 dB of its floor (not resolved) is written as max(value, floor)" if clamped else "; every point is more than 3 dB above its floor"), "not_resolved_points": clamped} dc.write_curve_csv(OUT / "d3_alias_vs_oversampling.csv", "curve_csv", ["オーバーサンプリング [倍]"] + [f"ゲイン {g:g} [dB]" for g in P["d3_gains"]], [factors] + fig_cols, [0] + [2] * len(fig_cols)) sc = P["d3_spectrum_case"] f0s = note_hz(sc["note"]) c1b = note_bin(sc["note"]) edges = np.arange(0.0, FS / 2 + sc["bin_hz"] / 2, sc["bin_hz"]) fbin = np.arange(N_FFT // 2 + 1) * DF which = np.minimum(np.searchsorted(edges, fbin, side="right") - 1, len(edges) - 2) spec_cols = [] for L in sc["factors"]: p = c3["spectra"][L] ref = float(p[c1b]) hold = np.zeros(len(edges) - 1) np.maximum.at(hold, which, p / ref) spec_cols.append(np.maximum(dc.db(np.sqrt(hold)), sc["display_floor_db"])) centres = (edges[:-1] + edges[1:]) / 2 dc.write_curve_csv(OUT / "d3_spectrum_g30_e6.csv", "spectrum_csv", ["周波数 [Hz]"] + [f"{L}倍 [dB]" for L in sc["factors"]], [centres] + spec_cols, [0] + [2] * len(spec_cols)) r3["spectrum_csv"] = {"file": "d3_spectrum_g30_e6.csv", "case": sc, "values": f"the largest bin of the coherent (rectangular-window) power spectrum in each {sc['bin_hz']:g} Hz band " f"(x = band centre), dB re the fundamental's bin; values below {sc['display_floor_db']:g} dB are written as " f"{sc['display_floor_db']:g} dB (the numerical floor is lower: see floor)"} fig, axes = plt.subplots(2, 1, figsize=(7.4, 6.4)) ax = axes[0] for i, (g, col) in enumerate(zip(P["d3_gains"], fig_cols)): ax.plot(factors, col, color=dc.SERIES[i], marker="o", ms=5, label=f"g = {g:g}") fl = max(r["floor_alias_power_db"] for r in rows if r["note"] == nm_fig) y_lo = 20 * np.floor(min(float(np.min(c)) for c in fig_cols) / 20) - 10 ax.set_xscale("log", base=2) ax.set_xticks(factors, [f"{int(v)}x" for v in factors]) ax.minorticks_off() ax.set_ylim(y_lo, 0) ax.set_xlabel("oversampling") ax.set_ylabel("alias power (dB re total)") ax.set_title(f"D3 tanh shaper at {f0s:.2f} Hz (E6): alias power (numerical floor {fl:.0f} dB)") ax.legend(loc="upper right") ax = axes[1] for i, (L, col) in enumerate(zip(sc["factors"], spec_cols)): ax.plot(centres, col, color=dc.SERIES[i], lw=0.8, label=f"{L}x") ax.set_xlim(0, FS / 2) ax.set_ylim(sc["display_floor_db"] - 5, 5) ax.set_xlabel("frequency (Hz)") ax.set_ylabel("dB re fundamental") ax.set_title(f"g = {sc['gain']:g}, {f0s:.2f} Hz: largest bin in each 10 Hz band") ax.legend(loc="upper right") fig.tight_layout() dc.save_figure(fig, OUT / "d3_aliasing.png") # ------------------------------------------------------------ record detail = { "evidence_id": EVIDENCE_ID, "article_slug": SLUG, "signal_seed": SEED, "random_numbers_used": False, "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"), "params_note": "source: the citation for a value taken from literature or a maker's document; " "null marks a condition chosen for the demonstration", "literature_d3": [SRC_ALIASING, SRC_OVERSAMPLING], "params": PARAMS, "method": { "arithmetic": "float64 throughout", "d1": r1["formula"], "d2": r2["model"], "d3": r3["measurement"], }, "results": results, } pk = lambda d: (f"L {d['inductance_h']:g} H、R {d['resistance_ohm'] / 1e3:g} kΩ、C {d['capacitance_f'] * 1e12:g} pF、" # noqa: E731 f"Rl {d['loss_resistance_ohm'] / 1e3:g} kΩ") cab = P["d2_cable_capacitance_total"] ledger = { "summary_editorial": "合成した信号と式で、アンプシミュレーターのレイテンシー、ギター用の入力インピーダンス、歪みで出るエイリアシングノイズを確かめた実証の記録です。", "subject_ja": "バッファーサイズとレイテンシー、パッシブのピックアップと入力インピーダンス、歪みで出るエイリアシングノイズとオーバーサンプリングの関係を計算で確かめる実証です。", "material_description_ja": "スクリプトが式から生成したサイン波（ギターの音の高さ）と、ピックアップと配線の回路モデルです。乱数は使っていません。", "settings": [ {"name": "レイテンシーの計算のサンプリングレート", "value": "、".join(f"{fs / 1000:g} kHz" for fs in rates)}, {"name": "バッファーサイズ", "value": "、".join(str(n) for n in P["d1_buffer_sizes_samples"]) + "サンプル"}, {"name": "音の速さ", "value": f"{P['d1_speed_of_sound']:g} m/s（レイテンシーを距離に直すため）"}, {"name": "往復のレイテンシーの見積もり", "value": "入力と出力のバッファーに、変換とドライバーのレイテンシー" f"{min(P['d1_overhead_sweep']):g}〜{max(P['d1_overhead_sweep']):g} msを足す（製品の値ではない）"}, {"name": "ピックアップの回路モデル", "value": "直列のLとR、並列のCと磁気損失の抵抗Rl（Paivaら2012の図12(a)と式30）"}, {"name": "基準のピックアップ", "value": pk(P["d2_pickup_baseline"])}, {"name": "もう1つのピックアップ", "value": pk(P["d2_pickup_L5H_R12k_C100pF"])}, {"name": "ケーブルの静電容量", "value": f"合計{cab['main'] * 1e12:g} pF（" + "と".join(f"{b * 1e12:g} pF" for b in cab["bounds"]) + "も計算）"}, {"name": "ギターのつまみ", "value": "ボリュームとトーンは回路モデルに入れていない"}, {"name": "入力インピーダンス", "value": "、".join(ohm_label(e["ohm"]) for e in P["d2_input_impedances"])}, {"name": "歪みのモデル", "value": "tanh(g·x)/tanh(g)、g = " + "、".join(f"{g:g}" for g in P["d3_gains"])}, {"name": "歪みに入れた音", "value": "、".join(f"{nm} {note_hz(nm):.2f} Hz" for nm in P["d3_notes"]) + "（1弦の開放、12、24フレット）"}, {"name": "歪みの計算のサンプリングレート", "value": f"{FS} Hz、64ビット浮動小数点"}, {"name": "オーバーサンプリング", "value": "、".join(f"{L}倍" for L in P["d3_oversampling_factors"]) + "（カイザー窓のFIRでアップサンプリングとダウンサンプリング）"}, {"name": "エイリアシングノイズの測り方", "value": f"{N_FFT}点のFFT、矩形窓。音をFFTの刻みにそろえ、倍音以外のビンの電力をエイリアシングノイズとする"}, ], "comparison_conditions_ja": ["同じ式と同じ入力で、バッファーサイズ、入力インピーダンス、オーバーサンプリングの倍率だけを変えて比べました。"], "procedure_ja": [ "バッファーサイズとサンプリングレートからバッファー1つ分の時間を求め、音が空気中を進む距離に直しました。", "入力と出力のバッファーに変換とドライバーのレイテンシーを足して、往復のレイテンシーを見積もりました。", "ピックアップを起電力、直列のインダクタンスと抵抗、並列の静電容量と損失の抵抗の回路モデルにし、ケーブルと入力インピーダンスをつないだ周波数特性を求めました。", "共振周波数、100 Hzに対する共振ピークの高さ、起電力に対する100 Hzの大きさを、場合ごとに求めました。", "ギターの音の高さのサイン波（FFTの刻みにそろえた周波数）をtanhの歪みに通し、オーバーサンプリングの倍率ごとにエイリアシングノイズの電力と最大の成分を求めました。", "ほぼ線形の歪みを同じ方法で測り、計算の誤差で決まる下限（ノイズフロア）を求めました。", ], "results": { "aliasing_max_dbfs": {"value": round(c3["row"](sc["gain"], sc["note"], 1)["strongest_alias_dbfs"], 2), "decimals": 2}, }, "signal_seed": SEED, } csv_kinds = {"d1_buffer_ms.csv": "curve_csv", "d1_round_trip_ms.csv": "curve_csv", "d2_pickup_loading.csv": "frequency_response_csv", "d3_alias_power_table.csv": "curve_csv", "d3_alias_vs_oversampling.csv": "curve_csv", "d3_spectrum_g30_e6.csv": "spectrum_csv"} figures.write_all(OUT, results) rec = dc.write_outputs(HERE, OUT, EVIDENCE_ID, csv_kinds, detail, ledger) print(rec) if __name__ == "__main__": main()