/* Try a somewhat realistic simulation of a piano note, by hand, with
 * Karplus-Strong.  The result sounds good, and sounds like a string
 * instrument, but not like a piano.
 *
 * Produces 48kHz 16-bit little-endian 2-channel signed output, like a DAT.
 *
 make -k ks-string && ./ks-string | sox -t raw -r 48000 -e signed-integer -b 16 -L -c 2 - guitars.wav
 *
 * The wrinkles on top of the basic Karplus-Strong algorithm, which is
 * one line of code, are:
 *
 * - There are three separate delay lines with slightly separate
 *   delays, in this case implemented by tweaking the smoothing filter
 *   to be a crude fractional-delay filter as well.  (That’s a thing,
 *   right?)
 *
 * - The output is not just a single tap on the delay line, but two
 *   taps about one fourteenth of the buffer apart, to suppress
 *   seventh harmonics, like the piano hammer striking the string one
 *   seventh of the way down.
 *
 * - I haven't gotten this working well yet, but the initial impulse
 *   ought to have a spectrum that omits the high frequencies, like
 *   probably everything above about 4kHz ought to be super weak,
 *   because the actual hammer contact with the strings takes place
 *   over a millisecond or two.
 */

#include <stdio.h>
#include <stdint.h>

typedef uint8_t u8;
typedef uint16_t u16;
typedef uint32_t u32;
typedef int32_t  s32;

enum {
  buflen = 219,             /* 48000 / 220 + 1, more or less */
  n_resonators = 3,
  decay_bits = 5,
};

int main() {
  /* We only use the low 24 bits of these in order to be able to use
     32-bit math for weighted sums.
   */
  s32 buf[n_resonators][buflen] = {{0}};
  u32 x = 34232;                /* random generator */
  s32 decay[4] = {0};           /* for smoothing the random noise */

  for (u32 i = 0; i < 48000 * 5; i++) {
    u16 n = i % buflen;
    u16 n1 = (i + 1) % buflen;
    u16 n2 = (i + 2) % buflen;
    s32 e = 0;
    if (i < (buflen - 1) * 4) {
      /* Attack! */
      x = x * 1103515245 + 12345;
      /* x is a 32-bit random number.  We want a 24-bit *signed*
         random number, but actually one limited to a quarter of the
         maximum of its range: so -2²¹ to 2²¹
      */
      decay[0] = (x >> 10) - (1 << 21);
      for (u8 j = 0; j < 3; j++) {
        decay[j+1] -= decay[j+1] >> decay_bits;
        decay[j+1] += decay[j] >> decay_bits;
      }
      e = decay[3] << 4;
      /* e = (i > 96 ? 192 - i : i) << 15; */
    }

    /* These three FIR filters on the feedback are intended to all be
     * low-pass, slightly below unity gain
     * and also give different fractional delays, by a hertz
     * or so at the fundamental. */
    buf[0][n] = ( 256 * e
                + 64 * buf[0][n]
                + 127 * buf[0][n1]
                + 64 * buf[0][n2]) / 256;
    buf[1][n] = (  256 * e
                 + 113 * buf[1][n]
                 + 102 * buf[1][n1]
                 +  40 * buf[1][n2]) / 256;
    buf[2][n] = (  256 * e
                 +   8 * buf[2][n]
                 + 106 * buf[2][n1]
                 + 142 * buf[2][n2]) / 256;
    u16 t = (i + 26) % buflen;
    s32 s
        = (85 * buf[0][n] + 85 * buf[0][t]
         + 85 * buf[1][n] + 85 * buf[1][t]
         + 85 * buf[2][n] + 85 * buf[2][t]) / 256;
    putchar(s >> 8); putchar(s >> 16);
    putchar(s >> 8); putchar(s >> 16);
  }
  return 0;
}