import math from pathlib import Path import sys import numpy as np from numba import njit from scipy.signal import resample_poly SOURCE = Path(__file__).resolve().parent.parent / "alto-sax-model-03" sys.path.insert(0, str(SOURCE)) import model as physics SAMPLE_RATE = physics.SAMPLE_RATE MODEL_RATE = physics.MODEL_SAMPLE_RATE def build_instrument(vent_radius=0.00065): geometry = physics.load_geometry() network = physics.build_network(67, geometry) old_positions = network["positions"][network["branch_nodes"]] nodes = np.unique(np.append(network["positions"], (0.090, 0.235))) edges = np.concatenate(([0.0], (nodes[:-1] + nodes[1:]) * 0.5, [geometry.length_m])) left = physics.radius_at(edges[:-1], geometry) right = physics.radius_at(edges[1:], geometry) volumes = math.pi * np.diff(edges) * (left**2 + left * right + right**2) / 3 volumes[0] += physics.missing_cone_volume(geometry, physics.Reed()) radii = physics.radius_at(nodes, geometry) lengths = np.diff(nodes) inertance = physics.AIR_DENSITY * lengths / (math.pi * radii[:-1] * radii[1:]) inertance[0] += physics.AIR_DENSITY * geometry.mouthpiece_throat_m / (math.pi * geometry.neck_radius_m**2) edge_radius = np.sqrt(radii[:-1] * radii[1:]) viscous = np.maximum(1.0, edge_radius * np.sqrt(2 * math.pi * 1000 * physics.AIR_DENSITY / physics.AIR_VISCOSITY) / 4) network["positions"] = nodes network["compliance"] = volumes / (physics.AIR_DENSITY * physics.SPEED_OF_SOUND**2) network["inertance"] = inertance network["resistance"] = 8 * physics.AIR_VISCOSITY * lengths / (math.pi * edge_radius**4) * viscous network["branch_nodes"] = np.searchsorted(nodes, np.insert(old_positions, 16, 0.090)) for name in ("branch_inertance", "branch_resistance", "load_inertance", "load_resistance", "load_compliance"): network[name] = np.insert(network[name], 16, network[name][15]) network["branch_keys"] = np.insert(network["branch_keys"], 16, 0) area = math.pi * vent_radius**2 characteristic = physics.AIR_DENSITY * physics.SPEED_OF_SOUND / area viscous = max(1.0, vent_radius * math.sqrt(2 * math.pi * 1000 * physics.AIR_DENSITY / physics.AIR_VISCOSITY) / 4) network["branch_inertance"][15:17] = physics.AIR_DENSITY * (0.003 + 0.82 * vent_radius) / area network["branch_resistance"][15:17] = 8 * physics.AIR_VISCOSITY * 0.003 / (math.pi * vent_radius**4) * viscous network["load_inertance"][15:17] = characteristic * 0.61 * vent_radius / physics.SPEED_OF_SOUND network["load_resistance"][15:17] = characteristic * 4 * 0.61**2 network["load_compliance"][15:17] = area * 0.003 / (physics.AIR_DENSITY * physics.SPEED_OF_SOUND**2) network["valve_inertance_scale"] = network["load_inertance"] * 0.5 network["valve_resistance_scale"] = network["load_resistance"] * 0.015 return network def key_apertures(network, written_midi): if not 58 <= written_midi <= 85: raise ValueError("The generic horn covers written Bb3 through C#6") base = written_midi if written_midi <= 73 else written_midi - 12 apertures = (network["branch_keys"] <= base).astype(np.float64) apertures[15:17] = 0.0 if written_midi > 73: apertures[15 if base < 69 else 16] = 1.0 apertures[-1] = 1.0 return apertures def input_impedance(network, written_midi, frequencies): frequency = 2j * math.pi * np.atleast_1d(frequencies) edge_y = 1 / (network["resistance"][:, None] + network["inertance"][:, None] * frequency) diagonal = network["compliance"][:, None] * frequency diagonal[:-1] += edge_y diagonal[1:] += edge_y for branch, aperture in enumerate(key_apertures(network, written_midi)): ratio = 1 / max(1e-6, aperture) - 1 valve = network["valve_resistance_scale"][branch] * ratio**2 + network["valve_inertance_scale"][branch] * max(1e-6, ratio) * frequency radiation = 1 / (1 / network["load_resistance"][branch] + 1 / (network["load_inertance"][branch] * frequency)) load_y = network["load_compliance"][branch] * frequency + 1 / (valve + radiation) chimney = network["branch_resistance"][branch] + network["branch_inertance"][branch] * frequency diagonal[network["branch_nodes"][branch]] += 1 / (chimney + 1 / load_y) schur = diagonal[-1].copy() for node in range(len(diagonal) - 2, -1, -1): schur = diagonal[node] - edge_y[node]**2 / schur return 1 / schur @njit(cache=True) def _simulate(noise, starts, ends, key_starts, pressures, pressure_ends, gaps, tongued, targets, compliance, inertance, resistance, branch_nodes, neck_inertance, neck_resistance, radiation_inertance, radiation_resistance, pocket_compliance, valve_scale, valve_loss, reed_parameters, tongue_strength, key_seconds): timestep = 1 / MODEL_RATE node_count = len(compliance) branch_count = len(branch_nodes) capacitor_g = 2 * compliance / timestep edge_g = 1 / (2 * inertance / timestep + resistance) neck_g = 1 / (2 * neck_inertance / timestep + neck_resistance) pocket_g = 2 * pocket_compliance / timestep radiation_g = 1 / radiation_resistance + timestep / (2 * radiation_inertance) capacitor_history = np.zeros(node_count) edge_history = np.zeros(node_count - 1) neck_history = np.zeros(branch_count) pocket_history = np.zeros(branch_count) radiation_history = np.zeros(branch_count) valve_flux = np.zeros(branch_count) apertures = targets[0].copy() valve_g = np.zeros(branch_count) valve_l = np.zeros(branch_count) outlet_g = np.zeros(branch_count) outlet_history = np.zeros(branch_count) load_g = np.zeros(branch_count) load_history = np.zeros(branch_count) branch_g = np.zeros(branch_count) branch_history = np.zeros(branch_count) pressure = np.zeros(node_count) inverse_diagonal = np.zeros(node_count) multipliers = np.zeros(node_count - 1) port_source = np.zeros(node_count) port_source[0] = 1 audio = np.zeros(len(noise)) energy = np.zeros((len(noise) - 1) // 882 + 1) mass, stiffness, damping, width, reed_area, contact_stiffness = reed_parameters denominator = 4 * mass / timestep**2 + 2 * damping / timestep + stiffness reed_position = 0.0 reed_velocity = 0.0 reed_force = 0.0 mouth_pressure = 0.0 lip_gap = gaps[0] tongue_contact = 0.0 breath_state = 0.0 last_radiating_flow = 0.0 radiation_state = 0.0 highpass_state = 0.0 previous_radiation = 0.0 maximum_residual = 0.0 peak_bore_pressure = 0.0 peak_radiating_flow = 0.0 event = 0 key_event = 0 key_alpha = 1 - math.exp(-timestep / key_seconds) pressure_alpha = 1 - math.exp(-timestep / 0.012) lip_alpha = 1 - math.exp(-timestep / 0.025) tongue_alpha = 1 - math.exp(-timestep / 0.0015) breath_alpha = 1 - math.exp(-2 * math.pi * 160 * timestep) microphone_alpha = 1 - math.exp(-2 * math.pi * 7200 * timestep) highpass_alpha = math.exp(-2 * math.pi * 35 * timestep) for sample in range(len(noise)): if event + 1 < len(starts) and sample >= starts[event + 1]: event += 1 if key_event + 1 < len(starts) and sample >= key_starts[key_event + 1]: key_event += 1 desired_pressure = 0.0 desired_gap = gaps[event] desired_tongue = 0.0 active = starts[event] <= sample < ends[event] if active: progress = (sample - starts[event]) / max(1, ends[event] - starts[event]) desired_pressure = pressures[event] + (pressure_ends[event] - pressures[event]) * progress desired_gap *= desired_pressure / max(1.0, pressures[event]) elif sample < starts[0]: if (starts[0] - sample) * timestep < 0.060: desired_pressure = pressures[0] desired_tongue = 1.0 else: since_end = (sample - ends[event]) * timestep if since_end < 0.080: desired_pressure = pressure_ends[event] * math.cos(0.5 * math.pi * since_end / 0.080)**2 if event + 1 < len(starts) and (starts[event + 1] - sample) * timestep < 0.060: desired_pressure = pressures[event + 1] desired_gap = gaps[event + 1] desired_tongue = 1.0 elapsed = (sample - starts[event]) * timestep if active and tongued[event] and elapsed < 0.008: desired_tongue = math.cos(0.5 * math.pi * elapsed / 0.008)**2 if event + 1 < len(starts) and tongued[event + 1]: to_next = (starts[event + 1] - sample) * timestep if 0 < to_next < 0.018: desired_tongue = max(desired_tongue, math.sin(0.5 * math.pi * min(1.0, (0.018 - to_next) / 0.012))**2) if desired_pressure <= 0: desired_tongue = 0.0 mouth_pressure += pressure_alpha * (desired_pressure - mouth_pressure) lip_gap += lip_alpha * (desired_gap - lip_gap) tongue_contact += tongue_alpha * (desired_tongue * tongue_strength - tongue_contact) breath_state += breath_alpha * (noise[sample] - breath_state) played_pressure = mouth_pressure * (1 + 0.0015 * breath_state) diagonal = capacitor_g.copy() diagonal[:-1] += edge_g diagonal[1:] += edge_g right_hand = -capacitor_history.copy() for edge in range(node_count - 1): right_hand[edge] -= edge_history[edge] right_hand[edge + 1] += edge_history[edge] for branch in range(branch_count): apertures[branch] += key_alpha * (targets[key_event, branch] - apertures[branch]) ratio = 1 / max(1e-6, apertures[branch]) - 1 valve_l[branch] = valve_scale[branch] * max(1e-6, ratio) valve_r = valve_loss[branch] * ratio**2 valve_g[branch] = 1 / (valve_l[branch] / timestep + valve_r) valve_history = valve_g[branch] * valve_flux[branch] / timestep divisor = 1 + valve_g[branch] / radiation_g[branch] outlet_g[branch] = valve_g[branch] / divisor outlet_history[branch] = (valve_history + valve_g[branch] / radiation_g[branch] * radiation_history[branch]) / divisor load_g[branch] = pocket_g[branch] + outlet_g[branch] load_history[branch] = pocket_history[branch] + outlet_history[branch] divisor = 1 + neck_g[branch] / load_g[branch] branch_g[branch] = neck_g[branch] / divisor branch_history[branch] = (neck_history[branch] + neck_g[branch] / load_g[branch] * load_history[branch]) / divisor diagonal[branch_nodes[branch]] += branch_g[branch] right_hand[branch_nodes[branch]] -= branch_history[branch] for node in range(1, node_count): multipliers[node - 1] = -edge_g[node - 1] / diagonal[node - 1] diagonal[node] += multipliers[node - 1] * edge_g[node - 1] for node in range(node_count): inverse_diagonal[node] = 1 / diagonal[node] port_response = physics._solve_tridiagonal(port_source, inverse_diagonal, edge_g, multipliers) free_pressure = physics._solve_tridiagonal(right_hand, inverse_diagonal, edge_g, multipliers) position_history = ((4 * mass / timestep**2 + 2 * damping / timestep - stiffness) * reed_position + 4 * mass / timestep * reed_velocity + reed_force) / denominator bore_pressure = free_pressure[0] source_flow = 0.0 new_position = reed_position new_velocity = reed_velocity tongue_force = 0.0 for iteration in range(48): difference = played_pressure - bore_pressure new_position = position_history + reed_area * difference / denominator for collision in range(5): penetration = max(0.0, new_position - lip_gap) contact = contact_stiffness * penetration**1.5 contact_derivative = 1.5 * contact_stiffness * math.sqrt(penetration) new_velocity = 2 / timestep * (new_position - reed_position) - reed_velocity tongue_penetration = max(0.0, lip_gap + 0.000010 - new_position) tongue_force = tongue_contact * (20000 * tongue_penetration - 0.12 * new_velocity) if tongue_penetration > 0 else 0.0 tongue_derivative = tongue_contact * (20000 + 0.12 * 2 / timestep) if tongue_penetration > 0 else 0.0 residual = denominator * (new_position - position_history) - reed_area * difference + contact - tongue_force new_position -= residual / (denominator + contact_derivative + tongue_derivative) penetration = max(0.0, new_position - lip_gap) contact_derivative = 1.5 * contact_stiffness * math.sqrt(penetration) new_velocity = 2 / timestep * (new_position - reed_position) - reed_velocity tongue_penetration = max(0.0, lip_gap + 0.000010 - new_position) tongue_force = tongue_contact * (20000 * tongue_penetration - 0.12 * new_velocity) if tongue_penetration > 0 else 0.0 tongue_derivative = tongue_contact * (20000 + 0.12 * 2 / timestep) if tongue_penetration > 0 else 0.0 position_derivative = -reed_area / (denominator + contact_derivative + tongue_derivative) opening = max(0.0, lip_gap - new_position) magnitude = abs(difference) root = math.sqrt(2 / physics.AIR_DENSITY) * magnitude / math.sqrt(magnitude + 5) direction = 1.0 if difference >= 0 else -1.0 jet_factor = 1 + 0.0015 * noise[sample] source_flow = direction * width * opening * root * jet_factor + reed_area * new_velocity derivative_flow = -width * opening * math.sqrt(2 / physics.AIR_DENSITY) * (0.5 * magnitude + 5) / (magnitude + 5)**1.5 * jet_factor if opening > 0: derivative_flow -= direction * width * root * position_derivative * jet_factor derivative_flow += reed_area * 2 / timestep * position_derivative residual = bore_pressure - free_pressure[0] - port_response[0] * source_flow if abs(residual) < 1e-6: break derivative = 1 - port_response[0] * derivative_flow bore_pressure -= max(-2000.0, min(2000.0, residual / max(0.25, derivative))) maximum_residual = max(maximum_residual, abs(bore_pressure - free_pressure[0] - port_response[0] * source_flow)) reed_force = reed_area * (played_pressure - bore_pressure) - contact_stiffness * max(0.0, new_position - lip_gap)**1.5 + tongue_force reed_position = new_position reed_velocity = new_velocity stored = 0.5 * mass * reed_velocity**2 + 0.5 * stiffness * reed_position**2 for node in range(node_count): pressure[node] = free_pressure[node] + port_response[node] * source_flow current = capacitor_g[node] * pressure[node] + capacitor_history[node] capacitor_history[node] = -capacitor_g[node] * pressure[node] - current stored += 0.5 * compliance[node] * pressure[node]**2 peak_bore_pressure = max(peak_bore_pressure, abs(pressure[0])) for edge in range(node_count - 1): difference = pressure[edge] - pressure[edge + 1] flow = edge_g[edge] * difference + edge_history[edge] edge_history[edge] = edge_g[edge] * ((2 * inertance[edge] / timestep - resistance[edge]) * flow + difference) stored += 0.5 * inertance[edge] * flow**2 radiating_flow = 0.0 for branch in range(branch_count): flow = branch_g[branch] * pressure[branch_nodes[branch]] + branch_history[branch] pocket_pressure = (flow - load_history[branch]) / load_g[branch] neck_pressure = pressure[branch_nodes[branch]] - pocket_pressure neck_history[branch] = neck_g[branch] * ((2 * neck_inertance[branch] / timestep - neck_resistance[branch]) * flow + neck_pressure) pocket_current = pocket_g[branch] * pocket_pressure + pocket_history[branch] pocket_history[branch] = -pocket_g[branch] * pocket_pressure - pocket_current emitted_flow = outlet_g[branch] * pocket_pressure + outlet_history[branch] external_pressure = (emitted_flow - radiation_history[branch]) / radiation_g[branch] radiation_flow = emitted_flow - external_pressure / radiation_resistance[branch] radiation_history[branch] = radiation_flow + timestep / (2 * radiation_inertance[branch]) * external_pressure valve_flux[branch] = valve_l[branch] * emitted_flow radiating_flow += emitted_flow stored += 0.5 * neck_inertance[branch] * flow**2 + 0.5 * pocket_compliance[branch] * pocket_pressure**2 stored += 0.5 * radiation_inertance[branch] * radiation_flow**2 + 0.5 * valve_l[branch] * emitted_flow**2 peak_radiating_flow = max(peak_radiating_flow, abs(radiating_flow)) radiation = physics.AIR_DENSITY / (4 * math.pi) * (radiating_flow - last_radiating_flow) / timestep last_radiating_flow = radiating_flow radiation_state += microphone_alpha * (radiation - radiation_state) highpass_state = highpass_alpha * (highpass_state + radiation_state - previous_radiation) previous_radiation = radiation_state audio[sample] = highpass_state if sample % 882 == 0: energy[sample // 882] = stored return audio, energy, maximum_residual, peak_bore_pressure, peak_radiating_flow def render_sequence(events, duration, seed=20261011, breath_scale=1.0, tongue_strength=1.0, vent_radius=0.00065, key_seconds=0.006): if not events or duration <= 0 or breath_scale < 0 or not 0 <= tongue_strength <= 1 or key_seconds <= 0: raise ValueError("A valid physical playing sequence is required") if any(right["time"] <= left["time"] for left, right in zip(events, events[1:])): raise ValueError("Events must be strictly time ordered") if any(event["duration"] <= 0 or event["time"] < 0 or event["time"] + event["duration"] > duration for event in events): raise ValueError("Every positive-duration event must fit inside the render") network = build_instrument(vent_radius) written = [event["writtenMidi"] for event in events] starts = np.asarray([round(event["time"] * MODEL_RATE) for event in events], dtype=np.int64) ends = np.asarray([round((event["time"] + event["duration"]) * MODEL_RATE) for event in events], dtype=np.int64) if np.any(ends[:-1] > starts[1:]): raise ValueError("The instrument is monophonic") pressure = np.asarray([event.get("pressurePa", 5500 if pitch in (74, 75) else 6500 if pitch > 73 else 5000) for event, pitch in zip(events, written)], dtype=np.float64) * breath_scale pressure_ends = np.asarray([event.get("pressureEndPa", initial / max(breath_scale, 1e-20)) for event, initial in zip(events, pressure)]) * breath_scale gaps = np.asarray([event.get("gapM", 0.00060 if pitch in (74, 75) else 0.00050 if pitch == 76 else 0.00045 if pitch > 73 else 0.00075) for event, pitch in zip(events, written)]) tongued = np.asarray([event.get("tongued", True) for event in events], dtype=np.bool_) key_starts = starts.copy() for index in range(1, len(events)): if tongued[index]: key_starts[index] = max(ends[index - 1], starts[index] - round(0.012 * MODEL_RATE)) reed = physics.Reed() damping = 2 * reed.damping_ratio * math.sqrt(reed.mass_kg * reed.stiffness_n_m) parameters = np.asarray((reed.mass_kg, reed.stiffness_n_m, damping, reed.width_m, reed.effective_area_m2, reed.contact_stiffness)) noise = np.random.default_rng(seed).normal(size=round(duration * MODEL_RATE)) audio, energy, residual, peak, peak_flow = _simulate(noise, starts, ends, key_starts, pressure, pressure_ends, gaps, tongued, np.asarray([key_apertures(network, pitch) for pitch in written]), network["compliance"], network["inertance"], network["resistance"], network["branch_nodes"], network["branch_inertance"], network["branch_resistance"], network["load_inertance"], network["load_resistance"], network["load_compliance"], network["valve_inertance_scale"], network["valve_resistance_scale"], parameters, tongue_strength, key_seconds) if not np.isfinite(audio).all() or not np.isfinite(energy).all() or residual > 0.01: raise FloatingPointError(f"The coupled horn solve failed: {residual} Pa") return resample_poly(audio, 1, physics.OVERSAMPLE), { "maximumCouplingResidualPa": float(residual), "peakBorePressurePa": float(peak), "peakRadiatingVolumeFlowM3PerSecond": float(peak_flow), "peakStoredEnergyJ": float(np.max(energy)), "finalStoredEnergyJ": float(energy[-1]), "energyTraceJ": energy.tolist(), "modelSampleRate": MODEL_RATE, "registerVentRadiusMm": vent_radius * 1000, "keyTimeConstantMs": key_seconds * 1000, "keyStateResets": 0, "tongueModel": "Compliant force and contact damping on the reed; no tongue-driven mouthpiece-gap changes.", "radiationModel": "Actual outlet flow after the key aperture; no crossfaded flow-summing masks.", "limits": "Generic geometry; mixed trapezoidal/backward-Euler lumped acoustics; prescribed key motion does mechanical work not included in the energy trace.", }