From a2df3452167af167d0adb8a5f2d87cb45b5fa9d7 Mon Sep 17 00:00:00 2001 From: Bliblank Date: Fri, 19 Jun 2026 13:50:49 -0500 Subject: [PATCH] sync --- scripts/string_model_de.py | 53 ++++++++++++++++++++++---------------- scripts/string_model_ss.py | 11 ++++++-- 2 files changed, 40 insertions(+), 24 deletions(-) diff --git a/scripts/string_model_de.py b/scripts/string_model_de.py index 2f5e95e..ecbfd5e 100644 --- a/scripts/string_model_de.py +++ b/scripts/string_model_de.py @@ -2,6 +2,7 @@ import matplotlib.pyplot as plt import numpy as np import math +import sounddevice as sd # https://people.cs.uchicago.edu/~ridg/stabil/pianostring.pdf # a differential equations model for a piano string @@ -104,23 +105,26 @@ import math # simulation parameters f1 = 50 # fundamental frequency f_e = 44100 # sampling frequency -N = 50 # number of string segments +N = 200 # number of string segments delta_t = 1/f_e # time step L = 0.5 # string length, meters delta_x = L/N # spatial step -H = int(f_e * 0.1) # length of simulation in time +H = int(f_e * 2) # length of simulation in time # string parameters -E = 200 * 10**9 # youngs modulus, steel = 200GPa +E = 100 * 10**9 # youngs modulus, steel = 200GPa rho = 8000 # density, steel, kg/m^3 -radius = 0.002 # meters +radius = 0.005 # meters kappa = radius/2 # radius of gyration, r/2 for a circular string S = math.pi*radius**2 # string cross sectional area, assuming circular mu = S*rho # linear mass density M_S = mu*L # string mass -T = 10000 # string tension, N +T = 500 # string tension, N c = math.sqrt(T/mu) # transverse wave velocity -stiffness = kappa**2 * delta_t**2 / delta_x**4 # string stiffness parameter +I = math.pi*radius**4 / 4 +k = math.sqrt(E*I/(rho*S)) +# stiffness = k**2 * delta_t**2 / delta_x**4 # string stiffness parameter +stiffness = 0 sigma = 1 # decay rate tau = 1/sigma # decay time omega = 1/f1 # angular frequency @@ -128,14 +132,14 @@ omega = 1/f1 # angular frequency # hammer parameters M_H = 0.5 # hammer mass, kg HSMR = M_H/M_S # hammer-mass string ratio -V_H_0 = 10 # initial hammer velocity at t=0, m/s +V_H_0 = 10000000000 # initial hammer velocity at t=0, m/s x_0 = L/2 # distance of hammer from agraffe alpha = x_0 / L # relative hammer striking position # empirical constants -b_1 = 1 # some constant -b_3 = 0.001 # some constant -K = 10 # hammer stiffness +b_1 = 0 # some constant +b_3 = 0 # some constant +K = 1 # hammer stiffness p = 1 # stiffness nonlinear exponent # derived components @@ -152,7 +156,7 @@ x = [0] * N # current string position x_initial = [0] * N # assuming string at rest last_x1 = [0] * N # string position from last timestep last_x2 = [0] * N # string position from two timesteps ago -x_sample = int(alpha*N) # location where we sample the string position for signal +x_sample = int(0.8*N) # location where we sample the string position for signal # hammer F_H_current = 0 # current force exterted by hammer @@ -175,7 +179,7 @@ def plot_current(): def spatial_window(i): if(i < x_hammer+i_H/2 and i > x_hammer-i_H/2): - return 1 + return abs(i-x_hammer-i_H/2) else: return 0 @@ -196,8 +200,6 @@ eta_current = eta_next for i in range(1, N-1): x[i] = last_x1[i+1] + last_x1[i-1] - x_initial[i] + ((delta_t**2)*N*F_H_current*spatial_window(i))/M_S -plot_current() - last_x2 = last_x1.copy() last_x1 = x.copy() F_H_current = K * abs(eta_current - x[x_hammer])**p @@ -217,25 +219,29 @@ for n in range(H): # time F_H_current = K * abs(push)**p else: F_H_current = 0 + F_H_current = 0 next_x = [0] * N for i in range(2, N-2): # space - term_1 = a_1*x[i] + a_2*last_x1[i] - term_2 = a_3*(x[i+1] + x[i-1]) + a_4*(x[i+2] + x[i-2]) - term_3 = a_5*(last_x1[i+1] + last_x1[i-1] + last_x2[i]) - term_4 = (delta_t**2 * N*F_H_current * spatial_window(i))/M_S + # term_1 = a_1*x[i] + a_2*last_x1[i] + # term_2 = a_3*(x[i+1] + x[i-1]) + a_4*(x[i+2] + x[i-2]) + # term_3 = a_5*(last_x1[i+1] + last_x1[i-1] + last_x2[i]) + # term_4 = (delta_t**2 * N*F_H_current * spatial_window(i))/M_S - next_x[i] = term_1 + term_2 + term_3 + term_4 + # next_x[i] = term_1 + term_2 + term_3 + term_4 + next_x[i] = 2*(1-r**2)*x[i] + r**2 * (x[i+1] + x[i-1]) - last_x1[i] + + x[0] = 0 + x[1] = 0 + x[N-2] = 0 + x[N-1] = 0 last_x2 = last_x1.copy() last_x1 = x.copy() x = next_x.copy() - if(n == 30): - plot_current() - x_out[n] = x[x_sample] # plotting @@ -245,3 +251,6 @@ plt.xlabel("t") plt.ylabel("x(t, x_sample)") plt.grid() plt.show() + +sd.play(x_out, f_e) +sd.wait() diff --git a/scripts/string_model_ss.py b/scripts/string_model_ss.py index 936fb14..a736961 100644 --- a/scripts/string_model_ss.py +++ b/scripts/string_model_ss.py @@ -2,6 +2,7 @@ import scipy.signal as sig import matplotlib.pyplot as plt import numpy as np +import sounddevice as sd import math # simple first order step response simulation @@ -32,8 +33,11 @@ import math f_1 = 440 # fundamental frequency def f_n(n): return f_1 * n +omega_1 = 2*math.pi*f_1 +b = 0.01 def omega_n(n): - return 2*math.pi*f_n(n) # a cooler option would be omega_n = c*n*omega_1*sqrt(1+B*n^2) to factor in string stiffness to its vibration mode + # return 2*math.pi*f_n(n) # a cooler option would be omega_n = c*n*omega_1*sqrt(1+B*n^2) to factor in string stiffness to its vibration mode + return 1.001*omega_1*math.sqrt(1+b*n**2) # x = [ q1, q1dot, q2, q2dot, q3, q3dot ]^T < --state vector omega_1 = omega_n(1) @@ -52,7 +56,7 @@ A = [ ] # isnt this formatting gorgeous B = [ [0], [0.707], [0], [0], [0], [-0.707] ] -c_1 = 0.001 +c_1 = 0.0001 c_2 = c_1 / 2 c_3 = c_1 / 3 C = [[-c_1*omega_1**2, -2*c_1*zeta_1*omega_1, -c_2*omega_2**2, -2*c_2*zeta_2*omega_2, -c_3*omega_3**2, -2*c_3*zeta_3*omega_3]] @@ -80,3 +84,6 @@ plt.xlabel("t") plt.ylabel("y") plt.grid() plt.show() + +sd.play(y, 44100) +sd.wait()