Modélisation du son des notes de guitare à l'aide de l'algorithme Karplus-Strong en python

Rencontrez la note de référence A de la première octave (440 Hz):





Cela semble douloureux, n'est-ce pas? Que dire de plus sur le fait que la même note sonne différemment sur différents instruments de musique. Pourquoi en est-il ainsi? Tout dépend de la présence d' harmoniques supplémentaires qui créent un timbre unique pour chaque instrument.



Mais nous sommes intéressés par une autre question: comment simuler ce timbre unique sur un ordinateur?

Algorithme standard de Karplus-Strong

Illustration tirée de ce site .



L'essence de l'algorithme est la suivante:



1) Créer un tableau de taille N à partir de nombres aléatoires (N est directement lié à la fréquence fondamentale du son).



2) Ajoutez à la fin de ce tableau la valeur calculée par la formule suivante:

$$

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

Où $$$y$Est notre tableau.



3) Nous effectuons le point 2 le nombre de fois requis.



Commençons à écrire le code:



1) Importez les bibliothèques requises.

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

2) Nous initialisons les variables.

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

3) Créez du bruit.

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

4) Créez un tableau pour stocker les valeurs et ajoutez du bruit au début.

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

5) Nous utilisons la formule.

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

6) Nous normalisons et traduisons dans le type de données souhaité.

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

7) Enregistrer dans un fichier.

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

8) Concevons tout comme une fonction. En fait, c'est tout le code.

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) Courons et obtenons:





Pour améliorer le son de la corde, améliorons légèrement la formule:

$$

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

Une sixième corde ouverte (82,41 Hz) ressemble à ceci:





La première chaîne ouverte (329,63 Hz) ressemble à ceci:





Ça sonne bien, n'est-ce pas?



Vous pouvez sélectionner ce coefficient à l'infini et trouver la moyenne entre le beau son et la durée, mais il vaut mieux aller directement à l'algorithme Advanced Karplus-Strong.

Un peu de Z-transform

Laisser être $$$x$ Est un tableau de valeurs d'entrée, et $$$y$- un tableau de valeurs de sortie. Chaque élément du tableau y est exprimé par la formule suivante:

$$

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

Si l'index est en dehors du tableau, la valeur est 0. C'est-à-dire $$$x(0-1)=0$... (Regardez le code précédent, il était implicitement utilisé).



Cette formule peut être écrite dans la transformée Z correspondante:

$$

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

Si la formule est comme ça:

$$

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

Autrement dit, chaque élément du tableau d'entrée dépend de l'élément précédent du même tableau (à l'exception de l'élément zéro, bien sûr). Ensuite, la transformation Z correspondante ressemble à ceci:

$$

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

Processus inverse: obtenez la formule pour chaque élément à partir de la transformation en Z. Par exemple,

$$

$$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 quelqu'un ne comprend pas, la formule est: $$$Y(z)*α*z^{-k}=α*y(n-k)$où $$$α$- n'importe quel nombre réel.



Si vous devez multiplier deux transformations Z l'une par l'autre, alors$$$z^{-a}*z^{-b}=z^{-a-b}.$

Algorithme Karplus-Strong étendu

Illustration tirée de ce site.



Voici une brève description de chaque fonctionnalité.

Partie I. Fonctions qui transforment le bruit initial

1) Filtre passe-bas dans le sens du prélèvement (filtre passe-bas)$$$H_p (z)$...

$$

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

Formule correspondante:

$$

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

Le code:

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

Vous devez toujours créer un autre tableau pour éviter les erreurs. Peut-être qu'il n'aurait pas pu être utilisé ici, mais dans le prochain filtre, vous ne pouvez pas vous en passer.



2) Filtre en peigne de sélection (filtre en peigne)$$$H_β (z)$...

$$

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

Formule correspondante:

$$

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

Le code:

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

Dans le premier paragraphe de la page 13 de ce document, ce qui suit est écrit (pas littéralement, mais avec la préservation du sens): le coefficient β imite la position de la corde pincée. Si$$$β=1/2$, cela signifie que le pincement a été fait au milieu de la corde. Si$$$β=1/10$ - le pincement a été fait sur un dixième de la corde du pont.

Deuxieme PARTIE. Fonctions liées à la partie principale de l'algorithme

Il y a ici un piège que nous devons contourner. Par exemple, le filtre String-dampling$$$H_d (z)$ écrit comme ceci: $$$H_d (z)=(1-S)+Sz^{-1}$... Mais l'image montre qu'il prend le sens là où il le donne. Autrement dit, il s'avère que les signaux d'entrée et de sortie pour ce filtre sont identiques. Cela signifie que chaque filtre ne peut pas être appliqué séparément, comme dans la section précédente, tous les filtres doivent être appliqués simultanément. Cela peut être fait, par exemple, en recherchant le produit de chaque filtre. Mais cette approche n'est pas rationnelle: lors de l'ajout ou du changement d'un filtre, vous devrez tout multiplier à nouveau. Il est possible de le faire, mais cela n'a aucun sens. Je voudrais changer le filtre en un clic, et ne pas tout multiplier encore et encore.

Puisque le signal de sortie du filtre est considéré comme l'entrée d'un autre filtre, je propose d'écrire chaque filtre comme une fonction distincte qui appelle la fonction du filtre précédent à l'intérieur de lui-même.

Je pense que l'exemple de code expliquera clairement ce que je veux dire.

1) Filtre de ligne de retard $$$z^{-N}.$

$$

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

Formule correspondante:

$$

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

Le code:

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

2) Filtre dampling de cordes $$$H_d (z)$...

$$

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

Dans l'algorithme d'origine $$$S=0.5.$

Formule correspondante:

$$

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

Le code:

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

Dans ce cas, ce filtre est le filtre One Zero String-dampling. Il existe d'autres options, vous pouvez les lire ici .



3) Filtre passe-tout à rigidité des cordes$$$H_s (z)$...

Peu importe à quel point je regardais, hélas, je n'ai rien trouvé de spécifique. Ici, le filtre est écrit en termes généraux. Mais cela ne fonctionne pas car le plus difficile est de trouver les bonnes chances. Il y a autre chose dans ce document à la page 14, mais je n'ai pas assez de base mathématique pour comprendre ce qui s'y passe et comment l'utiliser. Si quelqu'un peut, faites le moi savoir.



4) Filtre passe-tout d'accordage des cordes du premier ordre$$$H_ρ (z)$...

Page 6, en bas à gauche de ce document:

$$

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

Formule correspondante:

$$

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

Le code:

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

Il ne faut pas oublier que si vous ajoutez plus de filtres après ce filtre, vous devrez stocker la valeur passée, car elle ne sera plus stockée dans le tableau des échantillons.

Puisque la longueur du bruit initial est un entier, nous jetons la partie fractionnaire lors du comptage. Cela provoque des erreurs et des inexactitudes. Par exemple, si la fréquence d'échantillonnage est de 44100 et la longueur du bruit est de 133 et 134, les fréquences de signal correspondantes sont 331,57 Hz et 329,10 Hz. Et la fréquence des notes E de la première octave (la première corde ouverte) est de 329,63 Hz. Ici, la différence est en dixièmes, mais, par exemple, pour la 15e frette, la différence peut déjà être de plusieurs Hz. Pour réduire cette erreur, ce filtre existe. Il peut être omis si la fréquence d'échantillonnage est élevée (vraiment élevée: plusieurs centaines de milliers de Hz, voire plus) ou si la fréquence fondamentale est basse, comme par exemple pour les cordes de basse.

Il existe d'autres variantes, vous pouvez toutes les lire ici .



5) Nous utilisons nos fonctions.

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

Partie III. Filtre passe-bas de niveau dynamique $$$H_L (z).$

$$$\check{ω}=ω\frac{T}{2}=2πf\frac{T}{2}=π\frac{f}{F_s}$où $$$f$ - la fréquence fondamentale, $$$F_s$- fréquence d'échantillonnage.

Nous trouvons d'abord le tableau$$$y$ avec la formule suivante:

$$

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

Formule correspondante:

$$

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

Ensuite, nous appliquons la formule suivante:

$$

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

Le code:

# 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

Le paramètre L affecte la valeur de diminution du volume. Avec ses valeurs égales à 0,001, 0,01, 0,1, 0,32, le volume du signal diminue respectivement de 60, 40, 20 et 10 dB.



Concevons tout comme une fonction. En fait, c'est tout le code.

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)

Une sixième corde ouverte (82,41 Hz) ressemble à ceci:





Et la première chaîne ouverte (329,63 Hz) ressemble à ceci:





La première corde ne sonne pas très bien, pour le moins dire. Plus comme une cloche qu'une corde. J'essaie depuis très longtemps de comprendre ce qui ne va pas avec l'algorithme. Je pensais que c'était un filtre inutilisé. Après des jours d'expérimentation, j'ai réalisé que je devais augmenter la fréquence d'échantillonnage à au moins 100 000:





Ça sonne mieux, n'est-ce pas?



Des modules complémentaires tels que la lecture de glissando ou la simulation d'une chaîne sympathique peuvent être lus dans ce document (p. 11-12).



Voici un combat:





Séquence d'accords: CG # Am F. Strike: Six. Le délai entre deux pincements consécutifs de la corde est de 0,015 seconde; le délai entre deux coups sûrs consécutifs dans une bataille est de 0,205 seconde; le délai lui-même au combat est de 0,41 seconde. L'algorithme a changé la valeur de L à 0,2.



Merci d'avoir lu l'article. Bonne chance!