from dataclasses import dataclass, replace import json import math from pathlib import Path import numpy as np from numba import njit from scipy.signal import resample_poly SAMPLE_RATE = 44100 OVERSAMPLE = 2 MODEL_SAMPLE_RATE = SAMPLE_RATE * OVERSAMPLE AIR_DENSITY = 1.20 SPEED_OF_SOUND = 343.0 AIR_VISCOSITY = 1.84e-5 @dataclass(frozen=True) class Reed: mass_kg: float = 1.2e-5 stiffness_n_m: float = 1300.0 damping_ratio: float = 0.22 width_m: float = 0.014 gap_m: float = 0.00075 effective_area_m2: float = 8.5e-5 contact_stiffness: float = 3.0e7 @dataclass(frozen=True) class Player: pressure_pa: float = 5000.0 attack_seconds: float = 0.075 release_seconds: float = 0.12 tail_seconds: float = 0.30 vibrato_cents: float = 0.0 vibrato_rate_hz: float = 5.3 vibrato_start_seconds: float = 0.8 turbulence: float = 0.0015 tongue_rate_hz: float = 0.0 tongue_strength: float = 0.88 crescendo: float = 0.0 opening_ratio: float = 1.0 @dataclass(frozen=True) class Geometry: length_m: float = 1.00 neck_radius_m: float = 0.006 cone_slope: float = 0.028 bell_flare_m: float = 0.018 mouthpiece_scale: float = 1.0 mouthpiece_throat_m: float = 0.046 mesh_spacing_m: float = 0.010 hole_positions: tuple = () def midi_frequency(midi_note): return 440.0 * 2.0 ** ((midi_note - 69) / 12.0) def note_name(midi_note): names = ("C", "C#", "D", "Eb", "E", "F", "F#", "G", "G#", "A", "Bb", "B") return f"{names[midi_note % 12]}{midi_note // 12 - 1}" def load_geometry(): path = Path(__file__).resolve().parent / "geometry.json" if not path.exists(): return Geometry() parameters = json.loads(path.read_text())["geometry"] parameters["hole_positions"] = tuple(parameters["hole_positions"]) return Geometry(**parameters) def radius_at(position, geometry): fraction = np.maximum(0.0, (np.asarray(position) / geometry.length_m - 0.84) / 0.16) return geometry.neck_radius_m + geometry.cone_slope * np.asarray(position) + geometry.bell_flare_m * fraction**3 def missing_cone_volume(geometry, reed): virtual_length = geometry.neck_radius_m / geometry.cone_slope volume = math.pi * geometry.neck_radius_m**2 * virtual_length / 3.0 reed_volume = AIR_DENSITY * SPEED_OF_SOUND**2 * reed.effective_area_m2**2 / reed.stiffness_n_m return volume * geometry.mouthpiece_scale - reed_volume def build_network(written_midi, geometry=None, reed=None, register_vent=True): geometry = geometry or Geometry() reed = reed or Reed() if not 58 <= written_midi <= 85: raise ValueError("This instrument study covers written Bb3 through C#6") if missing_cone_volume(geometry, reed) <= 0.0: raise ValueError("Mouthpiece compliance must remain positive") base_note = written_midi if written_midi <= 73 else written_midi - 12 positions = geometry.hole_positions or tuple( min(geometry.length_m - 0.02, max(0.10, SPEED_OF_SOUND / (2.0 * midi_frequency(note - 9)) - geometry.neck_radius_m / geometry.cone_slope - 0.030)) for note in range(59, 74) ) if len(positions) != 15 or any(position <= 0.0 or position >= geometry.length_m for position in positions): raise ValueError("Fifteen positive tone-hole positions are required inside the bore") holes = [] for note, position in zip(range(59, 74), positions): radius = 0.0068 + 0.0040 * (position / geometry.length_m) holes.append((position, radius, 0.0045, note <= base_note, note)) vent_position = 0.235 if base_note < 69 else 0.090 holes.append((vent_position, 0.00065, 0.0030, written_midi > 73 and register_vent, 0)) nodes = np.unique(np.concatenate(( np.linspace(0.0, geometry.length_m, math.ceil(geometry.length_m / geometry.mesh_spacing_m) + 1), np.asarray([hole[0] for hole in holes]), ))) cell_edges = np.concatenate(([0.0], (nodes[:-1] + nodes[1:]) * 0.5, [geometry.length_m])) left_radii = radius_at(cell_edges[:-1], geometry) right_radii = radius_at(cell_edges[1:], geometry) volumes = math.pi * np.diff(cell_edges) * (left_radii**2 + left_radii * right_radii + right_radii**2) / 3.0 volumes[0] += missing_cone_volume(geometry, reed) compliance = volumes / (AIR_DENSITY * SPEED_OF_SOUND**2) radii = radius_at(nodes, geometry) lengths = np.diff(nodes) inertance = AIR_DENSITY * lengths / (math.pi * radii[:-1] * radii[1:]) inertance[0] += AIR_DENSITY * geometry.mouthpiece_throat_m / (math.pi * geometry.neck_radius_m**2) edge_radius = np.sqrt(radii[:-1] * radii[1:]) viscous_factor = np.maximum(1.0, edge_radius * np.sqrt(2.0 * math.pi * 1000.0 * AIR_DENSITY / AIR_VISCOSITY) / 4.0) resistance = 8.0 * AIR_VISCOSITY * lengths / (math.pi * edge_radius**4) * viscous_factor branch_nodes = [] branch_inertance = [] branch_resistance = [] load_inertance = [] load_resistance = [] load_compliance = [] branch_open = [] branch_keys = [] for position, radius, chimney, is_open, key in holes: area = math.pi * radius**2 characteristic = AIR_DENSITY * SPEED_OF_SOUND / area branch_nodes.append(int(np.searchsorted(nodes, position))) branch_inertance.append(AIR_DENSITY * (chimney + 0.82 * radius) / area) branch_loss = max(1.0, radius * math.sqrt(2.0 * math.pi * 1000.0 * AIR_DENSITY / AIR_VISCOSITY) / 4.0) branch_resistance.append(8.0 * AIR_VISCOSITY * chimney / (math.pi * radius**4) * branch_loss) load_inertance.append(characteristic * 0.61 * radius / SPEED_OF_SOUND) load_resistance.append(characteristic * 4.0 * 0.61**2) load_compliance.append(area * chimney / (AIR_DENSITY * SPEED_OF_SOUND**2)) branch_open.append(is_open) branch_keys.append(key) radius = float(radii[-1]) characteristic = AIR_DENSITY * SPEED_OF_SOUND / (math.pi * radius**2) branch_nodes.append(len(nodes) - 1) branch_inertance.append(AIR_DENSITY * geometry.mesh_spacing_m * 0.25 / (math.pi * radius**2)) branch_resistance.append(0.0) load_inertance.append(characteristic * 0.61 * radius / SPEED_OF_SOUND) load_resistance.append(characteristic * 4.0 * 0.61**2) load_compliance.append(0.0) branch_open.append(True) branch_keys.append(-1) return { "positions": nodes, "compliance": compliance, "inertance": inertance, "resistance": resistance, "branch_nodes": np.asarray(branch_nodes, dtype=np.int64), "branch_inertance": np.asarray(branch_inertance), "branch_resistance": np.asarray(branch_resistance), "load_inertance": np.asarray(load_inertance), "load_resistance": np.asarray(load_resistance), "load_compliance": np.asarray(load_compliance), "branch_open": np.asarray(branch_open, dtype=np.bool_), "branch_keys": np.asarray(branch_keys, dtype=np.int64), } def input_impedance(network, frequencies, reed=None): frequencies = np.atleast_1d(frequencies) complex_frequency = 2j * math.pi * frequencies edge_admittance = 1.0 / (network["resistance"][:, None] + network["inertance"][:, None] * complex_frequency) diagonal = network["compliance"][:, None] * complex_frequency diagonal[:-1] += edge_admittance diagonal[1:] += edge_admittance for branch, node in enumerate(network["branch_nodes"]): if network["branch_open"][branch]: load = 1.0 / (1.0 / network["load_resistance"][branch] + 1.0 / (network["load_inertance"][branch] * complex_frequency)) else: load = 1.0 / (network["load_compliance"][branch] * complex_frequency) diagonal[node] += 1.0 / (network["branch_resistance"][branch] + network["branch_inertance"][branch] * complex_frequency + load) if reed is not None: damping = 2.0 * reed.damping_ratio * math.sqrt(reed.mass_kg * reed.stiffness_n_m) diagonal[0] += reed.effective_area_m2**2 * complex_frequency / ( reed.mass_kg * complex_frequency**2 + damping * complex_frequency + reed.stiffness_n_m ) schur = diagonal[-1].copy() for node in range(len(diagonal) - 2, -1, -1): schur = diagonal[node] - edge_admittance[node]**2 / schur return 1.0 / schur @njit(cache=True) def _solve_tridiagonal(right_hand, inverse_diagonal, edge_g, multipliers): result = right_hand.copy() for node in range(1, len(result)): result[node] -= multipliers[node - 1] * result[node - 1] result[-1] *= inverse_diagonal[-1] for node in range(len(result) - 2, -1, -1): result[node] = (result[node] + edge_g[node] * result[node + 1]) * inverse_diagonal[node] return result def prepare_solver(network): timestep = 1.0 / MODEL_SAMPLE_RATE capacitor_g = 2.0 * network["compliance"] / timestep edge_g = 1.0 / (2.0 * network["inertance"] / timestep + network["resistance"]) neck_g = 1.0 / (2.0 * network["branch_inertance"] / timestep + network["branch_resistance"]) load_g = np.where(network["branch_open"], 1.0 / network["load_resistance"] + timestep / (2.0 * network["load_inertance"]), 2.0 * network["load_compliance"] / timestep, ) branch_g = neck_g / (1.0 + neck_g / load_g) diagonal = capacitor_g.copy() diagonal[:-1] += edge_g diagonal[1:] += edge_g np.add.at(diagonal, network["branch_nodes"], branch_g) multipliers = np.zeros(len(edge_g)) for node in range(1, len(diagonal)): multipliers[node - 1] = -edge_g[node - 1] / diagonal[node - 1] diagonal[node] += multipliers[node - 1] * edge_g[node - 1] inverse_diagonal = 1.0 / diagonal port_source = np.zeros(len(diagonal)) port_source[0] = 1.0 port_response = _solve_tridiagonal(port_source, inverse_diagonal, edge_g, multipliers) return capacitor_g, edge_g, neck_g, load_g, branch_g, multipliers, inverse_diagonal, port_response @njit(cache=True) def _simulate(duration, pressure_target, attack, release, vibrato_cents, vibrato_rate, vibrato_start, turbulence, tongue_rate, tongue_strength, crescendo, noise, reed_parameters, compliance, inertance, resistance, branch_nodes, branch_inertance, branch_resistance, load_inertance, load_resistance, load_compliance, branch_open, capacitor_g, edge_g, neck_g, load_g, branch_g, multipliers, inverse_diagonal, port_response): timestep = 1.0 / MODEL_SAMPLE_RATE count = len(noise) node_count = len(compliance) branch_count = len(branch_nodes) capacitor_history = np.zeros(node_count) edge_history = np.zeros(node_count - 1) neck_history = np.zeros(branch_count) load_history = np.zeros(branch_count) branch_history = np.zeros(branch_count) audio = np.zeros(count) energy = np.zeros((count - 1) // 882 + 1) pressure = np.zeros(node_count) mass, stiffness, damping, width, rest_gap, reed_area, contact_stiffness = reed_parameters reed_denominator = 4.0 * mass / timestep**2 + 2.0 * damping / timestep + stiffness reed_position = 0.0 reed_velocity = 0.0 reed_force = 0.0 mouth_pressure = 0.0 breath_state = 0.0 last_output_flow = 0.0 radiation_state = 0.0 highpass_state = 0.0 previous_radiation = 0.0 peak_bore_pressure = 0.0 maximum_source_residual = 0.0 contact_count = 0 pressure_alpha = 1.0 - math.exp(-timestep / 0.008) breath_alpha = 1.0 - math.exp(-2.0 * math.pi * 160.0 * timestep) radiation_alpha = 1.0 - math.exp(-2.0 * math.pi * 7200.0 * timestep) highpass_alpha = math.exp(-2.0 * math.pi * 35.0 * timestep) for sample in range(count): time = sample * timestep if time < attack: envelope = math.sin(0.5 * math.pi * time / attack)**2 elif time > duration - release: envelope = math.cos(0.5 * math.pi * min(1.0, (time - duration + release) / release))**2 else: envelope = 1.0 desired_pressure = pressure_target * envelope * (1.0 + crescendo * time / duration) mouth_pressure += pressure_alpha * (desired_pressure - mouth_pressure) breath_state += breath_alpha * (noise[sample] - breath_state) played_pressure = mouth_pressure * (1.0 + turbulence * breath_state) vibrato = 0.0 if time > vibrato_start: vibrato = math.sin(2.0 * math.pi * vibrato_rate * (time - vibrato_start)) * ( 1.0 - math.exp(-(time - vibrato_start) / 0.3)) * vibrato_cents / 1200.0 lip_gap = rest_gap * (1.0 + 2.0 * vibrato) if tongue_rate > 0.0 and time > attack and time < duration - release: phase = (time - attack) * tongue_rate % 1.0 closure = math.exp(-0.5 * ((phase - 0.92) / 0.040)**2) lip_gap *= 1.0 - tongue_strength * closure 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): branch_history[branch] = (neck_history[branch] + neck_g[branch] / load_g[branch] * load_history[branch]) / ( 1.0 + neck_g[branch] / load_g[branch]) right_hand[branch_nodes[branch]] -= branch_history[branch] free_pressure = _solve_tridiagonal(right_hand, inverse_diagonal, edge_g, multipliers) position_history = ((4.0 * mass / timestep**2 + 2.0 * damping / timestep - stiffness) * reed_position + 4.0 * mass / timestep * reed_velocity + reed_force) / reed_denominator bore_pressure = free_pressure[0] source_flow = 0.0 new_position = reed_position new_velocity = reed_velocity for iteration in range(48): pressure_difference = played_pressure - bore_pressure new_position = position_history + reed_area * pressure_difference / reed_denominator for collision_step in range(3): penetration = max(0.0, new_position - lip_gap) contact = contact_stiffness * penetration**1.5 derivative = 1.5 * contact_stiffness * math.sqrt(penetration) residual = reed_denominator * (new_position - position_history) - reed_area * pressure_difference + contact new_position -= residual / (reed_denominator + derivative) penetration = max(0.0, new_position - lip_gap) contact_derivative = 1.5 * contact_stiffness * math.sqrt(penetration) position_derivative = -reed_area / (reed_denominator + contact_derivative) opening = max(0.0, lip_gap - new_position) new_velocity = 2.0 / timestep * (new_position - reed_position) - reed_velocity pressure_magnitude = abs(pressure_difference) root = math.sqrt(2.0 / AIR_DENSITY) * pressure_magnitude / math.sqrt(pressure_magnitude + 5.0) direction = 1.0 if pressure_difference >= 0.0 else -1.0 jet_factor = 1.0 + turbulence * noise[sample] jet_flow = direction * width * opening * root * jet_factor source_flow = jet_flow + reed_area * new_velocity flow_derivative = -width * opening * math.sqrt(2.0 / AIR_DENSITY) * ( 0.5 * pressure_magnitude + 5.0) / (pressure_magnitude + 5.0)**1.5 * jet_factor if opening > 0.0: flow_derivative -= direction * width * root * position_derivative * jet_factor flow_derivative += reed_area * 2.0 / timestep * position_derivative residual = bore_pressure - free_pressure[0] - port_response[0] * source_flow if abs(residual) < 1e-6: break derivative = 1.0 - port_response[0] * flow_derivative change = max(-2000.0, min(2000.0, residual / max(0.25, derivative))) bore_pressure -= change maximum_source_residual = max(maximum_source_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 reed_position = new_position reed_velocity = new_velocity if reed_position > lip_gap: contact_count += 1 for node in range(node_count): pressure[node] = free_pressure[node] + port_response[node] * source_flow peak_bore_pressure = max(peak_bore_pressure, abs(pressure[0])) accumulated_energy = 0.5 * mass * reed_velocity**2 + 0.5 * stiffness * reed_position**2 accumulated_energy += contact_stiffness * max(0.0, reed_position - lip_gap)**2.5 / 2.5 for node in range(node_count): capacitor_current = capacitor_g[node] * pressure[node] + capacitor_history[node] capacitor_history[node] = -capacitor_g[node] * pressure[node] - capacitor_current accumulated_energy += 0.5 * compliance[node] * pressure[node]**2 for edge in range(node_count - 1): pressure_difference = pressure[edge] - pressure[edge + 1] flow = edge_g[edge] * pressure_difference + edge_history[edge] edge_history[edge] = edge_g[edge] * ((2.0 * inertance[edge] / timestep - resistance[edge]) * flow + pressure_difference) accumulated_energy += 0.5 * inertance[edge] * flow**2 output_flow = 0.0 for branch in range(branch_count): flow = branch_g[branch] * pressure[branch_nodes[branch]] + branch_history[branch] load_pressure = (flow - load_history[branch]) / load_g[branch] neck_pressure = pressure[branch_nodes[branch]] - load_pressure neck_history[branch] = neck_g[branch] * ((2.0 * branch_inertance[branch] / timestep - branch_resistance[branch]) * flow + neck_pressure) accumulated_energy += 0.5 * branch_inertance[branch] * flow**2 if branch_open[branch]: load_flow = flow - load_pressure / load_resistance[branch] load_history[branch] = load_flow + timestep / (2.0 * load_inertance[branch]) * load_pressure accumulated_energy += 0.5 * load_inertance[branch] * load_flow**2 output_flow += flow else: load_history[branch] = -load_g[branch] * load_pressure - flow accumulated_energy += 0.5 * load_compliance[branch] * load_pressure**2 radiation = AIR_DENSITY / (4.0 * math.pi) * (output_flow - last_output_flow) / timestep last_output_flow = output_flow radiation_state += radiation_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] = accumulated_energy return audio, energy, peak_bore_pressure, contact_count, maximum_source_residual def render_note(written_midi, duration=3.0, player=None, geometry=None, reed=None, seed=20261011): low_octave = written_midi in (74, 75) default_pressure = 5500.0 if low_octave else 6500.0 if written_midi > 73 else 5000.0 default_gap = 0.00060 if low_octave else 0.00045 if written_midi > 73 else 0.00075 player = player or Player(pressure_pa=default_pressure) geometry = geometry or load_geometry() reed = reed or Reed(gap_m=default_gap) reed = replace(reed, gap_m=reed.gap_m * player.opening_ratio) if duration <= player.attack_seconds + player.release_seconds + player.tail_seconds or player.pressure_pa < 0.0: raise ValueError("A positive duration and nonnegative breath pressure are required") network = build_network(written_midi, geometry, reed) solver = prepare_solver(network) damping = 2.0 * reed.damping_ratio * math.sqrt(reed.mass_kg * reed.stiffness_n_m) reed_parameters = np.asarray((reed.mass_kg, reed.stiffness_n_m, damping, reed.width_m, reed.gap_m, reed.effective_area_m2, reed.contact_stiffness)) noise = np.random.default_rng(seed).normal(size=round(duration * MODEL_SAMPLE_RATE)) audio, energy, peak_pressure, contacts, source_residual = _simulate( duration - player.tail_seconds, player.pressure_pa, player.attack_seconds, player.release_seconds, player.vibrato_cents, player.vibrato_rate_hz, player.vibrato_start_seconds, player.turbulence, player.tongue_rate_hz, player.tongue_strength, player.crescendo, noise, reed_parameters, 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["branch_open"], *solver, ) if not np.all(np.isfinite(audio)) or not np.all(np.isfinite(energy)): raise FloatingPointError("The coupled acoustic model produced nonfinite output") if source_residual > 0.01: raise FloatingPointError(f"The reed/bore pressure solve failed to converge: {source_residual} Pa") return resample_poly(audio, 1, OVERSAMPLE), { "writtenMidi": written_midi, "writtenNote": note_name(written_midi), "concertMidi": written_midi - 9, "concertNote": note_name(written_midi - 9), "expectedFrequencyHz": midi_frequency(written_midi - 9), "breathPressurePa": player.pressure_pa, "closingPressurePa": reed.stiffness_n_m * reed.gap_m / reed.effective_area_m2, "effectiveReedGapMm": reed.gap_m * 1000.0, "peakBorePressurePa": float(peak_pressure), "reedContactSamples": int(contacts), "maximumCouplingResidualPa": float(source_residual), "peakStoredEnergyJ": float(np.max(energy)), "finalStoredEnergyJ": float(energy[-1]), "mouthpieceVolumeCm3": missing_cone_volume(geometry, reed) * 1e6, "acousticPressureNodes": len(network["compliance"]), "openHoles": [int(key) for key, is_open in zip(network["branch_keys"], network["branch_open"]) if is_open and key >= 0], "modelSampleRate": MODEL_SAMPLE_RATE, "energyTraceJ": energy.tolist(), }