/* 
 * guttersynth.c
 * Part of libguttersynth
 * 
 * Copyright Tom Mudd 2019, Richard Knight 2021, 2022
 * 
 * Ported to C by Richard Knight
 * from https://github.com/tommmmudd/guttersynthesis by Tom Mudd
 * 
 */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <sys/param.h>
#include "guttersynth.h"


// TODO: freqs temp and Q temp, remove, not needed..
//

/*
 * Returns a random float within a range
 *
 * min: bottom end of range
 * max: top end of range
 * 
 * returns: random float
 */
static FLT rangerandom(FLT min, FLT max) 
{
    FLT random = ((FLT) rand()) / (FLT) RAND_MAX;
    FLT diff = max - min;
    FLT r = random * diff;
    return min + r;
}


/*
 * Calculate filter coefficients
 * 
 * s: gutter_state struct to apply operation to
 */
void gutter_calccoeffs(gutter_state* s) 
{
    int i, b, f;
    for (b = 0; b < s->bankCount; b++) {
        for (f = 0; f < s->filterCount; f++) {
            i = b * s->bankCount + f;
            s->gi.V[i] = (FLT) pow(10.0, 0.05);
            s->gi.K[i] = (FLT) tan(((FLT)M_PI * s->filterFreqs[i]) / s->gi.samplerate);
            s->gi.norm[i] = (1.0 / (1.0 + s->gi.K[i] / s->Q[i] + s->gi.K[i] * s->gi.K[i]));
            s->gi.a0[i] = s->gi.K[i] / s->Q[i] * s->gi.norm[i];
            s->gi.a1[i] = 0.0;
            s->gi.a2[i] = - s->gi.a0[i];
            s->gi.b1[i] = 2.0 * (s->gi.K[i] * s->gi.K[i] - 1.0) * s->gi.norm[i];
            s->gi.b2[i] = (1.0 - s->gi.K[i] / s->Q[i] + s->gi.K[i] * s->gi.K[i]) *s->gi.norm[i];
        }
    }
}


/*
 * Apply distortion
 */
static void gutter_distortion(gutter_state* s, FLT* audio_ptr) 
{
    FLT input = *audio_ptr;
    FLT output;
    switch (s->distortionMethod) {
        case 0:
            output = MAX(MIN(input, 1), -1);
            break;
        case 1:
            if (s->gi.finalY <= -1) {
                output = -0.666666667;
            } else if (input <= 1) {
                output = input - input * input * input / 3.0;
            } else {
                output = 0.666666667;
            }
            break;
        case 2:
            output = (FLT) atan(input);
            break;
        case 3:
            output = 0.75 * (((FLT)sqrt(input * 1.3 * (input * 1.3) + 1.0)) * 1.65 - 1.65) / input;
            break;
        case 4:
            output = (0.1076 * input * input * input + 3.029 * input) / (input * input + 3.124);
            break;
        case 5:
            output = 2.0 / (1.0 + ((FLT)exp(-1.0 * input)));
            break;
        default:
            output = input;
    }
    *audio_ptr = output;
}


/*
 * Randomise filter bank frequencies
 * 
 * s: gutter_state struct to apply operation to
 */
void gutter_randomisefilters(gutter_state* s)
{
    int i, b, f;
    for (b = 0; b < s->bankCount; b++) {
        for (f = 0; f < s->filterCount; f++) {
            i = b * s->bankCount + f;
            s->filterFreqs[i] = rangerandom(100, 4000);
            s->gains[i] = rangerandom(0, 1);
            s->Q[i] = rangerandom(20, 200);
        }
    }
    gutter_calccoeffs(s);
}


int gutter_setdistortionmethod(gutter_state* s, int method)
{
    if (method < 0 || method > 5) {
        return -1;
    }
    s->distortionMethod = method;
    return 0;
}

/*
 * Set frequency of a filter in a bank
 * 
 * s: gutter_state struct to apply operation to
 * bank: bank index
 * filter: filter index
 * value: frequency to set
 * 
 * returns: 0 on success, -1 if bank or filter is out of range
 */
int gutter_setfreq(gutter_state* s, int bank, int filter, FLT value) 
{
    if (bank >= s->bankCount || filter >= s->filterCount || bank < 0 || filter < 0) {
        return -1;
    }
    s->filterFreqs[bank*s->bankCount+filter] = value;
    return 0;
}


/*
 * Set Q of a filter in a bank
 * 
 * s: gutter_state struct to apply operation to
 * bank: bank index
 * filter: filter index
 * value: Q to set
 * 
 * returns: 0 on success, -1 if bank or filter is out of range
 */
int gutter_setq(gutter_state* s, int bank, int filter, FLT value)
{
    if (bank >= s->bankCount || filter >= s->filterCount || bank < 0 || filter < 0) {
        return -1;
    }
    s->Q[bank*s->bankCount+filter] = value;
    return 0;
}


/*
 * Set gain of a filter in a bank
 * 
 * s: gutter_state struct to apply operation to
 * bank: bank index
 * filter: filter index
 * value: gain to set
 * 
 * returns: 0 on success, -1 if bank or filter is out of range
 */
int gutter_setgain(gutter_state* s, int bank, int filter, FLT value)
{
    if (bank >= s->bankCount || filter >= s->filterCount || bank < 0 || filter < 0) {
        return -1;
    }
    s->gains[bank*s->bankCount+filter] = value;
    return 0;
}


/*
 * Set the Q of all filters in all banks
 * 
 * s: gutter_state struct to apply operation to
 * value: Q to set
 */
void gutter_setqall(gutter_state* s, FLT value)
{
    int i;
    for (i = 0; i < s->bankCount * s->filterCount; i++) {
        s->Q[i] = value;
    }
}

int gutter_setqbank(gutter_state* s, int bank, FLT value)
{
    if (bank < 0 || bank <= s->bankCount) {
        return -1;
    }
    int f;
    for (f = 0; f < s->filterCount; f++) {
        s->Q[bank * s->bankCount + f] = value;
    }
}

/*
 * Initialise with custom memory allocator and deallocator
 * 
 * bankCount: number of banks (typically 1 or 2)
 * filterCount: number of filters (typically up to 24)
 * samplerate: samplerate to calculate audio at
 * allocator: memory allocator function taking a size_t and returning void*
 * deallocator: memory deallocator function taking a void*
 * 
 * returns: gutter_state struct to be used in the synthesis session
 */
gutter_state* gutter_init_ca(int bankCount, int filterCount, int samplerate, 
        void* (*allocator)(size_t), void (*deallocator)(void*)) 
{
    int size2d, i, f, b;
    
    srand(time(0));
    
    gutter_state* s = (gutter_state*) allocator(sizeof(gutter_state));
    
    s->gi.samplerate = samplerate;
    s->bankCount = bankCount;
    s->filterCount = filterCount;
    s->gi.dealloc = deallocator;
    
    
    s->gi.dcblock.xm1 = 0;
    s->gi.dcblock.ym1 = 0;
    s->gi.dcblock.r = 0.995;
    
    
    size2d = sizeof(FLT) * bankCount * filterCount;
    s->gains = (FLT*) allocator(size2d);
    s->Q = (FLT*) allocator(size2d);
    s->filterFreqs = (FLT*) allocator(size2d);
    s->gi.prevX1 = (FLT*) allocator(size2d);
    s->gi.prevX2 = (FLT*) allocator(size2d);
    s->gi.prevY1 = (FLT*) allocator(size2d);
    s->gi.prevY2 = (FLT*) allocator(size2d);
    s->gi.V = (FLT*) allocator(size2d);
    s->gi.K = (FLT*) allocator(size2d);
    s->gi.norm = (FLT*) allocator(size2d);
    s->gi.a0 = (FLT*) allocator(size2d);
    s->gi.a1 = (FLT*) allocator(size2d);
    s->gi.a2 = (FLT*) allocator(size2d);
    s->gi.b1 = (FLT*) allocator(size2d);
    s->gi.b2 = (FLT*) allocator(size2d);
    s->gi.y = (FLT*) allocator(size2d);
    
    s->filtersOn = 1;
    s->smoothing = 1;
    s->distortionMethod = 2;
    s->enableAudioInput = 0;
    
    for (b = 0; b < bankCount; b++) {
        for (f = 0; f < filterCount; f++) {
            i = b * bankCount + f;
            FLT freq = ((FLT)(f + 1)) / 2.0 * 20.0 * (b + 1) * 1.2 + 80;
            s->filterFreqs[i] = freq;
            s->gi.y[i] = 0;
            s->gi.prevX1[i] = 0;
            s->gi.prevX2[i] = 0;
            s->gi.prevY1[i] = 0;
            s->gi.prevY2[i] = 0;
            s->gains[i] = 1;
            s->Q[i] = 1;
        }
    }

    gutter_randomisefilters(s);
    
    return s;
}


/*
 * Initialise with default memory allocator and deallocator (malloc and free)
 * 
 * bankCount: number of banks (typically 1 or 2)
 * filterCount: number of filters (typically up to 24)
 * samplerate: samplerate to calculate audio at
 * 
 * returns: gutter_state struct to be used in the synthesis session
 */
gutter_state* gutter_init(int bankCount, int filterCount, int samplerate)
{
    return gutter_init_ca(bankCount, filterCount, samplerate, &malloc, &free);
}


/*
 * Clean up state: free memory
 * 
 * s: gutter_state struct to apply operation to
 */
void gutter_cleanup(gutter_state* s) 
{
    void (*da)(void*) = s->gi.dealloc;
    
    da(s->gains);
    da(s->filterFreqs);
    da(s->gi.prevX1);
    da(s->gi.prevX2);
    da(s->gi.prevY1);
    da(s->gi.prevY2);
    da(s->gi.V);
    da(s->gi.K);
    da(s->gi.norm);
    da(s->gi.a0);
    da(s->gi.a1);
    da(s->gi.a2);
    da(s->gi.b1);
    da(s->gi.b2);
    da(s->gi.y);
    da(s->Q);
    da(s);
}


/*
 * Reset oscillator
 * 
 * s: gutter_state struct to apply operation to
 */
void gutter_reset(gutter_state* s) 
{
    s->gi.duffX = 0;
    s->gi.duffY = 0;
    s->gi.dx = 0;
    s->gi.dy = 0;
    s->gi.t = 0;
}


/*
 * Apply DC blocking to audio referencing current state
 *
 * s: gutter_state struct to apply operation to
 * input: audio to apply DC block to
 * 
 * returns: DC blocked audio
 */
static FLT dcblock(gutter_state* s, FLT input) {
    FLT y = input - s->gi.dcblock.xm1 + s->gi.dcblock.r * s->gi.dcblock.ym1;
    s->gi.dcblock.xm1 = input;
    s->gi.dcblock.ym1 = y;
    return y;
}

/*
 * Synthesise one sample with audio input
 *
 * s: gutter_state struct to apply operation to
 * audioInput: audio sample to feed into synthesis
 * 
 * returns: synthesised audio sample
 */
FLT gutter_process_input(gutter_state* s, FLT audioInput) 
{
    if (s->filtersOn) {
        int i, b, f;
        for (b = 0; b < s->bankCount; b++) {
            for (f = 0; f < s->filterCount; f++) {
                i = b * s->bankCount + f;
                s->gi.y[i] = s->gi.a0[i] * s->gi.duffX 
                        + s->gi.a1[i] * s->gi.prevX1[i] 
                        + s->gi.a2[i] * s->gi.prevX2[i] 
                        - s->gi.b1[i] * s->gi.prevY1[i] 
                        - s->gi.b2[i] * s->gi.prevY2[i];
                s->gi.prevX2[i] = s->gi.prevX1[i];
                s->gi.prevX1[i] = s->gi.duffX;
                s->gi.prevY2[i] = s->gi.prevY1[i];
                s->gi.prevY1[i] = s->gi.y[i];
                s->gi.finalY += s->gi.y[i] * s->gains[i] * s->singleGain;
            }
        }
    } else {
        s->gi.finalY = s->gi.duffX;
    }

    if (s->enableAudioInput) {
        s->gi.dy = (s->gi.finalY - s->gi.finalY * s->gi.finalY * s->gi.finalY 
                - s->c * s->gi.duffY + s->gamma * audioInput);
    } else {
        s->gi.dy = (s->gi.finalY - s->gi.finalY * s->gi.finalY * s->gi.finalY 
                - s->c * s->gi.duffY + s->gamma * ((FLT)sin(s->omega * s->gi.t)));

    }

    s->gi.duffY += s->gi.dy;
    s->gi.dx = s->gi.duffY;

    FLT out;
    s->gi.duffX = ((s->gi.finalY + s->gi.dx) - s->gi.duffX) / s->smoothing; // low pass filter
    //s->gi.finalY = dcblock(s, s->gi.finalY); // RK: experimental
    if (s->filtersOn) {
        gutter_distortion(s, &(s->gi.duffX));
        out = s->gi.finalY * 0.125;
    } else {
        s->gi.duffX = MAX(MIN(s->gi.duffX, 100), -100);
        if (abs(s->gi.duffX) > 99) {
            gutter_reset(s);
        }
        out = MAX(MIN(s->gi.duffX * s->singleGain, 1.0), -1.0);
    }
    
    s->gi.t += s->dt;
    
    if (isnan(s->gi.duffX)) {
        gutter_reset(s);
    }

    //return out;
    return dcblock(s, out);
}


/*
 * Synthesise specified number of samples with audio input
 *
 * s: gutter_state struct to apply operation to
 * audioInput: audio sample pointer to feed into synthesis, must be at least nsamps size
 * audioOutput: audio output pointer, must be at least nsamps size
 * nsamps: number of samples to synthesise
 */
void gutter_process_input_samples(gutter_state* s, FLT* audioInput, FLT* audioOutput, int nsamps) 
{
    int i;
    for (i = 0; i < nsamps; i++) {
        audioOutput[i] = gutter_process_input(s, audioInput[i]);
    }
}


/*
 * Synthesise one sample
 *
 * s: gutter_state struct to apply operation to
 * 
 * returns: synthesised audio sample
 */
FLT gutter_process(gutter_state* s)
{
    return gutter_process_input(s, 0);
}


/*
 * Synthesise specified number of samples
 *
 * s: gutter_state struct to apply operation to
 * audioOutput: audio output pointer, must be at least nsamps size
 * nsamps: number of samples to synthesise
 */
void gutter_process_samples(gutter_state* s, FLT* audioOutput, int nsamps)
{   
    int i;
    for (i = 0; i < nsamps; i++) {
        audioOutput[i] = gutter_process_input(s, 0);
    }
}