Modeling the sound of guitar notes using the Karplus-Strong algorithm in python

Meet the reference note A of the first octave (440 Hz):





Sounds painful, doesn't it? What else to say about the fact that the same note sounds differently on different musical instruments. Why is it so? It's all about the presence of additional harmonics that create a unique timbre for each instrument.



But we are interested in another question: how to simulate this unique timbre on a computer?

Standard Karplus-Strong algorithm

Illustration taken from this site .



The essence of the algorithm is as follows:



1) Create an array of size N from random numbers (N is directly related to the fundamental sound frequency).



2) Add to the end of this array the value calculated by the following formula:

$$

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

Where $$$y$Is our array.



3) We carry out point 2 the required number of times.



Let's start writing the code:



1) Import the required libraries.

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

2) We initialize the variables.

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

3) Create noise.

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

4) Create an array to store the values ​​and add noise at the beginning.

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

5) We use the formula.

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

6) We normalize and translate into the desired data type.

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

7) Save to file.

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

8) Let's design everything as a function. Actually, that's all the 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) Let's run and get:





To make the string sound better, let's slightly improve the formula:

$$

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

An open sixth string (82.41 Hz) sounds like this:





The open first string (329.63 Hz) sounds like this:





Sounds good, doesn't it?



You can endlessly select this coefficient and find the average between beautiful sound and duration, but it is better to go straight to the Advanced Karplus-Strong algorithm.

A little about Z-transform

Let be $$$x$ Is an array of input values, and $$$y$- an array of output values. Each element in y is expressed by the following formula:

$$

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

If the index is outside the array, then the value is 0. That is $$$x(0-1)=0$... (Look at the previous code, there it was implicitly used).



This formula can be written in the corresponding Z-transform:

$$

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

If the formula is like this:

$$

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

That is, each element of the input array depends on the previous element of the same array (except for the zero element, of course). Then the corresponding Z-transform looks like this:

$$

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

Reverse process: get the formula for each element from the Z-transform. For instance,

$$

$$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).$$

If someone did not understand, then the formula is: $$$Y(z)*α*z^{-k}=α*y(n-k)$where $$$α$- any real number.



If you need to multiply two Z-transformations by each other, then$$$z^{-a}*z^{-b}=z^{-a-b}.$

Extended Karplus-Strong algorithm

Illustration taken from this site.



Here's a quick summary of each feature.

Part I. Functions that transform the initial noise

1) Pick-direction lowpass filter (low pass filter)$$$H_p (z)$...

$$

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

Corresponding formula:

$$

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

The 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

You should always create another array to avoid errors. Maybe it could not have been used here, but in the next filter you cannot do without it.



2) Pick-position comb filter (comb filter)$$$H_β (z)$...

$$

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

Corresponding formula:

$$

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

The 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

In the first paragraph on page 13 of this document, the following is written (not literally, but with the preservation of the meaning): the coefficient β mimics the position of the plucked string. If$$$β=1/2$, it means that the pluck was made in the middle of the string. If$$$β=1/10$ - the pluck was made on one tenth of the string from the bridge.

Part II. Functions related to the main part of the algorithm

There is a trap here that we have to work around. For example, String-dampling filter$$$H_d (z)$ written like this: $$$H_d (z)=(1-S)+Sz^{-1}$... But the picture shows that he takes meaning from where he gives it. That is, it turns out that the input and output signals for this filter are one and the same. This means that each filter cannot be applied separately, as in the previous section, all filters must be applied simultaneously. This can be done, for example, by finding the product of each filter. But this approach is not rational: when adding or changing a filter, you will have to multiply everything again. It is possible to do this, but it makes no sense. I would like to change the filter in one click, and not multiply everything over and over again.

Since the output signal from the filter is considered the input for another filter, I propose to write each filter as a separate function that calls the function of the previous filter inside itself.

I think the example code will make it clear what I mean.

1) Delay Line filter $$$z^{-N}.$

$$

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

Corresponding formula:

$$

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

The code:

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

2) String-dampling filter $$$H_d (z)$...

$$

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

In the original algorithm $$$S=0.5.$

Corresponding formula:

$$

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

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

In this case, this filter is the One Zero String-dampling filter. There are other options, you can read about them here .



3) String-stiffness allpass filter$$$H_s (z)$...

No matter how much I looked, alas, I could not find anything specific. Here this filter is written in general terms. But that doesn't work as the hardest part is finding the right odds. There is something else in this document on page 14, but I do not have enough mathematical background to understand what is happening there and how to use it. If anyone can, let me know.



4) First-order string-tuning allpass filter$$$H_ρ (z)$...

Page 6, bottom left in this document:

$$

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

Corresponding formula:

$$

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

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

It should be remembered that if you add more filters after this filter, you will have to store the past value, because it will no longer be stored in the samples array.

Since the length of the initial noise is an integer, we throw out the fractional part when counting. This causes errors and inaccuracies. For example, if the sampling rate is 44100 and the noise length is 133 and 134, then the corresponding signal frequencies are 331.57 Hz and 329.10 Hz. And the frequency of the E notes of the first octave (the first open string) is 329.63 Hz. Here the difference is in tenths, but, for example, for the 15th fret, the difference may already be several Hz. This filter exists to reduce this error. It can be omitted if the sampling frequency is high (really high: several hundred thousand Hz, or even more) or the fundamental frequency is low, as, for example, for bass strings.

There are other variations, you can read about them all there .



5) We use our functions.

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

Part III. Dynamic Level Lowpass Filter $$$H_L (z).$

$$$\check{ω}=ω\frac{T}{2}=2πf\frac{T}{2}=π\frac{f}{F_s}$where $$$f$ - fundamental frequency, $$$F_s$- sampling frequency.

First we find the array$$$y$ with the following formula:

$$

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

Corresponding formula:

$$

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

Then we apply the following formula:

$$

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

The 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

The L parameter affects the volume decrease value. With its values ​​equal to 0.001, 0.01, 0.1, 0.32, the signal volume decreases by 60, 40, 20 and 10 dB, respectively.



Let's design everything as a function. Actually, that's all the 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)

An open sixth string (82.41 Hz) sounds like this:





And the open first string (329.63 Hz) sounds like this:





The first string does not sound very good, to put it mildly. More like a bell than a string. I have been trying for a very long time to understand what is wrong with the algorithm. Thought it was an unused filter. After days of experimenting, I realized that I needed to increase the sample rate to at least 100,000:





Sounds better, doesn't it?



Add-ons such as playing glissando or simulating a sympathetic string can be read in this document (pp. 11-12).



Here's a fight for you:





Chord Sequence: CG # Am F. Strike: Six. The delay between two consecutive plucking of the string is 0.015 seconds; the delay between two consecutive hits in a battle is 0.205 seconds; the delay itself in battle is 0.41 seconds. The algorithm has changed the value of L to 0.2.



Thanks for reading the article. Good luck!