Okay, I finally found my copy of Ken Steiglitz’s A DSP Primer (a great book, but sadly more expensive now than when I got my copy) and read through the implementation of the tunable plucked string instrument. A couple of things really need to be added: first of all, I was off considerably in my tuning because the averaging operation inserts a delay of 1/2 a sample, and the tuning was necessarily coarse because the frequency was entirely determined by the length of the delay buffer, which was an integer. To make that work, you need to introduce an all-pass filter, which has a pretty simple form. Without any additional explanation, here’s the source code that I tinkered together over a half an hour of reading and typing.

[sourcecode lang=”C”]
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sndfile.h>

/*
* ks.c
* an implementation of the Karplus-Strong algorithm, revised after
* reading Ken Steiglitz’s treatment in "A DSP Primer".
*/

#define SAMPLE_RATE (22050)

double freq = 440. ;

main(int argc, char *argv[])
{
SNDFILE *sf ;
SF_INFO sfinfo ;
double delay ;
int L ;
int i, j ;
int sample ;
double a ;
double w0 ;
double x0, y0, x1, y1 ;

freq = atof(argv[1]) ;

sfinfo.samplerate = SAMPLE_RATE ;
sfinfo.channels = 1 ;
sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16 ;

sf = sf_open(argv[2], SFM_WRITE, &sfinfo) ;

/* first of all, figure out the total delay…
* freq = SAMPLE_RATE / (L + 0.5)
* or L = SAMPLE_RATE / freq – 0.5 ;
*/

delay = SAMPLE_RATE / freq – 0.5 ;
fprintf(stderr, "… delay = %f\n", delay) ;
L = floor(delay) ;
delay -= L ;
fprintf(stderr, "… buffer size is = %d\n", L) ;
fprintf(stderr, "… fractional delay is = %f\n", delay) ;

fprintf(stderr, "… approximate value for a = %f\n", (1-delay)/(1+delay)) ;
w0 = 2.0 * M_PI * freq / SAMPLE_RATE ;
a = sin((1. – delay) * w0/2.) / (sin((1. + delay) * w0/2.)) ;
fprintf(stderr, "… exact value for a = %f\n", a) ;

/* okay, now generate one second of the plucked string sound.
*/

double *buf = calloc(L, sizeof(double)) ;

/* initialize with random numbers. */
for (i=0; i<L; i++)
buf[i] = 2.0 * drand48() – 1.0 ;

x0 = y0 = 0. ;
for (sample=i=0; sample < SAMPLE_RATE * 0.25 ; sample++) {
j = i + 1 ;
if (j >= L) j = 0 ;
x1 = (buf[i] + buf[j]) / 2. ;
/* implement the "allpass" filter. */
y1 = a * x1 + x0 – a * y0 ;
sf_write_double(sf, &y1, 1) ;
y0 = y1 ;
x0 = x1 ;
buf[i] = y1 ;
i = j ;
}

sf_close(sf) ;
}
[/sourcecode]

This just writes out 1/4 of a second of a note tuned at the frequency that you specify. I wrote a little python program, and computed an mp3 of a simple two octave run of notes.

Addendum: I synthesized the small set of three note chords just to test how them merged. They do sound a bit plinky and a little thin, but not terrible.