import hashlib import json import math from pathlib import Path import subprocess import numpy as np import soundfile as sf from scipy.signal import butter, resample_poly, sosfilt from horn import SAMPLE_RATE, build_instrument, input_impedance, render_sequence from analysis import _yin_frames, analyze_tone OUTPUT = Path(__file__).resolve().parent ARCHIVE = OUTPUT.parent / "alto-sax-model-03" PROBES = ( ("steady-c4", 3.2, [{"time": 0.15, "duration": 2.5, "writtenMidi": 69}]), ("gentle-tongue", 3.2, [ {"time": 0.15, "duration": 0.85, "writtenMidi": 69}, *[{"time": 1.0 + index * 0.125, "duration": 0.125, "writtenMidi": 69} for index in range(12)], ]), ("connected-keys", 3.15, [ {"time": round(0.15 + index * 0.45, 8), "duration": 0.45, "writtenMidi": pitch, "tongued": index == 0} for index, pitch in enumerate((69, 71, 73, 71, 69)) ]), ("octave-switch", 3.6, [ {"time": round(0.15 + index * 0.9, 8), "duration": 0.9, "writtenMidi": pitch} for index, pitch in enumerate((69, 81, 69)) ]), ("pressure-swell", 3.8, [ {"time": 0.15, "duration": 2.9, "writtenMidi": 69, "pressurePa": 4800.0, "pressureEndPa": 5600.0, "gapM": 0.00072}, ]), ) def digest(path): return hashlib.sha256(path.read_bytes()).hexdigest() def rms(signal): return float(np.sqrt(np.mean(signal**2))) def write_audio(filename, signal): if not np.isfinite(signal).all() or np.max(np.abs(signal)) >= 1: raise ValueError(f"Invalid playback signal: {filename}") destination = OUTPUT / f"{filename}.wav" sf.write(destination, signal, SAMPLE_RATE, subtype="PCM_24") subprocess.run([ "ffmpeg", "-v", "error", "-y", "-i", str(destination), "-c:a", "libmp3lame", "-b:a", "192k", str(OUTPUT / f"{filename}.mp3"), ], check=True) def note_measurements(signal, events): measurements = [] for event in events: if event["duration"] < 0.3: continue start = event["time"] + min(0.3, event["duration"] * 0.4) end = event["time"] + event["duration"] - 0.04 frequency = 440 * 2 ** ((event["writtenMidi"] - 9 - 69) / 12) measured = analyze_tone(signal, frequency, start_seconds=start, end_seconds=end) measured.update({"writtenMidi": event["writtenMidi"], "concertMidi": event["writtenMidi"] - 9, "analysisWindowSeconds": [start, end]}) measurements.append(measured) return measurements def transient_metrics(signal, events): highpass = sosfilt(butter(4, 3000, fs=SAMPLE_RATE, btype="highpass", output="sos"), signal) block_size = round(0.01 * SAMPLE_RATE) steady_levels = [] for event in events: start = round((event["time"] + min(0.15, event["duration"] * 0.4)) * SAMPLE_RATE) end = round((event["time"] + event["duration"] - min(0.03, event["duration"] * 0.2)) * SAMPLE_RATE) blocks = [rms(signal[index:index + block_size]) for index in range(start, end - block_size + 1, block_size)] if blocks: steady_levels.append(float(np.median(blocks))) steady = max(float(np.median(steady_levels)), 1e-20) transitions = [] for event in events[1:]: onset = event["time"] start = max(0, round((onset - 0.04) * SAMPLE_RATE)) end = min(len(signal), round((onset + 0.06) * SAMPLE_RATE)) peak = float(np.max(np.abs(signal[start:end]))) high_blocks = [rms(highpass[index:index + block_size]) for index in range(start, end - block_size + 1, block_size)] transitions.append({ "timeSeconds": onset, "peakOverSteadyDb": 20 * math.log10(max(peak, 1e-20) / steady), "maximum10msHighpassRmsAtMinus20DbfsSteady": 20 * math.log10(max(max(high_blocks), 1e-20) * 0.1 / steady), }) return { "steadyMedianBlockRms": steady, "peakOverSteadyDb": 20 * math.log10(max(float(np.max(np.abs(signal))), 1e-20) / steady), "transitions": transitions, "limits": "Waveform and highpass proxies, not listening or loudness judgments; probes differ from the archived band performance.", } def short_note_measurements(signal, events): measurements = [] for event in events: if event["duration"] >= 0.3: continue start = event["time"] + 0.03 end = event["time"] + event["duration"] - 0.015 section = signal[round(start * SAMPLE_RATE):round(end * SAMPLE_RATE)] pitches, confidence = _yin_frames(resample_poly(section, 1, 4), SAMPLE_RATE / 4, 512, 128) voiced = confidence > 0.85 frequency = float(np.median(pitches[voiced])) if voiced.any() else None expected = 440 * 2 ** ((event["writtenMidi"] - 9 - 69) / 12) measurements.append({"timeSeconds": event["time"], "estimatedFrequencyHz": frequency, "pitchErrorCents": 1200 * math.log2(frequency / expected) if frequency else None, "voicedFraction": float(np.mean(voiced)), "analysisWindowSeconds": [start, end]}) return measurements def main(): archived_files = (ARCHIVE / "tone-c4.wav", ARCHIVE / "model.py", ARCHIVE / "geometry.json") archive_before = {path.name: digest(path) for path in archived_files} rendered = [] for filename, duration, events in PROBES: signal, details = render_sequence(events, duration) details.update({"file": f"{filename}.mp3", "durationSeconds": duration, "events": events, "notes": note_measurements(signal, events), "shortNotes": short_note_measurements(signal, events), "transients": transient_metrics(signal, events), "tailToPeakEnergyRatio": details["finalStoredEnergyJ"] / max(details["peakStoredEnergyJ"], 1e-20)}) if filename == "pressure-swell": details["swellWindows"] = [ dict(analyze_tone(signal, 261.625565, start_seconds=start, end_seconds=end), analysisWindowSeconds=[start, end]) for start, end in ((0.65, 1.15), (1.4, 1.9), (2.3, 2.85)) ] rendered.append((filename, signal, details)) print(filename, [(note["concertMidi"], note["pitchErrorCents"]) for note in details["notes"]], flush=True) common_gain = 0.8 / max(float(np.max(np.abs(signal))) for _, signal, _ in rendered) for filename, signal, _ in rendered: write_audio(filename, signal * common_gain) archived, archived_rate = sf.read(ARCHIVE / "tone-c4.wav") if archived_rate != SAMPLE_RATE or archived.ndim != 1: raise ValueError("Expected the archived mono 44.1 kHz C4") comparison = [] start, end = round(0.8 * SAMPLE_RATE), round(2.5 * SAMPLE_RATE) for filename, signal in (("ab-old", archived), ("ab-new", rendered[0][1])): gain = 0.1 / rms(signal[start:end]) comparison.append((filename, signal * gain, gain)) headroom = min(1.0, 0.85 / max(float(np.max(np.abs(signal))) for _, signal, _ in comparison)) matching = [] for filename, signal, gain in comparison: signal = signal * headroom write_audio(filename, signal) matching.append({"file": f"{filename}.mp3", "totalGain": gain * headroom, "sustainRmsBeforeEncoding": rms(signal[start:end]), "peakBeforeEncoding": float(np.max(np.abs(signal)))}) network = build_instrument() frequencies = np.geomspace(80, 16000, 512) minimum_real = min(float(np.min(np.real(input_impedance(network, pitch, frequencies)))) for pitch in range(58, 86)) quiet, quiet_details = render_sequence( [{"time": 0.05, "duration": 0.3, "writtenMidi": 69}], 0.55, breath_scale=0) unforced, unforced_details = render_sequence( [{"time": 0.1, "duration": 1.0, "writtenMidi": 69}], 1.8, tongue_strength=0) unforced_pitch = analyze_tone(unforced, 261.625565, start_seconds=0.5, end_seconds=1.05) notes = [note for _, _, details in rendered for note in details["notes"]] short_notes = [note for _, _, details in rendered for note in details["shortNotes"]] swell_windows = rendered[-1][2]["swellWindows"] checks = { "zeroBreathIsSilent": not bool(quiet.any()), "zeroBreathStoredEnergyIsZero": quiet_details["peakStoredEnergyJ"] == 0, "positiveRealPassiveInputImpedance": minimum_real > 0, "minimumRealImpedancePaSecondsPerM3": minimum_real, "maximumCouplingResidualPa": max(details["maximumCouplingResidualPa"] for _, _, details in rendered), "couplingResidualBelow0Point01Pa": all(details["maximumCouplingResidualPa"] < 0.01 for _, _, details in rendered), "testedNotesWithin25Cents": all(note["pitchErrorCents"] is not None and abs(note["pitchErrorCents"]) < 25 for note in notes), "testedNotesVoicedAbove95Percent": all(note["voicedFraction"] > 0.95 for note in notes), "eightHzTonguedNotesStayPitched": len(short_notes) == 12 and all( note["pitchErrorCents"] is not None and abs(note["pitchErrorCents"]) < 25 and note["voicedFraction"] > 0.95 for note in short_notes), "swellStaysOnC4InThreeWindows": all(note["pitchErrorCents"] is not None and abs(note["pitchErrorCents"]) < 25 and note["voicedFraction"] > 0.95 for note in swell_windows), "energyDecaysAfterBreath": all(details["tailToPeakEnergyRatio"] < 0.02 for _, _, details in rendered), "noAcousticStateResets": all(details["keyStateResets"] == 0 for _, _, details in rendered), "c4OscillatesWithoutTongueForce": unforced_pitch["pitchErrorCents"] is not None and abs(unforced_pitch["pitchErrorCents"]) < 25, "archiveHashesUnchanged": archive_before == {path.name: digest(path) for path in archived_files}, "comparisonSustainRmsDifference": abs(matching[0]["sustainRmsBeforeEncoding"] - matching[1]["sustainRmsBeforeEncoding"]), } report = { "model": "Continuous reed, bore, chimney pockets and finite key apertures with compliant tongue contact", "status": "Unapproved standalone instrument prototype; no listening or realism claim", "recordedSamplesUsed": False, "pitchOscillatorsUsed": False, "audioPitchShiftingUsed": False, "sampleRate": SAMPLE_RATE, "commonProbePlaybackGain": common_gain, "sourceSha256": {name: digest(OUTPUT / name) for name in ("horn.py", "render_tests.py")}, "archiveSha256": archive_before, "comparison": {"method": "Identical 0.8-2.5 s unweighted sustain RMS, then shared peak headroom; no compression or limiting", "source": "../alto-sax-model-03/tone-c4.wav", "sharedHeadroomGain": headroom, "files": matching, "limits": "RMS matching is not perceptual loudness matching; MP3 encoding may change levels slightly"}, "checks": checks, "tests": [details for _, _, details in rendered], "tongueForceDisabledCheck": {"pitch": unforced_pitch, "engine": unforced_details}, "limits": [ "Designed generic geometry and provisional physical coefficients, not a measured alto saxophone", "Key aperture motion preserves states but prescribed key-muscle work is not included in the energy trace", "Tongue is one compliant closing-contact approximation, not a complete measured human articulation model", "Register changes also use prescribed breath and slow embouchure adjustment", "Mixed trapezoidal and backward-Euler lumped acoustics, reference-frequency wall loss, first-order radiation", "No upstream vocal tract, detailed bell directivity, or external tonehole interaction", "Pitch checks cover only these probes, not the full instrument range or fast music", "Earlier-merge diagnostics identify correlated spikes, not measured attribution to a single cause", "Numerical checks do not establish that this sounds better or realistic", ], } (OUTPUT / "report.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n") print(json.dumps(checks, indent=2), flush=True) failures = [name for name, value in checks.items() if isinstance(value, bool) and not value] if failures: raise AssertionError(f"Instrument checks failed: {failures}") if __name__ == "__main__": main()