いろいろ調整してみたら自励振動させられてないっぽい。管体部分でループ(ハウリングみたいな状態)していただけだった。
物理モデリング音源の特許(だいたい権利失効してる)の実施例にある、励振部をpythonで実装した。
一応音は出るが、パラメータの調整が甘くか弱い音しか鳴らない。
とりあえず音が出る状態でメモしておく。
■pythonソース
import pyaudio
import numpy as np
from collections import deque
import matplotlib.pyplot as plt
import scipy.signal as signal
fs_sound = 44100
fs_internal = 44100
c = 340
class DelayLine:
def __init__(self, n_delay):
n_delay = np.round(n_delay).astype(np.int64)
self._dl = deque([0] * n_delay)
pass
def update(self, x_in):
self._dl.append(x_in)
x_out = self._dl.popleft()
return x_out
pass
class Filter(DelayLine):
def __init__(self, coef_fir):
n_taps = coef_fir.size
super().__init__(n_taps)
self.coef_fir = coef_fir[::-1]
pass
def update(self, x_in):
self._dl.append(x_in)
self._dl.popleft()
x_out = self._dl @ self.coef_fir
return x_out
pass
def sound_play(samples):
p = pyaudio.PyAudio()
stream = p.open(format=pyaudio.paFloat32,
channels=1,
rate=fs_sound,
frames_per_buffer=1024,
output=True)
t_end = 1
n_samples = t_end * fs_sound
n_t = np.arange(0, fs_sound)
f = 440
if np.max(np.abs(samples)) > 1:
samples = samples / np.max(np.abs(samples)) * 0.9
pass
stream.write(samples.astype(np.float32).tobytes())
stream.close()
pass
def filter_analy(a, fs):
freq, res = signal.freqz(a)
tate = 3
yoko = 1
plt.subplot(tate, yoko, 1)
plt.plot(0.5 * fs * freq / np.pi, np.abs(res))
plt.subplot(tate, yoko, 2)
plt.plot(0.5 * fs * freq / np.pi, np.angle(res))
plt.subplot(tate, yoko, 3)
plt.plot(a[::-1], ".-")
plt.show()
pass
def main():
fs = fs_internal
fc = 2000
freqs = np.array([0, fc, fc + 1000, fs / 2])
a = signal.firwin2(numtaps=500, freq=freqs / (fs / 2), gain=[1, 1, 0, 0], nfreqs=fs)
a = signal.minimum_phase(a, "homomorphic")
fl1 = Filter(a)
fl_press2reed = Filter(a)
fc = 30
freqs = np.array([0, fc, fc + 50, fs / 2])
a = signal.firwin2(numtaps=2001, freq=freqs / (fs / 2), gain=[0, 0, 1, 1], nfreqs=fs)
a = signal.minimum_phase(a, "homomorphic")
fl_outhpf = Filter(a)
dl_len = 0.8
n_dl = fs_internal * dl_len / c
print(n_dl)
dl1 = DelayLine(n_dl)
dl2 = DelayLine(n_dl)
dl_s = DelayLine(1)
t_simulation = 2
n_sim = fs_internal * t_simulation
xs_input = np.zeros(n_sim)
in_max = 3
in_min = 0.15
n_on1 = 10000
n_on2 = 20000
n_off1 = 30000
n_off2 = 80000
for n in range(n_sim):
if n_on1 <= n and n < n_on2:
x = in_min + (in_max - in_min) / (n_on2 - n_on1) * (n - n_on1)
pass
elif n_on2 <= n and n < n_off1:
x = in_max
pass
elif n_off1 <= n and n < n_off2:
x = in_min + (in_max - in_min) * (1 - 1 / (n_off2 - n_off1) * (n - n_off1))
pass
else:
x = in_min
pass
xs_input[n] = x
pass
xs_output = np.zeros(n_sim)
xs_reed = np.zeros(n_sim)
xs_volvel = np.zeros(n_sim)
k_refl_reed = 0.1
do1 = 0.0
do2 = 0.0
do_s = 0.0
def press2reed(x):
out = x
if x < -1:
out = -1
pass
elif x > 1:
out = 1
pass
else:
out = x
pass
return out
pass
def press2vel(x):
return x
pass
def volvel2press(x):
return x
pass
for n in range(n_sim):
x_n = xs_input[n]
s_press = do_s - x_n
s_press_reed = s_press - 0.8
s_reed_area = press2reed(fl_press2reed.update(s_press_reed))
xs_reed[n] = s_reed_area
s_velocity = press2vel(s_press)
s_volvel = s_reed_area * s_velocity
s_reedout = volvel2press(s_volvel)
xs_volvel[n] = s_volvel
s_forward_d = (1 + k_refl_reed) * s_reedout + (-k_refl_reed) * do2
do1 = dl1.update(s_forward_d)
x_n = fl_outhpf.update(do1)
s_backward_1 = fl1.update(x_n)
do2 = dl2.update(s_backward_1)
s_backward_2 = (1 - k_refl_reed) * do2 + k_refl_reed * s_reedout
do_s = dl_s.update(s_backward_2)
xs_output[n] = x_n
pass
plt.plot(xs_output, label="out")
plt.plot(xs_input, "--", label="in")
plt.plot(xs_reed, label="reed")
plt.plot(xs_volvel, label="volvel")
plt.legend()
plt.show()
sound_play(xs_output)
pass
if __name__ == '__main__':
main()