trumpet physical model (???)

trumpet.py

#!/usr/bin/env python3

import numpy as np
import math
import sys

speed_of_sound = 34300 # centimeters/second
rate = 24000 # samples/second

# seconds to render
render_time = 2

# f3
fundamental = 174.61 # pitch in hertz tuning of instrument lowest note

output_gain = 0.5
#noise = 0.025
noise = 0.01
vib_rate = 2
vib_depth = 0.5
ramp_time = .5
max_pressure = 0.9
note_time = 1.5

# time in seconds of a single sample
dt = 1 / rate

# the physical wavelength size of a single sample
sample_wavelength = speed_of_sound / rate

# wavelength of the fundamental tuning in cm
fundamental_wavelength = speed_of_sound / fundamental

# length of bore will be 1/2th of that
# this gonna be a even bore
bore_length = fundamental_wavelength / 2

# number of samples across the bore length
bore_sample_length = bore_length / sample_wavelength

# separate out the int part from the fractional part
delay_samples = bore_sample_length * 2
num_segs, fractional_delay = divmod(delay_samples, 1)
num_segs = int(num_segs)

# num samples to render
render_samples = render_time * rate

# create the waveguide for the round trip
shape = (int(num_segs),)
delay = np.zeros(shape=shape)

print(f'rate               = {rate}hz')
print(f'sample wavelen     = {sample_wavelength}cm')
print(f'bore length cm     = {bore_length}cm')
print(f'bore sample length = {bore_sample_length}cm')
print(f'number wg segs   = {num_segs}')
print(f'fractional delay = {fractional_delay}')

def run_delay(delay, input_sample):
    n = len(delay)
    output_sample = delay[n - 1]

    # propogate the delay line
    for i in range(n - 1):
        delay[n - 1 - i] = delay[n - 2 - i]

    delay[0] = input_sample

    # TODO: implement all pass for fractional delay
    return output_sample

x = 0
v = 0
k = (fundamental * 8 * math.pi)**2
ne = 10
nk = 2
b = 0.1
nb = 5

input_pressure_influence = 1
coupling                 = 0.5
reed_gain = 1
reflection = 0.5

delay_output = 0
def run(i):
    global x, v
    global delay_output
    t = i * dt
    j = min(1, t / ramp_time) if t < note_time \
    else max(0, 1 - (t - note_time) / ramp_time)
    input_pressure = j * max_pressure
    input_pressure *= 1 + math.sin(i * dt * 2 * math.pi * vib_rate) * vib_depth
    input_pressure += np.random.rand() * noise * input_pressure

    feedback = reflection * delay_output
    reed_input = delay_output - feedback
    reed_output = input_pressure * x**2 * reed_gain
    delay_input = feedback + reed_output
    delay_output = run_delay(delay, delay_input)
    output = delay_input * output_gain

    a = -k * (x + nk * x ** (2 * ne + 1)) - b * (v + nb * x ** 2) * fundamental
    a += input_pressure * input_pressure_influence * k - coupling * reed_input * k
    v += a * dt
    x += v * dt
    x = min(1, max(-1, x))

    return output, input_pressure, x

# render it now
print(f'bout to render {int(render_samples)} samples...')
buf1 = np.empty(shape=(int(render_samples),))
buf2 = np.empty(shape=(int(render_samples),))
buf3 = np.empty(shape=(int(render_samples),))
samps = range(int(render_samples))
for i in samps:
    buf1[i], buf2[i], buf3[i] = run(i)

print ("DONE")

# plot it
import matplotlib.pyplot as plt
plt.plot(samps, buf3, '-', samps, buf1, '--', samps, buf2, '--')
plt.legend(['reed x', 'bore out', 'pressure in'])
plt.figtext(.6, .5, f'bore length = {bore_length:.2f}cm')
plt.figtext(.6, .4, f'speed of sound = {speed_of_sound}cm/s')
plt.ylabel('sound pressure level')
plt.xlabel('time')
plt.xticks([0, render_samples/2, render_samples], ['0 ms', f'{render_time*500} ms', f'{render_time*1000} ms'])
plt.show()

# output it
import wave
out_buf = buf1
with wave.open('out.wav', 'wb') as f:
    f.setnchannels(1) # mono
    f.setsampwidth(1) # 8 bit
    f.setframerate(rate)
    data = [int((x/2+1)*128) for x in out_buf]
    f.writeframesraw( bytes(data) )