/*
    Simple Karplus-Strong class in C++ by Pieter Suurmond, 2004.
    The specification and the private variables resides in headerfile ks.hxx.
*/

using namespace std;                // Really necessary here?
#include "ks.hxx"                   // Prototypes and private vars.

KS::KS(double length)               // Constructor, expects (total) effective delay length in samples.
{
    if (length < 8.0)                               // Requested total wavelength, (W in ks.hxx).
        throw "Sorry, less than 8 is not allowed.";
    this->del_length = (int)(length - 0.5);         // Determine course-tuning parameter (L in ks.hxx).
    this->c1 = length - 0.5 - (double)del_length;   // Determine fine-tuning parameter(s) (C in ks.hxx).
    this->c0 = 1.0 - c1;                            // We don't want process() to calc this for any sample.
    this->max_length = 1;                   // Find nearest power of 2 greater than the necessary length,
    while (max_length <= del_length)        // so we can use bitmasks for wrapping.
        max_length <<= 1;                   // (If we need a delay of 128 samples, we must allocate 256!)
    this->delay = new double[max_length];   // Allocate memory for large delay.
                                            // It will throw an exeception on failure.
    // Are the above doubles cleared automatically? What does the C++ standard say about it?
    this->bmask = max_length - 1;           // Always use this bitmask for ALL reading and writing!
    this->n = 0;                            // Not really necessary to initialise index n.
    this->lp_fixed = 0.0;                   // Also clear these small intermediate signal buffers.
    this->lp_tuned = 0.0;
    if      (del_length > 512) this->gain = 0.49; // Somewhat less for longer delay lines.
    else if (del_length > 256) this->gain = 0.499;
    else if (del_length > 128) this->gain = 0.4999;
    else if (del_length >  64) this->gain = 0.49999;
    else                       this->gain = 0.499999;
    this->energy_accumulated = 0.0;
}

KS::~KS()                           // Destructor.
{
    delete [] delay;                // Release the buffer we allocated.
}


/* Same diagram as in ks.hxx, now
   marked with signals a, b and c.       (c0=C)
                                            |
                                            v
                   ----------     --------> x ----    (gain)
                  |          |   |                |     |
                  |    -1    v   |    -1          v     v
     x[n] ---> + ---> z  --> + ----> z  --> x --> + --> x -----> y[n]
               ^  a              b          ^              c |
               |                            |                |
               |                         (c1=1-C)            |
               |       -(L)                                  |
                ----- z    <---------------------------------
*/
double KS::process(double x)        // Processes one single incoming sample.
{
    double  a, b, c;                // Variables for easier reading, see diagram above.

    a = x + delay[(n - del_length) & bmask];    // Read from buffer, 100 samples ago.
    b = (a + lp_fixed);                         // Lowpass filter with fixed coefficients. 
    lp_fixed = a;                               // Remember for lowpass filtering the next sample.
    c = ((c0 * b) + 
         (c1 * lp_tuned)) * gain;               // One multiplication too many (can be optimised!).
    lp_tuned = b;                               // Remember for lowpass filtering the next sample.
    delay[n++ & bmask] = c;                     // Write in buffer; increment time; return result.

    // Non-linear extension of the simple model: string starts to loosen after a while...
    if ((energy_accumulated > 500) && (del_length < max_length-1))
        {
        if (c1 < 0.9)       // Relief string tension bit by bit (i.e. interpolating) 
            {               // so that we won't hear sudden 'clicks': 
            c1 += 0.1;      // +0.1 sample detune after collecting 500 units of energy.
            c0 = 1.0 - c1;
            }
        else
            {
            c1 = 0.0;
            c0 = 1.0 /* - c1 */;
            del_length++;
            }
        energy_accumulated = 0.0;
        }
    else
        energy_accumulated += c * c;

    return c;
}