# Karplus Strong pluck generation

I want to use a simple implementation of the Karplus-Strong algorithm to generate pluck sounds and chords.

My existing code seems correct but is slow when more than a few plucks are being generated. Is there a way to improve the performance of this. In particular is there a better way to use numpy here, or would something like a collections.deque be faster.

import numpy as np

sample_rate = 44100
damping = 0.999

def generate(f, vol, nsamples):
"""Generate a Karplus-Strong pluck.

Arguments:
f -- the frequency, as an integer
vol -- volume, a float in [0.0, 1.0]
nsamples -- the number of samples to generate. To generate a pluck t
seconds long, one needs t * sample_rate samples. nsamples is an int.

Return value:
A numpy array of floats with length nsamples.
"""

N = sample_rate // f
buf = np.random.rand(N) * 2 - 1
samples = np.empty(nsamples, dtype=float)

for i in range(nsamples):
samples[i] = buf[i % N]
avg = damping * 0.5 * (buf[i % N] + buf[(1 + i) % N])
buf[i % N] = avg

return samples * vol

def generate_chord(f, nsamples):
"""Generate a chord of several plucks

Arguments:
f -- a (usually short) list of integer frequencies
nsamples -- the length of the chord, a chord of length t seconds needs
t * sample_rate, an integer.

Return value:
A numpy array
"""
samples = np.zeros(nsamples, dtype=float)
for note in f:
samples += generate(note, 0.5, nsamples)
return samples

if __name__ == "__main__":
import matplotlib.pyplot as plt
from scipy.io.wavfile import write

strum = generate_chord(  # A Major
[220, 220 * 5 // 4, 220 * 6 // 4, 220 * 2], sample_rate * 5)
plt.plot(strum)
plt.show()

write("pluck.wav", sample_rate, strum)


This produces a sound file in pluck.wav with a waveform like this: • Could you add an example that is too slow? Right now, this code runs in less than a second on my machine... – Graipher Feb 16 '17 at 15:37

I was trying to improve your generate function and stumbled on a potential bug. You continuously update the buffer by adding two following values. This should be vectorizable easily by acting on the whole buffer at the same time, by doing something like this:

for i in range(0, nsamples - N, N):
samples[i: i + N] = buf[:]
buf = damping * 0.5 * (buf + (np.roll(buf, -1)))
# fractional buffer
i += N
k = nsamples - i
if k:
samples[i:] = buf[:k]


This differs from your algorithm for every time step (i.e. one pass on the buf) only in the last element. This is because this code acts on the whole array at $t_{i}$ to generate the buffer at $t_{i+1}$, whereas your code uses the value of buf already from $t_{i+1}$.

If you can live with this difference (I can not check if it sounds any different, right now), the code execution time goes down from about 0.71 seconds to 0.20 seconds.

# Congratulations

You mostly follows standards which is good. You also have good documentation which is great. All in all the code reads well and is easily understandable.

Style-wise, I just have some nitpicks:

• from PEP8, you should name your constants using all caps (SAMPLE_RATE, DAMPING);
• you should not abreviate that much your variable names frequency and volume are to be prefered to f and vol;
• you may want to name things to avoid comments: A_major = [220, 220 * 5 // 4, 220 * 6 // 4, 220 * 2].

# Summation

In pure Python world the following code:

value = 0
for element in array:
value += compute(element)
return value


is prefered written as:

return sum(compute(element) for element in array)


As it is both cleaner and faster. In numpy world, it is pretty much the same. You can write:

def generate_chord(frequencies, samples):
plucks = (generate(note, 0.5, samples) for note in frequencies)
return np.sum(plucks, axis=0)


and get the same result.

# Interface

When reading your docstrings, I wondered why you would ask your users (or yourself) to pass in a number of samples if what is more natural to them is to use a duration in seconds. Do not expose implementation details such as your sample rate to the users and let them use what is more natural to them.

I also wonder why you hardcode the volume in generate_chord and let it variable in generate. If the aim is to not have to choose it for each test, you can still use a parameter with default value instead of hardcoding it. It will make things more obvious.

Lastly, I don't really understand the need for integral frequencies as the computations can be performed just as well with floating point values. Just make sure to truncate the computation of N.

I would use the following code as a base for improvements:

import numpy as np

SAMPLE_RATE = 44100
DAMPING = 0.999

def generate(frequency, volume, duration):
"""Generate a Karplus-Strong pluck.

Arguments:
frequency -- the frequency, as a float
volume -- volume, a float in [0.0, 1.0]
duration -- the length of the generated pluck, in seconds.

Return value:
A numpy array of floats with length nsamples.
"""

samples_count = duration * SAMPLE_RATE
N = int(SAMPLE_RATE / frequency)
buf = np.random.rand(N) * 2 - 1
samples = np.empty(samples_count, dtype=float)

for i in range(samples_count):
samples[i] = buf[i % N]
avg = DAMPING * 0.5 * (buf[i % N] + buf[(1 + i) % N])
buf[i % N] = avg

return samples * volume

def generate_chord(frequencies, duration, volume=0.5):
"""Generate a chord of several plucks

Arguments:
frequencies -- a (usually short) list of frequencies
describing a chord.
duration -- the length of the generated pluck, in seconds.

Return value:
A numpy array
"""

plucks = (generate(note, volume, duration) for note in frequencies)
return np.sum(plucks, axis=0)

if __name__ == "__main__":
import scipy.io.wavfile
import matplotlib.pyplot as plt

A_major = [220, 220 * 5 // 4, 220 * 6 // 4, 220 * 2]
strum = generate_chord(A_major, duration=5)

plt.plot(strum)
plt.show()

scipy.io.wavfile.write(filename, SAMPLE_RATE, strum)


# Beyond frequencies

Now, I’m not musician. And even if I were, I would need to look up frequencies of required note as I don't think your average user would know them by heart. It would be best to provide named constants to remember that for us.

Reading at https://en.wikipedia.org/wiki/Musical_note one can easily define something like:

C = 2 ** (-9/12) * 440
C_sharp = 2 ** (-8/12) * 440
D_flat = 2 ** (-8/12) * 440
D = 2 ** (-7/12) * 440
D_sharp = 2 ** (-6/12) * 440
E_flat = 2 ** (-6/12) * 440
E = 2 ** (-5/12) * 440
F = 2 ** (-4/12) * 440
F_sharp = 2 ** (-3/12) * 440
G_flat = 2 ** (-3/12) * 440
G = 2 ** (-2/12) * 440
G_sharp = 2 ** (-1/12) * 440
A_flat = 2 ** (-1/12) * 440
A = 2 ** (0/12) * 440
A_sharp = 2 ** (1/12) * 440
B_flat = 2 ** (1/12) * 440
B = 2 ** (2/12) * 440


But this does not account for every octave. Let's fix that and simplify the reading using an enum:

import enum

def octave(order):
frequency_of_A = 440 * 2 ** (order - 4)

class Octave(enum.Enum):
C = -9
C_sharp = -8
D_flat = -8
D = -7
D_sharp = -6
E_flat = -6
E = -5
F = -4
F_sharp = -3
G_flat = -3
G = -2
G_sharp = -1
A_flat = -1
A = 0
A_sharp = 1
B_flat = 1
B = 2

def __float__(self):
return 2 ** (self.value/12) * frequency_of_A

return Octave

SubSubContra = octave(-1)
SubContra = octave(0)
Contra = octave(1)
Great = octave(2)
Small = octave(3)
OneLined = octave(4)
TwoLined = octave(5)
ThreeLined = octave(6)
FourLined = octave(7)
FiveLined = octave(8)
SixLined = octave(9)


Now I can easily define A Major using:

A_major = [Small.A, OneLined.C_sharp, OneLined.E, OneLined.A]


# Proposed improvements

import enum
import numpy as np

SAMPLE_RATE = 44100
DAMPING = 0.999

def octave(order):
"""Create an enum of all the notes in the octave
of the given order.
"""
frequency_of_A = 440 * 2 ** (order - 4)

class Octave(enum.Enum):
C = -9
C_sharp = -8
D_flat = -8
D = -7
D_sharp = -6
E_flat = -6
E = -5
F = -4
F_sharp = -3
G_flat = -3
G = -2
G_sharp = -1
A_flat = -1
A = 0
A_sharp = 1
B_flat = 1
B = 2

def __float__(self):
return 2 ** (self.value/12) * frequency_of_A

return Octave

SubSubContra = octave(-1)
SubContra = octave(0)
Contra = octave(1)
Great = octave(2)
Small = octave(3)
OneLined = octave(4)
TwoLined = octave(5)
ThreeLined = octave(6)
FourLined = octave(7)
FiveLined = octave(8)
SixLined = octave(9)

def generate(frequency, volume, duration):
"""Generate a Karplus-Strong pluck.

Arguments:
frequency -- the frequency, as a float
volume -- volume, a float in [0.0, 1.0]
duration -- the length of the generated pluck, in seconds.

Return value:
A numpy array of floats with length nsamples.
"""

samples_count = duration * SAMPLE_RATE
N = int(SAMPLE_RATE / frequency)
buf = np.random.rand(N) * 2 - 1
samples = np.empty(samples_count, dtype=float)

for i in range(samples_count):
samples[i] = buf[i % N]
avg = DAMPING * 0.5 * (buf[i % N] + buf[(1 + i) % N])
buf[i % N] = avg

return samples * volume

def generate_chord(frequencies, duration, volume=0.5):
"""Generate a chord of several plucks

Arguments:
frequencies -- a (usually short) list of frequencies
describing a chord.
duration -- the length of the generated pluck, in seconds.

Return value:
A numpy array
"""

notes = map(float, frequencies)
plucks = (generate(note, volume, duration) for note in notes)
return np.sum(plucks, axis=0)

if __name__ == "__main__":
import scipy.io.wavfile
import matplotlib.pyplot as plt

A_major = [Small.A, OneLined.C_sharp, OneLined.E, OneLined.A]
strum = generate_chord(A_major, duration=5)

plt.plot(strum)
plt.show()

scipy.io.wavfile.write(filename, SAMPLE_RATE, strum)

• Thanks. Actually I think damping shouldn't be a global variable at all. The "enum" of named notes is clever. I'm actually looking at just tuning, so the frequencies are not at 440*2**(n/12) but I think it can be adjusted. – James K Feb 16 '17 at 20:48
for i in range(nsamples):
samples[i] = buf[i % N]
avg = damping * 0.5 * (buf[i % N] + buf[(1 + i) % N])
buf[i % N] = avg


There's some definite possibility of optimisation by extracting repeated values.

j, nextj = 0, 1
dampAvg = damping * 0.5
for i in range(nsamples):
sample = buf[j]
samples[i] = sample
buf[j] = dampAvg * (sample + buf[nextj])
j = nextj
nextj = nextj + 1
if nextj == N:
nextj = 0


The next level of optimisation would be to observe that each element of buf is copied to samples and then modified. So if there's a fast way (I assume there is, but I don't know numpy) of appending a copy of buf to samples, you could rewrite as something along the lines of

samples = []
off = 0
while off < nsamples:
if nsamples - off <= N:
copy(buf, 0, samples, off, nsamples - off)
break
copy(buf, 0, samples, off, N)
off += N

for i in range(nsamples - 1):
buf[i] = dampAvg * (buf[i] + buf[i + 1])
buf[nsamples - 1] = dampAvg * (buf[nsamples - 1] + buf)


The step after that would be to look for a vectorised implementation of the averaging loop...