i think ive got it

scripts/string_model_fdm.py

@@ -0,0 +1,112 @@
+
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+import sounddevice as sd
+
+sample_rate = 44100
+seconds = 5
+
+N = 100 # number of string segments
+I = int(sample_rate * seconds) # number of samples to simulate
+
+# string parameters
+rho = 8000 # density, steel, kg/m^3
+radius = 0.001 # meters
+S = math.pi*radius**2 # string cross sectional area, assuming circular
+mu = S*rho # linear mass density
+T = 1200 # string tension, N
+c = math.sqrt(T/mu) # transverse wave velocity
+kappa = 0.001 # stiffness coefficient
+sigma = 0.5 # damping coefficient
+L = 0.5 # length of string
+strike_position = 0.2 # x of impulse location
+impulse_width = 0.02 # x of impulse width
+impulse_velocity = 1000.0 # x/t of impulse magnitude
+sample_position = 0.1 # percentage along L of sampling for audio
+
+dx = L / N # delta x
+dt = 1/sample_rate
+#dt = 0.2 * dx / c
+if(dx**2 < (c*dt)**2 + 4*(kappa*dt/(dx**2))**2):
+    print("warning: possibly unstable due to not enough segments, increase N")
+
+# derived constants
+r1 = c * dt/dx 
+r2 = (c * dt/dx) ** 2
+s1 = kappa * dt/dx**2
+s2 = (kappa * dt/dx**2) ** 2
+
+# string grid
+x = np.linspace(0, L, N + 1) # linspace my beloved
+
+# state vectors
+y_last = np.zeros(N + 1)
+y_current = np.zeros(N + 1)
+y_next = np.zeros(N + 1)
+# output
+y_sample = np.zeros(I)
+
+# vector of initial velocity across the string as a result of the impulse
+v0 = impulse_velocity * np.exp(-((x - strike_position) ** 2)/(2 * impulse_width ** 2))
+
+# initial conditions:
+# y(x, 0) = 0 and dy/dt(x, 0) = v0(x)
+# boundary conditions: y(0, t) = y(L, t) = 0 as well as all dy/dts thereafter
+
+def show_plot():
+    plt.plot(np.arange(0, N+1, 1), y_current)
+    plt.grid()
+    plt.show()
+
+# first iteration
+for n in range(2, N-2):
+    y_current[n] = dt * v0[n] + 0.5 * r1**2 * (y_last[n-1] - 2*y_last[n] + y_last[n+1])
+    # y_xx = y_last[n+1] - 2*y_last[n] + y_last[n-1]
+    # y_xxxx = y_last[n-2] - 4*y_last[n-1] + 6*y_last[n] - 4*y_last[n+1] + y_last[n+2]
+
+    # term1 = (2 - 2*sigma*dt) * y_current[n]
+    # term2 = (2*sigma*dt - 1) * y_last[n]
+    # term3 = r2 * y_xx - s2 * y_xxxx
+    # y_current[n] = dt * v0[n] + 0.5 * r1**2 *(term1 + term2 + term3)
+
+n_sample = int(sample_position * L * N)
+y_sample[0] = y_last[n_sample]
+y_sample[1] = y_current[n_sample]
+
+# rest of the iterations
+for i in range(2, I):
+
+    for n in range(2, N-2):
+        # simple wave equation
+        # y_next[n] = 2*y_current[n] - y_last[n] + r**2 * (y_current[n-1] - 2*y_current[n] + y_current[n+1])
+        
+        # stiff wave equation
+        y_xx = y_current[n+1] - 2*y_current[n] + y_current[n-1]
+        y_xxxx = y_current[n-2] - 4*y_current[n-1] + 6*y_current[n] - 4*y_current[n+1] + y_current[n+2]
+
+        term1 = (2 - 2*sigma*dt) * y_current[n]
+        term2 = (2*sigma*dt - 1) * y_last[n]
+        term3 = r2 * y_xx - s2 * y_xxxx
+        y_next[n] = term1 + term2 + term3
+
+    y_next[0] = 0
+    y_next[1] = 0
+    y_next[N-1] = 0
+    y_next[N-2] = 0
+
+    # y_sample[i] = math.tanh(y_next[n_sample])
+    y_sample[i] = y_next[n_sample]
+
+    y_last = y_current.copy()
+    y_current = y_next.copy()
+
+    if(i % 1000 == 0):
+        print(i/I * 100, "% complete")
+
+plt.plot(np.arange(0, I, 1), y_sample)
+plt.grid()
+plt.show()
+
+sd.play(y_sample, sample_rate)
+sd.wait()

scripts/string_model_ss.py

@@ -30,22 +30,23 @@ import math
 # where A = diagonal matrix of A_n and B = [ 0 b_1 0 b_2 ... b_n]^T
 
 # lets start with 3 nodes where the string is tuned to 440hz (we'll get to arbitrary modes evantually)
-f_1 = 440 # fundamental frequency
+f_1 = 100 # fundamental frequency
 def f_n(n):
     return f_1 * n
-omega_1 = 2*math.pi*f_1
-b = 0.01
+omega_0 = 2*math.pi*f_1
+b = 1.5
+c = 2
 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 1.001*omega_1*math.sqrt(1+b*n**2)
+    return c*n*omega_0*math.sqrt(1+b*(n-1)**2)
 
 # x = [ q1, q1dot, q2, q2dot, q3, q3dot ]^T < --state vector
 omega_1 = omega_n(1)
 omega_2 = omega_n(2)
 omega_3 = omega_n(3)
-zeta_1 = 0.0001 # i guessed
-zeta_2 = 2 * zeta_1
-zeta_3 = 3 * zeta_1
+zeta_1 = 0.001 # i guessed
+zeta_2 = 1.5 * zeta_1
+zeta_3 = 2 * zeta_1
 A = [
         [          0,                 1,           0,                 0,           0,                 0],
         [-omega_1**2, -2*zeta_1*omega_1,           0,                 0,           0,                 0],
@@ -56,7 +57,7 @@ A = [
     ] # isnt this formatting gorgeous
 
 B = [ [0], [0.707], [0], [0], [0], [-0.707] ]
-c_1 = 0.0001
+c_1 = 0.002
 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]]
@@ -77,6 +78,11 @@ sys = sig.StateSpace(A, B, C, D)
 # step Response
 t, y = sig.impulse(sys, T=t)
 
+for i in range(t.size):
+    y[i] = math.tanh(y[i])/2
+    if(y[i] > 1):
+        print("what", y[i])
+
 # plotting
 plt.plot(t, y)
 plt.title("Step Response")