sync

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()