autotune~.c

/* 
 autotune.c
 
 An auto-tuning PD External, based on 
 autotalent an auto-tuning LADSPA plugin and an older port by Maxus Germanus
 
 Free software by Thomas A. Baran.
 http://web.mit.edu/tbaran/www/autotune.html
 VERSION 0.2

 Ivica Ico Bukvic <ico.bukvic.net>
 VERSION 0.9.1
 changes:
     *added ability to specify FFT size
     *added ability to specify hop (window) size--requires further refinement
      to allow for more nimble pitch readjustment
     *added pull 2 value where pitch pull is relative to the closest pitch
     *circumvented segfaults due to NANs--may require a better fix

 VERSION 0.9
 changes:
     *compatible with the latest version of autotalent
 
 This program is free software; you can redistribute it and/or modify        
 it under the terms of the GNU General Public License as published by        
 the Free Software Foundation; either version 2 of the License, or           
 (at your option) any later version.                                         
 
 This program is distributed in the hope that it will be useful,             
 but WITHOUT ANY WARRANTY; without even the implied warranty of              
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               
 GNU General Public License for more details.                                
 
 You should have received a copy of the GNU General Public License           
 along with this program; if not, write to the Free Software                 
 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "m_pd.h"

#ifndef MAYER_H
#define MAYER_H
#define REAL float
void mayer_realfft(int n, REAL *real);
void mayer_realifft(int n, REAL *real);
#endif

#define PI (float)3.14159265358979323846
#define L2SC (float)3.32192809488736218171

#ifndef CLIP
#define CLIP(a, lo, hi) ( (a)>(lo)?( (a)<(hi)?(a):(hi) ):(lo) )
#endif

// Variables for FFT routine
typedef struct
{
	int nfft;        // size of FFT
	int numfreqs;    // number of frequencies represented (nfft/2 + 1)
	float* fft_data; // array for writing/reading to/from FFT function
} fft_vars;

// Constructor for FFT routine
fft_vars* fft_con(int nfft)
{
	fft_vars* membvars = (fft_vars*) malloc(sizeof(fft_vars));

	membvars->nfft = nfft;
	membvars->numfreqs = nfft/2 + 1;

	membvars->fft_data = (float*) calloc(nfft, sizeof(float));

	return membvars;
}

// Destructor for FFT routine
void fft_des(fft_vars* membvars)
{
	free(membvars->fft_data);

	free(membvars);
}

// Perform forward FFT of real data
// Accepts:
//   x - pointer to struct of FFT variables
//   input - pointer to an array of (real) input values, size nfft
//   output_re - pointer to an array of the real part of the output,
//     size nfft/2 + 1
//   output_im - pointer to an array of the imaginary part of the output,
//     size nfft/2 + 1
void fft_forward(fft_vars* membvars, float* input, float* output_re, float* output_im)
{
	int ti;
	int nfft;
	int hnfft;
	int numfreqs;

	nfft = membvars->nfft;
	hnfft = nfft/2;
	numfreqs = membvars->numfreqs;

	for (ti=0; ti<nfft; ti++) {
		membvars->fft_data[ti] = input[ti];
	}

	mayer_realfft(nfft, membvars->fft_data);

	output_im[0] = 0;
	for (ti=0; ti<hnfft; ti++) {
		output_re[ti] = membvars->fft_data[ti];
		output_im[ti+1] = membvars->fft_data[nfft-1-ti];
	}
	output_re[hnfft] = membvars->fft_data[hnfft];
	output_im[hnfft] = 0;
}

// Perform inverse FFT, returning real data
// Accepts:
//   x - pointer to struct of FFT variables
//   input_re - pointer to an array of the real part of the output,
//     size nfft/2 + 1
//   input_im - pointer to an array of the imaginary part of the output,
//     size nfft/2 + 1
//   output - pointer to an array of (real) input values, size nfft
void fft_inverse(fft_vars* membvars, float* input_re, float* input_im, float* output)
{
	int ti;
	int nfft;
	int hnfft;
	int numfreqs;

	nfft = membvars->nfft;
	hnfft = nfft/2;
	numfreqs = membvars->numfreqs;

	for (ti=0; ti<hnfft; ti++) {
		membvars->fft_data[ti] = input_re[ti];
		membvars->fft_data[nfft-1-ti] = input_im[ti+1];
	}
	membvars->fft_data[hnfft] = input_re[hnfft];

	mayer_realifft(nfft, membvars->fft_data);

	for (ti=0; ti<nfft; ti++) {
		output[ti] = membvars->fft_data[ti];
	}
}

// DONE WITH FFT CODE



void *autotune_class;

typedef struct _autotune 	// Data structure for this object PD
{
    t_object x_obj;
    float x_f;
		
	// parameters

	float fTune;
	float fFixed;
	float fPull;
	float fA;
	float fBb;
	float fB;
	float fC;
	float fDb;
	float fD;
	float fEb;
	float fE;
	float fF;
	float fGb;
	float fG;
	float fAb;
	float fAmount;
	//float fGlide;
	float fSmooth;
	float fScwarp;
	float fLfoamp;
	float fShift;
	float fLforate;
	float fLfoshape;
	float fLfosymm;
	float fLfoquant;
	float fFcorr;
	float fFwarp;
	float fMix;
	float fPitch;
	float fConf;

	fft_vars* fx; 		// member variables for fft routine
	unsigned long fs; 			// Sample rate
	unsigned long cbsize; 		// size of circular buffer
	unsigned long corrsize; 	// cbsize/2 + 1
	unsigned long cbiwr;
	unsigned long cbord;
	float* cbi; // circular input buffer
	float* cbf; // circular formant correction buffer
	float* cbo; // circular output buffer
	float* cbwindow; // hann of length N/2, zeros for the rest
	float* acwinv; // inverse of autocorrelation of window
	float* hannwindow; // length-N hann
	int noverlap;
	float* ffttime;
	float* fftfreqre;
	float* fftfreqim;
	
	// VARIABLES FOR LOW-RATE SECTION
	float aref; 				// A tuning reference (Hz)
	float inpitch; 				// Input pitch (semitones)
	float conf; 				// Confidence of pitch period estimate (between 0 and 1)
	float outpitch; 			// Output pitch (semitones)
	float vthresh; 				// Voiced speech threshold
	float pperiod; 				// Pitch period (seconds)
	float pitch; 				// Pitch (semitones)
	float pitchpers; 			// Pitch persist (semitones)
	float pmax; 				// Maximum allowable pitch period (seconds)
	float pmin; 				// Minimum allowable pitch period (seconds)
	unsigned long nmax; 		// Maximum period index for pitch prd est
	unsigned long nmin; 		// Minimum period index for pitch prd est
	float lrshift; 				// Shift prescribed by low-rate section
	int ptarget; 				// Pitch target, between 0 and 11
  	float sptarget; 			// Smoothed pitch target
	//float sptarget; 			// Smoothed pitch target
	//int wasvoiced; 				// 1 if previous frame was voiced
	//float persistamt; 			// Proportion of previous pitch considered during next voiced period
	//float glidepersist;
	float lfophase;
	
	// VARIABLES FOR PITCH SHIFTER
	//float phprd; 				// phase period
	float phprdd;				// default (unvoiced) phase period
	float inphinc; 				// input phase increment
	float outphinc; 			// input phase increment
	float phincfact; 			// factor determining output phase increment
	float phasein;
	float phaseout;
	float* frag; 				// windowed fragment of speech
	unsigned long fragsize; 	// size of fragment in samples	
	float clockinterval;
	void *periodout, *confout; 	// floatout for pitch
	void *clock;

	// VARIABLES FOR FORMANT CORRECTOR
	unsigned int ford;
	float falph;
	float flamb;
	float* fk;
	float* fb;
	float* fc;
	float* frb;
	float* frc;
	float* fsig;
	float* fsmooth;
	float fhp;
	float flp;
	float flpa;
	float** fbuff;
	float* ftvec;
	float fmute;
	float fmutealph;
	
} t_autotune;	



//prototypes for methods
void *autotune_new(t_symbol *s, int argc, t_atom *argv);
void autotune_free(t_autotune *x);	
void autotune_init(t_autotune *x, unsigned long sr);
t_int *autotune_perform(t_int *w);
void autotune_dsp(t_autotune *x, t_signal **sp, short *count);
void autotune_assist(t_autotune *x, void *b, long m, long a, char *s);
//void autotune_setparam(t_autotune *x, t_symbol *m, short argc, t_atom *argv);
void autotune_list(t_autotune *x, t_symbol *m, short argc, t_atom *argv);
void autotune_mix(t_autotune *x, t_floatarg f);
void autotune_shift(t_autotune *x, t_floatarg f);
void autotune_pull(t_autotune *x, t_floatarg f);
void autotune_pull_pitch(t_autotune *x, t_floatarg f);
void autotune_tune(t_autotune *x, t_floatarg f);
void autotune_correction(t_autotune *x, t_floatarg f);
void autotune_smooth(t_autotune *x, t_floatarg f);
void autotune_lfo_depth(t_autotune *x, t_floatarg f);
void autotune_lfo_rate(t_autotune *x, t_floatarg f);
void autotune_lfo_shape(t_autotune *x, t_floatarg f);
void autotune_lfo_symmetry(t_autotune *x, t_floatarg f);
void autotune_formant_correction(t_autotune *x, t_floatarg f);
void autotune_formant_warp(t_autotune *x, t_floatarg f);
void autotune_scale_warp(t_autotune *x, t_floatarg f);
void autotune_confidence(t_autotune *x, t_floatarg f);
void autotune_processclock(t_autotune *x);




void autotune_tilde_setup(void)
{
	autotune_class = class_new(gensym("autotune~"), (t_newmethod)autotune_new, 
							   (t_method)autotune_free ,sizeof(t_autotune),0,A_GIMME,0);
	CLASS_MAINSIGNALIN(autotune_class, t_autotune, x_f );
	class_addmethod(autotune_class,(t_method)autotune_dsp, gensym("dsp"),0);
	class_addmethod(autotune_class,(t_method)autotune_assist, gensym("assist"),0);
	class_addmethod(autotune_class,(t_method)autotune_list,gensym("list"),A_GIMME,0);
	class_addmethod(autotune_class,(t_method)autotune_mix,gensym("mix"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_shift,gensym("shift"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_pull,gensym("pull"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_pull_pitch,gensym("pullpitch"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_tune,gensym("tune"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_correction,gensym("correct"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_smooth,gensym("smooth"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_lfo_depth,gensym("lfodepth"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_lfo_rate,gensym("lforate"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_lfo_shape,gensym("lfoshape"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_lfo_symmetry,gensym("lfosym"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_formant_correction,gensym("fcorr"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_formant_warp,gensym("warp"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_scale_warp,gensym("scwarp"),A_FLOAT,0);
	class_addmethod(autotune_class,(t_method)autotune_confidence,gensym("confidence"),A_FLOAT,0);
	post("autotune~ v.0.9.2");
	post("Ivica Ico Bukvic 2016");
}


int autotune_power_of_two (int x)
{
	 while (((x % 2) == 0) && x > 1) /* While x is even and > 1 */
	 	x /= 2;
	 return (x == 1);
}


// Create - Contruction of signal inlets and outlets
void *autotune_new(t_symbol *s, int argc, t_atom *argv)
{
    unsigned long sr;
    int fft_window = 0;
    int fft_hop = 0;

    t_autotune *x = (t_autotune *)pd_new(autotune_class);

    if (argc && argv->a_type == A_FLOAT) // got FFT window
    {
        fft_window = (int)atom_getint(argv);

        if (fft_window < 0) fft_window = 0;
        if (!autotune_power_of_two(fft_window)) fft_window = 0;
        argv++;
        argc--;
    }

    if (argc && argv->a_type == A_FLOAT) // got hop number
    {
        fft_hop = (int)atom_getint(argv);

        if (fft_hop < 2 || fft_hop > 32) fft_hop = 0;
        argv++;
        argc--;
    }

    if (fft_window != 0)
    	x->cbsize = fft_window;

    if (fft_hop != 0)
    	x->noverlap = fft_hop;
	
	if(sys_getsr()) sr = sys_getsr();
	else sr = 44100;
	
	autotune_init(x , sr);
	
	// second argument = number of signal inlets
    //dsp_setup((t_object *)x, 1);
	//inlet_new(&x->x_obj, &x->x_obj.ob_pd,gensym("signal"));
	//inlet_new(&x->x_obj, &x->x_obj.ob_pd,gensym("signal"), gensym("signal"));
	
	//floatinlet_new (&x->x_obj, &x->fAmount);
	//inlet_new (&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("correct"));

	//floatinlet_new (&x->x_obj, &x->fGlide);

	//floatinlet_new (&x->x_obj, &x->fSmooth);
	//inlet_new (&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("smooth"));

	//floatinlet_new (&x->x_obj, &x->fMix);
	//inlet_new (&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("mix"));

	//floatinlet_new (&x->x_obj, &x->fShift);
	//inlet_new (&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("shift"));

	//floatinlet_new (&x->x_obj, &x->fTune);
	//inlet_new (&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("tune"));

	//floatinlet_new (&x->x_obj, &x->fPersist);
	//symbolinlet_new (&x->x_obj, &x->fAmount);
	
	outlet_new(&x->x_obj, gensym("signal"));
	//x->f_out = outlet_new(&x->x_obj, &s_float);
	x->confout = outlet_new(&x->x_obj, &s_float);
	x->periodout = outlet_new(&x->x_obj, &s_float);
	//x->confout = outlet_new(x, "float");
	//x->periodout = outlet_new(x, "float");
    x->clock = clock_new(x,(t_method)autotune_processclock);
	x->clockinterval = 10.;
    
    return (x);									// Return the pointer
}

// Destroy
void autotune_free(t_autotune *x)
{
	//dsp_free((t_object *)x);		// Always call dsp_free first in this routine
	unsigned int ti;
	clock_unset(x->clock);
	clock_free(x->clock); 
	fft_des(x->fx);
	freebytes(x->cbi,0);
	freebytes(x->cbf,0);
	freebytes(x->cbo,0);
	//freebytes(x->cbonorm,0);
	freebytes(x->cbwindow,0);
	freebytes(x->hannwindow,0);
	freebytes(x->acwinv,0);
	freebytes(x->frag,0);
	freebytes(x->ffttime,0);
	freebytes(x->fftfreqre,0);
	freebytes(x->fftfreqim,0);
	free(x->fk);
	free(x->fb);
	free(x->fc);
	free(x->frb);
	free(x->frc);
	free(x->fsmooth);
	free(x->fsig);
	for (ti=0; ti<x->ford; ti++) {
		free(x->fbuff[ti]);
	}
	free(x->fbuff);
	free(x->ftvec);
}

void autotune_init(t_autotune *x,unsigned long sr)
{
	unsigned long ti;

	x->fs = sr;
	x->aref = 440;
	x->fTune = x->aref;
	
	if (x->cbsize == 0)
	{
		if (x->fs >=88200) {
			x->cbsize = 4096;
		}
		else {
			x->cbsize = 2048;
		}
	}
	x->corrsize = x->cbsize / 2 + 1;
	
	x->pmax = 1/(float)70;  // max and min periods (ms)
	x->pmin = 1/(float)2400; // eventually may want to bring these out as sliders
	
	x->pperiod = x->pmax;
	
	x->nmax = (unsigned long)(x->fs * x->pmax);
	if (x->nmax > x->corrsize) {
		x->nmax = x->corrsize;
	}
	x->nmin = (unsigned long)(x->fs * x->pmin);
	
	x->cbi = (float*) calloc(x->cbsize, sizeof(float));
	x->cbf = (float*) calloc(x->cbsize, sizeof(float));
	x->cbo = (float*) calloc(x->cbsize, sizeof(float));
	//x->cbonorm = (float*) calloc(x->cbsize, sizeof(float));
	
	x->cbiwr = 0;
	x->cbord = 0;

	x->lfophase = 0;

	// Initialize formant corrector
	x->ford = 7; // should be sufficient to capture formants
	x->falph = pow(0.001, (float) 80 / (x->fs));
	x->flamb = -(0.8517*sqrt(atan(0.06583*x->fs))-0.1916); // or about -0.88 @ 44.1kHz
	x->fk = calloc(x->ford, sizeof(float));
	x->fb = calloc(x->ford, sizeof(float));
	x->fc = calloc(x->ford, sizeof(float));
	x->frb = calloc(x->ford, sizeof(float));
	x->frc = calloc(x->ford, sizeof(float));
	x->fsig = calloc(x->ford, sizeof(float));
	x->fsmooth = calloc(x->ford, sizeof(float));
	x->fhp = 0;
	x->flp = 0;
	x->flpa = pow(0.001, (float) 10 / (x->fs));
	x->fbuff = (float**) malloc((x->ford)*sizeof(float*));
	for (ti=0; ti<x->ford; ti++) {
		x->fbuff[ti] = calloc(x->cbsize, sizeof(float));
	}
	x->ftvec = calloc(x->ford, sizeof(float));
	x->fmute = 1;
	x->fmutealph = pow(0.001, (float)1 / (x->fs));
	
	// Standard raised cosine window, max height at N/2
	x->hannwindow = (float*) calloc(x->cbsize, sizeof(float));
	for (ti=0; ti<x->cbsize; ti++) {
		x->hannwindow[ti] = -0.5*cos(2*PI*ti/(x->cbsize - 1)) + 0.5;
	}
	
	// Generate a window with a single raised cosine from N/4 to 3N/4
	x->cbwindow = (float*) calloc(x->cbsize, sizeof(float));
	for (ti=0; ti<(x->cbsize / 2); ti++) {
		x->cbwindow[ti+x->cbsize/4] = -0.5*cos(4*PI*ti/(x->cbsize - 1)) + 0.5;
	}
	
	if (x->noverlap == 0)
		x->noverlap = 4;

	//fprintf(stderr,"%d %d\n", x->cbsize, x->noverlap);
	
	x->fx = fft_con(x->cbsize);
	
	x->ffttime = (float*) calloc(x->cbsize, sizeof(float));
	x->fftfreqre = (float*) calloc(x->corrsize, sizeof(float));
	x->fftfreqim = (float*) calloc(x->corrsize, sizeof(float));
	
	
	// ---- Calculate autocorrelation of window ----
	x->acwinv = (float*) calloc(x->cbsize, sizeof(float));
	for (ti=0; ti<x->cbsize; ti++) {
		x->ffttime[ti] = x->cbwindow[ti];
	}
	fft_forward(x->fx, x->cbwindow, x->fftfreqre, x->fftfreqim);
	for (ti=0; ti<x->corrsize; ti++) {
		x->fftfreqre[ti] = (x->fftfreqre[ti])*(x->fftfreqre[ti]) + (x->fftfreqim[ti])*(x->fftfreqim[ti]);
		x->fftfreqim[ti] = 0;
	}
	fft_inverse(x->fx, x->fftfreqre, x->fftfreqim, x->ffttime);
	for (ti=1; ti<x->cbsize; ti++) {
		x->acwinv[ti] = x->ffttime[ti]/x->ffttime[0];
		if (x->acwinv[ti] > 0.000001) {
			x->acwinv[ti] = (float)1/x->acwinv[ti];
		}
		else {
			x->acwinv[ti] = 0;
		}
	}
	x->acwinv[0] = 1;
	// ---- END Calculate autocorrelation of window ----	
	
	x->lrshift = 0;
	x->ptarget = 0;
	x->sptarget = 0;
	//x->sptarget = 0;
	//x->wasvoiced = 0;
	//x->persistamt = 0;
	
	//x->glidepersist = 100; // 100 ms glide persist
	
	x->vthresh = 0.7;  //  The voiced confidence (unbiased peak) threshold level
	
	// Pitch shifter initialization
	x->phprdd = 0.01; // Default period
	//x->phprd = x->phprdd;
	x->inphinc = (float)1/(x->phprdd * x->fs);
	x->phincfact = 1;
	x->phasein = 0;
	x->phaseout = 0;
	x->frag = (float*) calloc(x->cbsize, sizeof(float));
	x->fragsize = 0;
}

//*********************//
// DSP Methods PD //
//*********************//

void autotune_dsp(t_autotune *x, t_signal **sp, short *count)
{
	clock_delay(x->clock, 0.);
	
	if(x->fs != sp[0]->s_sr)  autotune_init(x, sp[0]->s_sr);
	
	dsp_add(autotune_perform, 4, x, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
}


//*************************//
// Perform Routine PD//
//*************************//
t_int *autotune_perform(t_int *w)
{
	t_autotune *x = (t_autotune *)(w[1]); // object is first arg 
	t_float *in = (t_float *)(w[2]);
	t_float *out = (t_float *)(w[3]);
	unsigned long SampleCount = (unsigned long)(w[4]);
	
	// copy struct variables to local
	
	/*float fA = x->fA;
	float fBb = x->fBb;
	float fB = x->fB;
	float fC = x->fC;
	float fDb = x->fDb;
	float fD = x->fD;
	float fEb = x->fEb;
	float fE = x->fE;
	float fF = x->fF;
	float fGb = x->fGb;
	float fG = x->fG;
	float fAb = x->fAb;*/
	//float fGlide = x->fGlide;
	//float fPersist = x->glidepersist;
	
	int iNotes[12];
	int iPitch2Note[12];
	int iNote2Pitch[12];
	int numNotes;

	float fAmount = x->fAmount;
	float fSmooth = x->fSmooth * 0.8;
	float fTune = x->fTune;
	iNotes[0] = (int) x->fA;
	iNotes[1] = (int) x->fBb;
	iNotes[2] = (int) x->fB;
	iNotes[3] = (int) x->fC;
	iNotes[4] = (int) x->fDb;
	iNotes[5] = (int) x->fD;
	iNotes[6] = (int) x->fEb;
	iNotes[7] = (int) x->fE;
	iNotes[8] = (int) x->fF;
	iNotes[9] = (int) x->fGb;
	iNotes[10] = (int) x->fG;
	iNotes[11] = (int) x->fAb;
	float fFixed = x->fFixed;
	float fPull = x->fPull;
	float fShift = x->fShift;
	int iScwarp = x->fScwarp;
	float fLfoamp = x->fLfoamp;
	float fLforate = x->fLforate;
	float fLfoshape = x->fLfoshape;
	float fLfosymm = x->fLfosymm;
	int iLfoquant = x->fLfoquant;
	int iFcorr = x->fFcorr;
	float fFwarp = x->fFwarp;
	float fMix = x->fMix;
	
	//x->aref = (float)440*pow(2,fTune/12);
	
	unsigned long int lSampleIndex;
	
	unsigned long N = x->cbsize;
	unsigned long Nf = x->corrsize;
	unsigned long fs = x->fs;

	float pmax = x->pmax;
	float pmin = x->pmin;
	unsigned long nmax = x->nmax;
	unsigned long nmin = x->nmin;
	
	//float pperiod = x->pperiod;
	//float pitch = x->pitch;
	
		//
	
	volatile long int ti;
	volatile long int ti2;
	volatile long int ti3;
	volatile long int ti4;
	volatile float tf;
	volatile float tf2;
	volatile float tf3;

	// Variables for cubic spline interpolator
	volatile float indd;
	volatile int ind0;
	volatile int ind1;
	volatile int ind2;
	volatile int ind3;
	volatile float vald;
	volatile float val0;
	volatile float val1;
	volatile float val2;
	volatile float val3;

	volatile int lowersnap;
	volatile int uppersnap;
	volatile float lfoval;

	volatile float pperiod;
	volatile float inpitch;
	volatile float conf;
	volatile float outpitch;
	volatile float aref;
	volatile float fa;
	volatile float fb;
	volatile float fc;
	volatile float fk;
	volatile float flamb;
	volatile float frlamb;
	volatile float falph;
	volatile float foma;
	volatile float f1resp;
	volatile float f0resp;
	volatile float flpa;
	volatile int ford;

	// Some logic for the semitone->scale and scale->semitone conversion
	// If no notes are selected as being in the scale, instead snap to all notes
	ti2 = 0;
	for (ti=0; ti<12; ti++) {
		if (iNotes[ti]>=0) {
			iPitch2Note[ti] = ti2;
			iNote2Pitch[ti2] = ti;
			ti2 = ti2 + 1;
		}
		else {
			iPitch2Note[ti] = -1;
		}
	}
	numNotes = ti2;
	while (ti2<12) {
		iNote2Pitch[ti2] = -1;
		ti2 = ti2 + 1;
	}
	if (numNotes==0) {
		for (ti=0; ti<12; ti++) {
			iNotes[ti] = 1;
			iPitch2Note[ti] = ti;
			iNote2Pitch[ti] = ti;
		}
		numNotes = 12;
	}
	iScwarp = (iScwarp + numNotes*5)%numNotes;

	ford = x->ford;
	falph = x->falph;
	foma = (float)1 - falph;
	flpa = x->flpa;
	flamb = x->flamb;
	tf = pow((float)2,fFwarp/2)*(1+flamb)/(1-flamb);
	frlamb = (tf - 1)/(tf + 1);

	x->aref = (float)fTune;

	N = x->cbsize;
	Nf = x->corrsize;
	fs = x->fs;

	pmax = x->pmax;
	pmin = x->pmin;
	nmax = x->nmax;
	nmin = x->nmin;

	aref = x->aref;
	pperiod = x->pmax;
	inpitch = x->inpitch;
	conf = x->conf;
	outpitch = x->outpitch;

	
	//******************//
	//  MAIN DSP LOOP   //
	//******************//
	for (lSampleIndex = 0; lSampleIndex < SampleCount; lSampleIndex++)  
	{
		
		// load data into circular buffer
		tf = (float) *(in++);
		ti4 = x->cbiwr;
		//fprintf(stderr,"ti4=%d N=%d\n", ti4, N);
		x->cbi[ti4] = tf;
		/*x->cbiwr++;
		if (x->cbiwr >= N) {
			x->cbiwr = 0;
		}*/
		
		if (iFcorr>=1) {
			// Somewhat experimental formant corrector
			//  formants are removed using an adaptive pre-filter and
			//  re-introduced after pitch manipulation using post-filter
			// tf is signal input
			fa = tf - x->fhp; // highpass pre-emphasis filter
			x->fhp = tf;
			fb = fa;
			for (ti=0; ti<(long)ford; ti++) {
				x->fsig[ti] = fa*fa*foma + x->fsig[ti]*falph;
				fc = (fb-x->fc[ti])*flamb + x->fb[ti];
				x->fc[ti] = fc;
				x->fb[ti] = fb;
				fk = fa*fc*foma + x->fk[ti]*falph;
				x->fk[ti] = fk;
				tf = fk/(x->fsig[ti] + 0.000001);
				tf = tf*foma + x->fsmooth[ti]*falph;
				x->fsmooth[ti] = tf;
				x->fbuff[ti][ti4] = tf;
				fb = fc - tf*fa;
				fa = fa - tf*fc;
			}
			x->cbf[ti4] = fa;
			// Now hopefully the formants are reduced
			// More formant correction code at the end of the DSP loop
		}
		else {
			x->cbf[ti4] = tf;
		}

		//fprintf(stderr,"x->cbf[ti4]=%f\n", x->cbf[ti4]);

	    // Input write pointer logic
	    x->cbiwr++;
	    if (x->cbiwr >= N) {
	      x->cbiwr = 0;
	    }

		
		// ********************//
		// * Low-rate section *//
		// ********************//

		//fprintf(stderr,"overlap=%d outpitch=%f inpitch=%f\n", (x->cbiwr)%(N/x->noverlap), outpitch, inpitch);
		//fprintf(stderr,"outpitch=%f inpitch=%f\n", outpitch, inpitch);
		
		// Every N/noverlap samples, run pitch estimation / correction code
		if ((x->cbiwr)%(N/x->noverlap) == 0) {
			
			//fprintf(stderr,"ti4=%d N=%d\n", ti4, N);
			// ---- Obtain autocovariance ---- //
			
			// Window and fill FFT buffer
			ti2 = (long) x->cbiwr;
			for (ti=0; ti<(long)N; ti++) {
				x->ffttime[ti] = (float)(x->cbi[(ti2-ti)%N]*x->cbwindow[ti]);
			}
			
			// Calculate FFT
			fft_forward(x->fx, x->ffttime, x->fftfreqre, x->fftfreqim);
			
			// Remove DC
			x->fftfreqre[0] = 0;
			x->fftfreqim[0] = 0;
			
			// Take magnitude squared
			for (ti=1; ti< (long) Nf; ti++) {
				x->fftfreqre[ti] = (x->fftfreqre[ti])*(x->fftfreqre[ti]) + (x->fftfreqim[ti])*(x->fftfreqim[ti]);
				x->fftfreqim[ti] = 0;
			}
			
			// Calculate IFFT
			fft_inverse(x->fx, x->fftfreqre, x->fftfreqim, x->ffttime);
			
			// Normalize
			for (ti=1; ti<(long)N; ti++) {
				x->ffttime[ti] = x->ffttime[ti] / x->ffttime[0];
			}
			x->ffttime[0] = 1;
			
			//  ---- END Obtain autocovariance ----
			
			
			//  ---- Calculate pitch and confidence ----
			
			// Calculate pitch period
			//   Pitch period is determined by the location of the max (biased)
			//   peak within a given range
			//   Confidence is determined by the corresponding unbiased height
			tf2 = 0;
			pperiod = pmin;
			for (ti=nmin; ti<(long)nmax; ti++) {
				ti2 = ti-1;
				ti3 = ti+1;
				if (ti2<0) {
					ti2 = 0;
				}
				if (ti3>(long)Nf) {
					ti3 = Nf;
				}
				tf = x->ffttime[ti];
				
				if (tf>x->ffttime[ti2] && tf>=x->ffttime[ti3] && tf>tf2) {
					tf2 = tf;
					ti4 = ti;
					//conf = tf*x->acwinv[ti];
					//pperiod = (float)ti/fs;
				}
			}
			if (tf2>0) {
				conf = tf2*x->acwinv[ti4];
				if (ti4>0 && ti4<(long)Nf) {
					// Find the center of mass in the vicinity of the detected peak
					tf = x->ffttime[ti4-1]*(ti4-1);
					tf = tf + x->ffttime[ti4]*(ti4);
					tf = tf + x->ffttime[ti4+1]*(ti4+1);
					tf = tf/(x->ffttime[ti4-1] + x->ffttime[ti4] + x->ffttime[ti4+1]);
					pperiod = tf/fs;
				}
				else {
					pperiod = (float)ti4/fs;
				}
			}
			
			// Convert to semitones
			tf = (float) -12*log10((float)aref*pperiod)*L2SC;
			//fprintf(stderr,"tf=%f aref=%f pperiod=%f\n", tf, aref, pperiod);
			//post("pperiod=%f conf=%f\n", pperiod, conf);
			float pp_test = x->pperiod/(x->pperiod - pperiod);
			if (pp_test < 0.5 || pp_test > 2)
				pp_test = 1;
			else
				pp_test = 0;
			if (conf>=x->vthresh && tf == tf) { // second check is for NANs
				inpitch = tf;
				x->inpitch = tf; // update pitch only if voiced
				x->pperiod = pperiod;
			}
			x->conf = conf;

			x->fPitch = inpitch;
			x->fConf = conf;

			//x->pitch = pitch;
			//x->pperiod = pperiod;
			//x->conf = conf;
			
			//  ---- END Calculate pitch and confidence ----
			
			
			/*
			//  ---- Determine pitch target ----
			
			// If voiced
			if (conf>=x->vthresh) {
				// TODO: Scale sliders
				// Determine pitch target
				tf = -1;
				tf2 = 0;
				tf3 = 0;
				for (ti=0; ti<12; ti++) {
					switch (ti) {
						case 0:
							tf2 = fA;
							break;
						case 1:
							tf2 = fBb;
							break;
						case 2:
							tf2 = fB;
							break;
						case 3:
							tf2 = fC;
							break;
						case 4:
							tf2 = fDb;
							break;
						case 5:
							tf2 = fD;
							break;
						case 6:
							tf2 = fEb;
							break;
						case 7:
							tf2 = fE;
							break;
						case 8:
							tf2 = fF;
							break;
						case 9:
							tf2 = fGb;
							break;
						case 10:
							tf2 = fG;
							break;
						case 11:
							tf2 = fAb;
							break;
					}
					// 	  if (ti==x->ptarget) {
					// 	    tf2 = tf2 + 0.01; // add a little hysteresis
					// 	  }
					tf2 = tf2 - (float)fabs( (pitch-(float)ti)/6 - 2*floorf(((pitch-(float)ti)/12 + 0.5)) ); // like a likelihood function
					if (tf2>=tf) {                                                                           // that we're maximizing
						tf3 = (float)ti;                                                                       // to find the target pitch
						tf = tf2;
					}
				}
				x->ptarget = tf3;
				
				// Glide persist
				if (x->wasvoiced == 0) {
					x->wasvoiced = 1;
					tf = x->persistamt;
					x->sptarget = (1-tf)*x->ptarget + tf*x->sptarget;
					x->persistamt = 1;
				}
				
				// Glide on circular scale