Modelando el sonido de las notas de guitarra usando el algoritmo Karplus-Strong en python

Conoce la nota de referencia A de la primera octava (440 Hz):





Suena doloroso, ¿no? ¿Qué más decir sobre el hecho de que la misma nota suena de manera diferente en diferentes instrumentos musicales? ¿Por que es esto entonces? Se trata de la presencia de armónicos adicionales que crean un timbre único para cada instrumento.



Pero nos interesa otra pregunta: ¿cómo simular este timbre único en una computadora?

Algoritmo estándar de Karplus-Strong

Ilustración tomada de este sitio .



La esencia del algoritmo es la siguiente:



1) Cree una matriz de tamaño N a partir de números aleatorios (N está directamente relacionado con la frecuencia fundamental del sonido).



2) Agregue al final de esta matriz el valor calculado por la siguiente fórmula:

$$

$$y(n)=\frac{y(n-N)+y(n-N-1)}{2},$$

Dónde $$$y$ es nuestra matriz.



3) Realizamos el punto 2 el número de veces requerido.



Comencemos a escribir el código:



1) Importe las bibliotecas requeridas.

import numpy as np
import scipy.io.wavfile as wave

2) Inicializamos las variables.

frequency = 82.41     #     
duration = 1          #    
sample_rate = 44100   #  

3) Crea ruido.

#  ,  frequency, ,        frequency .
#      sample_rate/length .
#  length = sample_rate/frequency.
noise = np.random.uniform(-1, 1, int(sample_rate/frequency))   

4) Cree una matriz para almacenar los valores y agregar ruido al principio.

samples = np.zeros(int(sample_rate*duration))
for i in range(len(noise)):
    samples[i] = noise[i]

5) Usamos la fórmula.

for i in range(len(noise), len(samples)):
    #   i   ,      .
    #  ,  i   ,       .
    samples[i] = (samples[i-len(noise)]+samples[i-len(noise)-1])/2

6) Normalizamos y traducimos al tipo de datos deseado.

samples = samples / np.max(np.abs(samples))  
samples = np.int16(samples * 32767)     

7) Guardar en archivo.

wave.write("SoundGuitarString.wav", 44100, samples)

8) Diseñemos todo en función. De hecho, ese es todo el código.

import numpy as np
import scipy.io.wavfile as wave
 
def GuitarString(frequency, duration=1., sample_rate=44100, toType=False):
    #  ,  frequency, ,        frequency .
    #      sample_rate/length .
    #  length = sample_rate/frequency.
    noise = np.random.uniform(-1, 1, int(sample_rate/frequency))      #  
 
    samples = np.zeros(int(sample_rate*duration))
    for i in range(len(noise)):
        samples[i] = noise[i]
    for i in range(len(noise), len(samples)):
        #   i   ,      .
        #  ,  i   ,       .
        samples[i] = (samples[i-len(noise)]+samples[i-len(noise)-1])/2
 
    if toType:
        samples = samples / np.max(np.abs(samples))  #   -1  1
        return np.int16(samples * 32767)             #     int16
    else:
        return samples
 
 
frequency = 82.41
sound = GuitarString(frequency, duration=4, toType=True)
wave.write("SoundGuitarString.wav", 44100, sound)

9) Corramos y obtengamos:





Para hacer que la cuerda suene mejor, mejoremos ligeramente la fórmula:

$$

$$y(n)=0.996\frac{y(n-N)+y(n-N-1)}{2}$$

Una sexta cuerda abierta (82,41 Hz) suena así:





La primera cuerda abierta (329,63 Hz) suena así:





Suena bien, ¿no?



Puede seleccionar infinitamente este coeficiente y encontrar el promedio entre un sonido hermoso y la duración, pero es mejor ir directamente al algoritmo avanzado Karplus-Strong.

Un poco sobre la transformación Z

Permitir $$$x$ es una matriz de valores de entrada, y$$$y$ es una matriz de valores de salida. Cada elemento de y se expresa mediante la siguiente fórmula:

$$

$$y(n)=x(n)+x(n-1).$$

Si el índice está fuera de la matriz, entonces el valor es 0. Es decir $$$x(0-1)=0$ . (Mire el código anterior, allí se usó implícitamente).



Esta fórmula se puede escribir en la transformada Z correspondiente:

$$

$$H(z)=1+z^{-1}.$$

Si la fórmula es así:

$$

$$y(n)=x(n)+x(n-1)-y(n-1).$$

Es decir, cada elemento de la matriz de entrada depende del elemento anterior de la misma matriz (excepto el elemento cero, por supuesto). Entonces la transformación Z correspondiente se ve así:

$$

$$H(z)=\frac{1+z^{-1}}{1+z^{-1}}.$$

El proceso inverso: de la transformada Z, obtenga la fórmula para cada elemento. Por ejemplo,

$$

$$H(z)=\frac{1+z^{-1}}{1-z^{-1}}.$$

$$

$$H(z)=\frac{Y(z)}{X(z)} =\frac{1+z^{-1}}{1-z^{-1}}.$$

$$

$$Y(z)*(1-z^{-1} )=X(z)*(1+z^{-1} ).$$

$$

$$Y(z)*1-Y(z)*z^{-1}= X(z)*1+X(z)*z^{-1}.$$

$$

$$y(n)-y(n-1)=x(n)+x(n-1).$$

$$

$$y(n)=x(n)+x(n-1)+y(n-1).$$

Si alguien no entiende, entonces la fórmula es: $$$Y(z)*α*z^{-k}=α*y(n-k)$dónde $$$α$- cualquier número real.



Si necesita multiplicar dos transformaciones Z entre sí, entonces$$$z^{-a}*z^{-b}=z^{-a-b}.$

Algoritmo extendido de Karplus-Strong

Ilustración tomada de este sitio.



Aquí hay una descripción rápida de cada función.

Parte I. Funciones que transforman el ruido inicial

1) Filtro de paso bajo de dirección de selección (filtro de paso bajo)$$$H_p (z)$...

$$

$$H_p (z)=\frac{1-p}{1-pz^{-1}},p ∈ [0,1).$$

Fórmula correspondiente:

$$

$$y(n)=(1-p)x(n)+py(n-1).$$

El código:

buffer = np.zeros_like(noise)
buffer[0] = (1 - p) * noise[0]
for i in range(1, N):
    buffer[i] = (1-p)*noise[i] + p*buffer[i-1]
noise = buffer

Siempre debe crear otra matriz para evitar errores. Quizás no podría haber sido usado aquí, pero en el siguiente filtro no puede prescindir de él.



2) Filtro de peine de posición de recogida (filtro de peine)$$$H_β (z)$...

$$

$$H_β (z)=1-z^{-int(β*N+1/2)},β∈(0,1).$$

Fórmula correspondiente:

$$

$$y(n)=x(n)-x(n-int(β*N+1/2)).$$

El código:

pick = int(beta*N+1/2)
if pick == 0:
    pick = N   #      
buffer = np.zeros_like(noise)
for i in range(N):
    if i-pick < 0:
        buffer[i] = noise[i]
    else:
        buffer[i] = noise[i]-noise[i-pick]
noise = buffer

En el primer párrafo de la página 13 de este documento, se escribe lo siguiente (no literalmente, pero con la preservación del significado): el coeficiente β imita la posición de la cuerda pulsada. Si$$$β=1/2$, entonces esto significa que el punteo se hizo en el medio de la cuerda. Si$$$β=1/10$ - el punteo se hizo en una décima parte de la cuerda desde el puente.

Parte II. Funciones relacionadas con la parte principal del algoritmo

Hay una trampa aquí que tenemos que solucionar. Por ejemplo, filtro amortiguador de cuerdas$$$H_d (z)$ escrito así: $$$H_d (z)=(1-S)+Sz^{-1}$... Pero la figura muestra que toma sentido de donde lo da. Es decir, resulta que las señales de entrada y salida de este filtro son una y la misma. Esto significa que cada filtro no se puede aplicar por separado, como en la sección anterior, todos los filtros deben aplicarse simultáneamente. Esto se puede hacer, por ejemplo, encontrando el producto de cada filtro. Pero este enfoque no es racional: al agregar o cambiar un filtro, tendrá que multiplicar todo nuevamente. Es posible hacer esto, pero no tiene sentido. Me gustaría cambiar el filtro con un solo clic y no multiplicar todo una y otra vez.

Dado que la señal de salida del filtro se considera la entrada para otro filtro, propongo escribir cada filtro como una función separada que llama a la función del filtro anterior dentro de sí mismo.

Creo que el código de ejemplo dejará claro lo que quiero decir.

1) Filtro de línea de retardo $$$z^{-N}.$

$$

$$H(z)=z^{-N}.$$

Fórmula correspondiente:

$$

$$y(n)=x(n-N).$$

El código:

#    ,     samples  0.
#    n-N<0   0,    .
def DelayLine(n):
    return samples[n-N]

2) Filtro amortiguador de cuerdas $$$H_d (z)$...

$$

$$H_d (z)=(1-S)+Sz^{-1},S∈[0,1].$$

En el algoritmo original $$$S=0.5.$

Fórmula correspondiente:

$$

$$y(n)=(1-S)x(n)+Sx(n-1).$$

El código:

# String-dampling filter.
# H(z) = 0.996*((1-S)+S*z^(-1)).    S = 0.5. S ∈ [0, 1]
# y(n)=0.996*((1-S)*x(n)+S*x(n-1))
def StringDampling_filter(n):
    return 0.996*((1-S)*DelayLine(n)+S*DelayLine(n-1))

En este caso, este filtro es el filtro de amortiguación One Zero String. Hay otras opciones, puedes leer sobre ellas aquí .



3) Filtro de paso todo de rigidez de cuerdas$$$H_s (z)$...

No importa cuánto busqué, por desgracia, no pude encontrar nada específico. Aquí este filtro está escrito en términos generales. Pero eso no funciona, ya que la parte más difícil es encontrar las probabilidades correctas. Hay algo más en este documento en la página 14, pero no tengo suficientes antecedentes matemáticos para comprender lo que está sucediendo allí y cómo usarlo. Si alguien puede, hágamelo saber.



4) Filtro de paso todo de ajuste de cuerdas de primer orden$$$H_ρ (z)$...

Página 6, parte inferior izquierda de este documento:

$$

$$H_ρ (z)=\frac{C+z^{-1}}{1+Cz^{-1}},C∈(-1,1).$$

Fórmula correspondiente:

$$

$$y(n)=Cx(n)+x(n-1)-Cy(n-1).$$

El código:

# First-order string-tuning allpass filter
# H(z) = (C+z^(-1))/(1+C*z^(-1)). C ∈ (-1, 1)
# y(n) = C*x(n)+x(n-1)-C*y(n-1)
def FirstOrder_stringTuning_allpass_filter(n):
    #        ,    ,  
    #    ,          samples.
    return C*(StringDampling_filter(n)-samples[n-1])+StringDampling_filter(n-1)

Debe recordarse que si agrega más filtros después de este filtro, tendrá que almacenar el valor pasado, porque ya no se almacenará en la matriz de muestras.

Dado que la longitud del ruido inicial es un número entero, descartamos la parte fraccionaria al contar. Esto provoca errores e inexactitudes. Por ejemplo, si la frecuencia de muestreo es 44100 y la longitud del ruido es 133 y 134, entonces las frecuencias de señal correspondientes son 331,57 Hz y 329,10 Hz. Y la frecuencia de las notas E de la primera octava (la primera cuerda al aire) es 329,63 Hz. Aquí la diferencia está en décimas, pero, por ejemplo, para el traste 15, la diferencia puede ser ya de varios Hz. Este filtro existe para reducir este error. No se puede utilizar si la frecuencia de muestreo es alta (muy alta: varios cientos de miles de Hz, o incluso más) o la frecuencia fundamental es baja, como, por ejemplo, para las cuerdas de bajo.

Hay otras variaciones, puedes leer sobre todas ellas allí .



5) Usamos nuestras funciones.

def Modeling(n):
    return FirstOrder_stringTuning_allpass_filter(n)
 
for i in range(N, len(samples)):
    samples[i] = Modeling(i)

Parte III. Filtro de paso bajo de nivel dinámico $$$H_L (z).$

$$$\check{ω}=ω\frac{T}{2}=2πf\frac{T}{2}=π\frac{f}{F_s}$dónde $$$f$ - frecuencia fundamental, $$$F_s$- frecuencia de muestreo.

Primero encontramos la matriz$$$y$ con la siguiente fórmula:

$$

$$H(z)=\frac{\check{ω}}{1+\check{ω}}\frac{1+z^{-1}}{1-\frac{1-\check{ω}}{1+\check{ω}} z^{-1}}$$

Fórmula correspondiente:

$$

$$y(n)=\frac{\check{ω}}{1+\check{ω}}(x(n)+x(n-1))+\frac{1-\check{ω}}{1+\check{ω}}y(n-1)$$

Luego aplicamos la siguiente fórmula:

$$

$$x(n)=L^{\frac{4}{3}}x(n)+(1-L)y(n),L∈(0,1)$$

El código:

# Dynamic-level lowpass filter. L ∈ (0, 1/3)
w_tilde = np.pi*frequency/sample_rate
buffer = np.zeros_like(samples)
buffer[0] = w_tilde/(1+w_tilde)*samples[0]
for i in range(1, len(samples)):
    buffer[i] = w_tilde/(1+w_tilde)*(samples[i]+samples[i-1])+(1-w_tilde)/(1+w_tilde)*buffer[i-1]
samples = (L**(4/3)*samples)+(1.0-L)*buffer

El parámetro L afecta el valor de disminución de volumen. Con sus valores iguales a 0.001, 0.01, 0.1, 0.32, el volumen de la señal disminuye en 60, 40, 20 y 10 dB, respectivamente.



Diseñemos todo en función. De hecho, ese es todo el código.

import numpy as np
import scipy.io.wavfile as wave
 
 
def GuitarString(frequency, duration=1., sample_rate=44100, p=0.9, beta=0.1, S=0.5, C=0.1, L=0.1, toType=False):
    N = int(sample_rate/frequency)            #    
 
    noise = np.random.uniform(-1, 1, N)   #  
 
    # Pick-direction lowpass filter (  ).
    # H(z) = (1-p)/(1-p*z^(-1)). p ∈ [0, 1)
    # y(n) = (1-p)*x(n)+p*y(n-1)
    buffer = np.zeros_like(noise)
    buffer[0] = (1 - p) * noise[0]
    for i in range(1, N):
        buffer[i] = (1-p)*noise[i] + p*buffer[i-1]
    noise = buffer
 
    # Pick-position comb filter ( ).
    # H(z) = 1-z^(-int(beta*N+1/2)). beta ∈ (0, 1)
    # y(n) = x(n)-x(n-int(beta*N+1/2))
    pick = int(beta*N+1/2)
    if pick == 0:
        pick = N   #      
    buffer = np.zeros_like(noise)
    for i in range(N):
        if i-pick < 0:
            buffer[i] = noise[i]
        else:
            buffer[i] = noise[i]-noise[i-pick]
    noise = buffer
 
    #    .
    samples = np.zeros(int(sample_rate*duration))
    for i in range(N):
        samples[i] = noise[i]
 
    #    ,     samples  0.
    #    n-N<0   0,    .
    def DelayLine(n):
        return samples[n-N]
 
    # String-dampling filter.
    # H(z) = 0.996*((1-S)+S*z^(-1)).    S = 0.5. S ∈ [0, 1]
    # y(n)=0.996*((1-S)*x(n)+S*x(n-1))
    def StringDampling_filter(n):
        return 0.996*((1-S)*DelayLine(n)+S*DelayLine(n-1))
 
    # First-order string-tuning allpass filter
    # H(z) = (C+z^(-1))/(1+C*z^(-1)). C ∈ (-1, 1)
    # y(n) = C*x(n)+x(n-1)-C*y(n-1)
    def FirstOrder_stringTuning_allpass_filter(n):
        #        ,    ,  
        #    ,          samples.
        return C*(StringDampling_filter(n)-samples[n-1])+StringDampling_filter(n-1)
 
    def Modeling(n):
        return FirstOrder_stringTuning_allpass_filter(n)
 
    for i in range(N, len(samples)):
        samples[i] = Modeling(i)
 
    # Dynamic-level lowpass filter. L ∈ (0, 1/3)
    w_tilde = np.pi*frequency/sample_rate
    buffer = np.zeros_like(samples)
    buffer[0] = w_tilde/(1+w_tilde)*samples[0]
    for i in range(1, len(samples)):
        buffer[i] = w_tilde/(1+w_tilde)*(samples[i]+samples[i-1])+(1-w_tilde)/(1+w_tilde)*buffer[i-1]
    samples = (L**(4/3)*samples)+(1.0-L)*buffer
 
    if toType:
        samples = samples/np.max(np.abs(samples))   #   -1  1
        return np.int16(samples*32767)              #     int16
    else:
        return samples
 
 
frequency = 82.51
sound = GuitarString(frequency, duration=4, toType=True)
wave.write("SoundGuitarString.wav", 44100, sound)

Una sexta cuerda abierta (82,41 Hz) suena así:





Y la primera cuerda abierta (329,63 Hz) suena así:





La primera cuerda no suena muy bien, por decirlo suavemente. Más parecido a una campana que a una cuerda. He intentado durante mucho tiempo comprender qué está mal con el algoritmo. Pensé que era un filtro sin usar. Después de días de experimentar, me di cuenta de que necesitaba aumentar la frecuencia de muestreo a al menos 100.000:





Suena mejor, ¿no?



En este documento se pueden leer complementos como tocar glissando o simular una cuerda simpática (págs. 11-12).



Aquí hay una pelea para ti:





Secuencia de acordes: CG # Am F. Strike: Six. La demora entre dos punteos consecutivos de la cuerda es de 0.015 segundos; el retraso entre dos golpes consecutivos en una batalla es de 0,205 segundos; el retraso en sí mismo en la batalla es de 0,41 segundos. El algoritmo ha cambiado el valor de L a 0,2.



Gracias por leer el articulo. ¡Buena suerte!