ladspa-missing-fundamental.c

/* Doesn't actually do a missing-fundamental effect
 * I couldn't work out how to do it.
 * Hopefully this is illustrative of how to use LADSPA and/or FFTW3.
 * Not that the latter is used properly.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define _USE_MATH_DEFINES
#include <math.h>

#include <ladspa.h>
#include <fftw3.h>

/* Unique plugin ID: luna0
 * 00 = zero
 * 01 = 'a'
 * 02 = 'b'
 * etc
*/

#define ID_STEREO 1221140100

/* ports */
#define FILTER_CUTOFF    0
#define N_HARMONICS      1
#define HARMONIC_ROLLOFF 2
#define L_INPUT_PORT     3
#define R_INPUT_PORT     4
#define L_OUTPUT_PORT    5
#define R_OUTPUT_PORT    6
#define PORTCOUNT_STEREO (R_OUTPUT_PORT+1)

typedef struct{
    float z1;
    float z2;
} LPF_State;

typedef struct{
    float a0, a1; /* a2 = a0 for LP biquad */
    float b1, b2;
} LPF_Coeff;

void calculateLPFCoeff(LPF_Coeff* coeff, const LADSPA_Data freq, const LADSPA_Data rate){
    const float K     = tanf(M_PI * freq / rate);
    const float K2    = K * K;
    const float Q     = 0.70710678118f; /* butterworth (1/sqrt(2)) */
    const float denom = Q + K2;
    const float norm  = 1.0f / (1.0f + K / denom);

    const float a0 = K2 * norm;
    coeff->a0 = a0;
    coeff->a1 = 2.0f * a0;
    coeff->b1 = 2.0f * (K2 - 1.0f) * norm;
    coeff->b2 = (1.0f - K / denom) * norm;
}

void calculateLPF(LPF_State *const state, const LPF_Coeff *const coeff, const LADSPA_Data *input, double *output, const unsigned long sample_count){
    unsigned long sample = 0;
    for(; sample < sample_count; sample++){
        const float in  = *(input++);
        const float out = in * coeff->a0 + state->z1;
        *(output++)     = out;
        state->z1 = in * coeff->a1 + state->z2 - coeff->b1 * out;
        state->z2 = in * coeff->a0 - coeff->b2 * out; /* a2 == a0 */
    }
}

typedef struct {
    LADSPA_Data  sample_rate;

    LPF_State    l_lpf;
    LPF_State    r_lpf;
    /* we can't generate a plan because we don't know the buffer size
     * really, we should set a small window size and then overlap and
     * do STFT, but this is where I got bored with this.
    fftw_plan    fft_plan;
    fftw_plan    ifft_plan;
    */
    
    /* ports */
    LADSPA_Data* filter_cutoff;
    LADSPA_Data* n_harmonics;
    LADSPA_Data* harmonic_rolloff;
    LADSPA_Data* l_input;
    LADSPA_Data* r_input;
    LADSPA_Data* l_output;
    LADSPA_Data* r_output;
} MissingFundamental;

LADSPA_Handle instantiateMissingFundamental(const LADSPA_Descriptor* descriptor, unsigned long sample_rate){
    MissingFundamental* missing_fundamental;

    missing_fundamental =
        (MissingFundamental*) malloc(sizeof(MissingFundamental));

    missing_fundamental->sample_rate =
        (LADSPA_Data) sample_rate;

    return missing_fundamental;
}

void cleanupMissingFundamental(LADSPA_Handle instance){
    free(instance);
}

void activateMissingFundamental(LADSPA_Handle instance){
}

void connectPortToMissingFundamental(LADSPA_Handle instance, unsigned long port, LADSPA_Data* data_location){
    MissingFundamental* mf = (MissingFundamental*) instance;

    switch(port){
        case FILTER_CUTOFF:
            mf->filter_cutoff = data_location;
            break;
        case N_HARMONICS:
            mf->n_harmonics = data_location;
            break;
        case HARMONIC_ROLLOFF:
            mf->harmonic_rolloff = data_location;
            break;
        case L_INPUT_PORT:
            mf->l_input  = data_location;
            break;
        case R_INPUT_PORT:
            mf->r_input  = data_location;
            break;
        case L_OUTPUT_PORT:
            mf->l_output = data_location;
            break;
        case R_OUTPUT_PORT:
            mf->r_output = data_location;
            break;
    }
}

void transformFFT(const fftw_complex* fft, fftw_complex* ifft, const unsigned long sample_count){
    for(unsigned long bin = 0; bin < sample_count; bin++){
        if(bin > 4 && bin < sample_count/2){
            ifft[bin][0] = fft[bin][0];
            ifft[bin][1] = fft[bin][1];
        }
        else{
            ifft[bin][0] = ifft[bin][1] = 0.0;
        }
    }
}

void runMissingFundamental(LADSPA_Handle instance, unsigned long sample_count){
    MissingFundamental* mf = (MissingFundamental*) instance;

    const LADSPA_Data sample_rate = mf->sample_rate;

    LADSPA_Data* filter_cutoff    = mf->filter_cutoff;
    LADSPA_Data* n_harmonics      = mf->n_harmonics;
    LADSPA_Data* harmonic_rolloff = mf->harmonic_rolloff;
    LADSPA_Data* l_output         = mf->l_output;
    LADSPA_Data* r_output         = mf->r_output;
    LADSPA_Data* l_input_start    = mf->l_input;
    LADSPA_Data* l_input          = mf->l_input;
    LADSPA_Data* r_input_start    = mf->r_input;
    LADSPA_Data* r_input          = mf->r_input;

    double *const l_buffer_start   =
        (double*) fftw_malloc(sample_count * sizeof(double));
    double *const r_buffer_start   =
        (double*) fftw_malloc(sample_count * sizeof(double));
    double* l_buffer         = l_buffer_start;
    double* r_buffer         = r_buffer_start;

    fftw_complex* fft_buffer  =
        (fftw_complex*) fftw_malloc(sample_count * sizeof(fftw_complex));
    fftw_complex* ifft_buffer =
        (fftw_complex*) fftw_malloc(sample_count * sizeof(fftw_complex));

    /* generating new FFTW3 plans is probably extremely slow *
     * but this does work okay on my system                  */
    fftw_plan fft_plan = fftw_plan_dft_r2c_1d(
            sample_count,
            l_buffer_start,
            fft_buffer,
            FFTW_ESTIMATE);
    fftw_plan ifft_plan = fftw_plan_dft_c2r_1d(
            sample_count,
            ifft_buffer,
            l_buffer_start,
            FFTW_ESTIMATE);

    LPF_State*   l_lpf            = &mf->l_lpf;
    LPF_State*   r_lpf            = &mf->r_lpf;
    LPF_Coeff    lpf_coeff        = {0};
    
    calculateLPFCoeff(&lpf_coeff, *filter_cutoff, sample_rate);
    
    calculateLPF(l_lpf, &lpf_coeff, l_input, l_buffer, sample_count);
    calculateLPF(r_lpf, &lpf_coeff, r_input, r_buffer, sample_count);
    
    fftw_execute(fft_plan);
    transformFFT(fft_buffer, ifft_buffer, sample_count);
    fftw_execute(ifft_plan);

    fftw_execute_dft_r2c(fft_plan, r_buffer_start, fft_buffer);
    transformFFT(fft_buffer, ifft_buffer, sample_count);
    fftw_execute_dft_c2r(ifft_plan, ifft_buffer, r_buffer_start);

    const double descaling = 1.0 / sample_count; /* un-scaling for FFTW3 */

    unsigned long sample = 0;
    for(sample   = 0,
        l_input  = l_input_start,
        r_input  = r_input_start,
        l_buffer = l_buffer_start,
        r_buffer = r_buffer_start;
        sample < sample_count;
        sample++){
        *(l_output++) = *(l_input++)*0.2 + (*l_buffer++)*0.3*descaling;
        *(r_output++) = *(r_input++)*0.2 + (*r_buffer++)*0.3*descaling;
    }

    fftw_free(l_buffer_start);
    fftw_free(r_buffer_start);
    fftw_free(ifft_buffer);
    fftw_free(fft_buffer);
}

LADSPA_Descriptor *descriptor = NULL; /* must be visible to _init and _fini */
void _init(){
    LADSPA_PortDescriptor  *port_descs = NULL;
    LADSPA_PortRangeHint   *port_hints = NULL;
    char                  **port_names = NULL;

    descriptor =
        (LADSPA_Descriptor*) malloc(sizeof(LADSPA_Descriptor));
    if(descriptor){
        descriptor->UniqueID   = ID_STEREO;
        descriptor->Label      = strdup("MissingFundamental");
        descriptor->Name       = strdup("Missing Fundamental");
        descriptor->Maker      = strdup("Luna");
        descriptor->Copyright  = strdup("AGPLv3");
        descriptor->Properties = LADSPA_PROPERTY_HARD_RT_CAPABLE;
        descriptor->PortCount  = PORTCOUNT_STEREO;

        descriptor->instantiate  = instantiateMissingFundamental;
        descriptor->connect_port = connectPortToMissingFundamental;
        descriptor->activate     = activateMissingFundamental;
        descriptor->run          = runMissingFundamental;
        descriptor->cleanup      = cleanupMissingFundamental;
        
        descriptor->run_adding          = NULL;
        descriptor->set_run_adding_gain = NULL;
        descriptor->deactivate          = NULL;

        port_descs =
            (LADSPA_PortDescriptor*) calloc(PORTCOUNT_STEREO,
                                            sizeof(LADSPA_PortDescriptor));
        port_descs[FILTER_CUTOFF   ] = LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL;
        port_descs[N_HARMONICS     ] = LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL;
        port_descs[HARMONIC_ROLLOFF] = LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL;
        port_descs[L_INPUT_PORT    ] = LADSPA_PORT_INPUT  | LADSPA_PORT_AUDIO;
        port_descs[R_INPUT_PORT    ] = LADSPA_PORT_INPUT  | LADSPA_PORT_AUDIO;
        port_descs[L_OUTPUT_PORT   ] = LADSPA_PORT_OUTPUT | LADSPA_PORT_AUDIO;
        port_descs[R_OUTPUT_PORT   ] = LADSPA_PORT_OUTPUT | LADSPA_PORT_AUDIO;
        descriptor->PortDescriptors = port_descs;

        port_names =
            (char**) calloc(PORTCOUNT_STEREO, sizeof(char*));
        port_names[FILTER_CUTOFF   ] = strdup("Filter Cutoff");
        port_names[N_HARMONICS     ] = strdup("Number of harmonics");
        port_names[HARMONIC_ROLLOFF] = strdup("Harmonic rolloff");
        port_names[L_INPUT_PORT    ] = strdup("Input");
        port_names[R_INPUT_PORT    ] = strdup("Input");
        port_names[L_OUTPUT_PORT   ] = strdup("Output");
        port_names[R_OUTPUT_PORT   ] = strdup("Output");
        descriptor->PortNames = (const char *const *)port_names;

        port_hints =
            (LADSPA_PortRangeHint*) calloc(PORTCOUNT_STEREO, sizeof(LADSPA_PortRangeHint));
        
        port_hints[FILTER_CUTOFF].HintDescriptor =
            LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE;
        port_hints[FILTER_CUTOFF].LowerBound = 20.0;
        port_hints[FILTER_CUTOFF].UpperBound = 11025.0;
        
        port_hints[N_HARMONICS].HintDescriptor =
            LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE;
        port_hints[N_HARMONICS].LowerBound = 0.0;
        port_hints[N_HARMONICS].UpperBound = 1024.0;

        port_hints[HARMONIC_ROLLOFF].HintDescriptor =
            LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE;
        port_hints[HARMONIC_ROLLOFF].LowerBound = 0.0;
        port_hints[HARMONIC_ROLLOFF].UpperBound = 8.0;

        port_hints[L_INPUT_PORT ].HintDescriptor =
        port_hints[R_INPUT_PORT ].HintDescriptor =
        port_hints[L_OUTPUT_PORT].HintDescriptor =
        port_hints[R_OUTPUT_PORT].HintDescriptor =
            0;

        descriptor->PortRangeHints = port_hints;
    }
}

void _fini(){
    if(descriptor){
        free((void*)descriptor->Label);
        free((void*)descriptor->Name);
        free((void*)descriptor->Maker);
        free((void*)descriptor->Copyright);
        free((void*)descriptor->PortDescriptors);
        for(unsigned long port = 0; port < descriptor->PortCount; port++){
            free((void*) descriptor->PortNames[port]);
        }
        free((void*)descriptor->PortRangeHints);
        free(descriptor);
        descriptor = NULL;
    }
}

const LADSPA_Descriptor* ladspa_descriptor(unsigned long index){
    if(index == 0){
        return descriptor;
    }
    else{
        return NULL;
    }
}