/*

dryUp~ - a resonance attenuation external for Pd

Copyright 2009 William Brent

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 3 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, see <http://www.gnu.org/licenses/>.

*/

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

static t_class *dryUp_tilde_class;

typedef struct _dryUp_tilde 
{
    t_object x_obj;
    float sr;
    float n;
    int dsp_tick;
    int buffer_limit;
    int overlap;
    int overlap_chunk;
    int window;
    int window_half;
    float amp_scalar;
    float mystery_scalar;
    float mystery_scalar_reciprocal;
    float fader;
    t_float *prev_re;
    t_float *prev_im;
    t_float *hann;

	t_sample *signal_buf;
	t_sample *signal_buf_windowed;
	t_float *signal_R;
	t_float *signal_I;
	
	t_float *normaled;
	t_float *real_chain;
	t_float *imaginary_chain;
	
	t_sample *nonoverlapped_output;
	t_sample *final_output;
	
    float x_f;
    
} dryUp_tilde;


/* ------------------------ dryUp~ -------------------------------- */

void dryUp_tilde_fader(dryUp_tilde *x, t_floatarg f)
{
	x->fader = f;
}

void dryUp_tilde_amp_scalar(dryUp_tilde *x, t_floatarg a)
{
	x->amp_scalar = a/(float)x->window;
	
	post("amp scalar: %f", x->amp_scalar);
}

void dryUp_tilde_mystery_scalar(dryUp_tilde *x, t_floatarg ms)
{
	if(ms > 0)
	{
		x->mystery_scalar = ms;
		x->mystery_scalar_reciprocal = 1.0/ms;
		post("mystery scalar: %f", x->mystery_scalar);
	}
	else
		post("mystery scalar must be greater than 0.");	
}

void dryUp_tilde_window(dryUp_tilde *x, t_floatarg w)
{
	int i, isPow2, window, window_half, overlap_chunk;
	
	isPow2 = (int)w && !( ((int)w-1) & (int)w );
	
	if( !isPow2 )
		error("requested window size is not a power of 2.");
	else
	{
		window = w;
		window_half = w/2.0;
		overlap_chunk = w/x->overlap;
		
		x->hann = (t_float *)t_resizebytes(x->hann, x->window * sizeof(t_float), window * sizeof(t_float));	
		x->prev_re = (t_float *)t_resizebytes(x->prev_re, x->window * sizeof(t_float), window * sizeof(t_float));
		x->prev_im = (t_float *)t_resizebytes(x->prev_im, x->window * sizeof(t_float), window * sizeof(t_float));
		x->signal_buf = (t_sample *)t_resizebytes(x->signal_buf, x->window * sizeof(t_sample), window * sizeof(t_sample));
		x->signal_buf_windowed = (t_sample *)t_resizebytes(x->signal_buf_windowed, x->window * sizeof(t_sample), window * sizeof(t_sample));
		x->signal_R = (t_float *)t_resizebytes(x->signal_R, (x->window_half+1) * sizeof(t_float), (window_half+1) * sizeof(t_float));
		x->signal_I = (t_float *)t_resizebytes(x->signal_I, (x->window_half+1) * sizeof(t_float), (window_half+1) * sizeof(t_float));
		x->normaled = (t_float *)t_resizebytes(x->normaled, (x->window_half+1) * sizeof(t_float), (window_half+1) * sizeof(t_float));
		x->real_chain = (t_float *)t_resizebytes(x->real_chain, (x->window_half+1) * sizeof(t_float), (window_half+1) * sizeof(t_float));
		x->imaginary_chain = (t_float *)t_resizebytes(x->imaginary_chain, (x->window_half+1) * sizeof(t_float), (window_half+1) * sizeof(t_float));
		x->nonoverlapped_output = (t_sample *)t_resizebytes(x->nonoverlapped_output, (x->window*x->overlap) * sizeof(t_sample), (window*x->overlap) * sizeof(t_sample));	
		x->final_output = (t_sample *)t_resizebytes(x->final_output, (x->overlap_chunk + x->n) * sizeof(t_sample), (overlap_chunk + x->n) * sizeof(t_sample));
		
		x->window = window;
		x->window_half = window_half;
		x->overlap_chunk = overlap_chunk;
		
		// initialize hann window
		for(i=0; i<x->window; i++)
			x->hann[i] = 0.5 * (1 - cos(2*M_PI*i/x->window));
	
		for(i=0; i<x->window; i++)
		{
			x->signal_buf[i] = 0.0;
			x->signal_buf_windowed[i] = 0.0;
			x->prev_re[i] = 0.0;
			x->prev_im[i] = 0.0;
		};
	
		// init output buffer
		for(i=0; i<(x->overlap_chunk + x->n); i++)
			x->final_output[i] = 0.0;
				
		// initialize nonoverlapped_output
		for(i=0; i<(x->window*x->overlap); i++)
			x->nonoverlapped_output[i] = 0.0;
		
		// init (window_half+1 size arrays)
		for(i=0; i<=x->window_half; i++)
		{
			x->signal_R[i] = 0.0;
			x->signal_I[i] = 0.0;
			x->normaled[i] = 0.0;
			x->real_chain[i] = 0.0;
			x->imaginary_chain[i] = 0.0;
		};	
		
		x->amp_scalar = (2.0/3.0)/window;

		x->dsp_tick = 0; 
		x->buffer_limit = x->window/x->n/x->overlap;
	
		post("window size: %i, overlap: %i, amp scalar: %f", x->window, x->overlap, x->amp_scalar);
	}
}

void dryUp_tilde_overlap(dryUp_tilde *x, t_floatarg o)
{
	int i, isPow2, overlap_chunk;
	
	isPow2 = (int)o && !( ((int)o-1) & (int)o );
	
	if( !isPow2 )
		error("requested window size is not a power of 2.");
	else
	{
		overlap_chunk = x->window/x->overlap;
		
		x->nonoverlapped_output = (t_sample *)t_resizebytes(x->nonoverlapped_output, (x->window*x->overlap) * sizeof(t_sample), (x->window*o) * sizeof(t_sample));	
		x->final_output = (t_sample *)t_resizebytes(x->final_output, (x->overlap_chunk + x->n) * sizeof(t_sample), (overlap_chunk + x->n) * sizeof(t_sample));
		
		x->overlap_chunk = x->window/o;
		x->overlap = o;
		
		// init output buffer
		for(i=0; i<(x->overlap_chunk + x->n); i++)
			x->final_output[i] = 0.0;
				
		// initialize nonoverlapped_output
		for(i=0; i<(x->window*x->overlap); i++)
			x->nonoverlapped_output[i] = 0.0;
		
		x->dsp_tick = 0; 
		x->buffer_limit = x->window/x->n/x->overlap;
	
		post("window size: %i, overlap: %i, amp scalar: %f", x->window, x->overlap, x->amp_scalar);
	}
}

void *dryUp_tilde_new(t_floatarg window, t_floatarg overlap)
{
	int i, j, isPow2;
	
    dryUp_tilde *x = (dryUp_tilde *)pd_new(dryUp_tilde_class);
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("float"), gensym("fader"));
	outlet_new(&x->x_obj, gensym("signal"));

	if(window)
	{
		isPow2 = (int)window && !( ((int)window-1) & (int)window );
		if(isPow2)
			x->window = window;
		else
		{
			error("dryUp~: window size must be a power of 2.");
			x->window = 4096;
		}
	}
	else
		x->window = 4096;

	if(overlap)
	{
		isPow2 = (int)overlap && !( ((int)overlap-1) & (int)overlap );
		if(isPow2)
			x->overlap = overlap;
		else
		{
			error("dryUp~: overlap must be a power of 2.");
			x->overlap = 4;
		}
	}
	else
		x->overlap = 4;
		
	
	x->sr = 44100;
	x->n = 64;
	x->dsp_tick = 0; 
	x->window_half = x->window/2.0;
	x->overlap_chunk = x->window/x->overlap;
	x->amp_scalar = (2.0/3.0)/x->window;
	x->mystery_scalar = 100.0;
	x->mystery_scalar_reciprocal = 1.0/x->mystery_scalar;

	x->buffer_limit = x->window/x->n/x->overlap; // divide by overlap since we want to push out buffered audio to a window when the main buffer has been updated by 1/overlap the window sizse.
	
	x->fader = 0;

	x->hann = (t_float *)getbytes(0);
	x->hann = (t_float *)t_resizebytes(x->hann, 0, x->window * sizeof(t_float));
	
	x->prev_re = (t_float *)getbytes(0);
	x->prev_re = (t_float *)t_resizebytes(x->prev_re, 0, x->window * sizeof(t_float));
	
	x->prev_im = (t_float *)getbytes(0);
	x->prev_im = (t_float *)t_resizebytes(x->prev_im, 0, x->window * sizeof(t_float));

	x->signal_buf = (t_sample *)getbytes(0);
	x->signal_buf = (t_sample *)t_resizebytes(x->signal_buf, 0, x->window * sizeof(t_sample));

	x->signal_buf_windowed = (t_sample *)getbytes(0);
	x->signal_buf_windowed = (t_sample *)t_resizebytes(x->signal_buf_windowed, 0, x->window * sizeof(t_sample));
	
	x->signal_R = (t_float *)getbytes(0);
	x->signal_R = (t_float *)t_resizebytes(x->signal_R, 0, (x->window_half+1) * sizeof(t_float));

	x->signal_I = (t_float *)getbytes(0);
	x->signal_I = (t_float *)t_resizebytes(x->signal_I, 0, (x->window_half+1) * sizeof(t_float));

	x->normaled = (t_float *)getbytes(0);
	x->normaled = (t_float *)t_resizebytes(x->normaled, 0, (x->window_half+1) * sizeof(t_float));
	
	x->real_chain = (t_float *)getbytes(0);
	x->real_chain = (t_float *)t_resizebytes(x->real_chain, 0, (x->window_half+1) * sizeof(t_float));

	x->imaginary_chain = (t_float *)getbytes(0);
	x->imaginary_chain = (t_float *)t_resizebytes(x->imaginary_chain, 0, (x->window_half+1) * sizeof(t_float));
	
	x->nonoverlapped_output = (t_sample *)getbytes(0);
	x->nonoverlapped_output = (t_sample *)t_resizebytes(x->nonoverlapped_output, 0, (x->window*x->overlap) * sizeof(t_sample));
	
	x->final_output = (t_sample *)getbytes(0);
	x->final_output = (t_sample *)t_resizebytes(x->final_output, 0, (x->overlap_chunk + x->n) * sizeof(t_sample));
	
	
	
	// initialize hann window
 	for(i=0; i<x->window; i++)
 		x->hann[i] = 0.5 * (1 - cos(2*M_PI*i/x->window));

 	for(i=0; i<x->window; i++)
 	{
 		x->signal_buf[i] = 0.0;
 		x->signal_buf_windowed[i] = 0.0;
 		x->prev_re[i] = 0.0;
 		x->prev_im[i] = 0.0;
	};

	// init output buffer
	for(i=0; i<(x->overlap_chunk + x->n); i++)
		x->final_output[i] = 0.0;
 			
	// initialize nonoverlapped_output
 	for(i=0; i<(x->window*x->overlap); i++)
 		x->nonoverlapped_output[i] = 0.0;
	
	// init (window_half+1 size arrays)
	for(i=0; i<=x->window_half; i++)
	{
		x->signal_R[i] = 0.0;
		x->signal_I[i] = 0.0;
		x->normaled[i] = 0.0;
		x->real_chain[i] = 0.0;
		x->imaginary_chain[i] = 0.0;
	};	
	
    post("dryUp~: window size: %i, overlap: %i", x->window, x->overlap);
    
    return (x);
}


t_int *dryUp_tilde_perform(t_int *w)
{
    int i, j, k, n, window, window_half, overlap, overlap_chunk;
    float amp_scalar;

    dryUp_tilde *x = (dryUp_tilde *)(w[1]);

    t_sample *in = (t_float *)(w[2]);
    t_sample *out = (t_float *)(w[3]);
    
    n = w[4];
	window = x->window;
	window_half = x->window_half;
	overlap = x->overlap;
	overlap_chunk = x->overlap_chunk;
	
	amp_scalar = x->amp_scalar * x->mystery_scalar;  // mysterious 100 scalar	

	// shift previous contents back
	for(i=0; i<(window-n); i++)
		x->signal_buf[i] = x->signal_buf[n+i];
		
	// buffer most recent block
	for(i=0; i<n; i++)
		x->signal_buf[window-n+i] = in[i] * x->mystery_scalar_reciprocal;  // mysterious 0.01 scalar
		
	
	if(x->dsp_tick==x->buffer_limit)
	{	
		x->dsp_tick = 0;
		// window the signal
		for(i=0; i<window; i++)
			x->signal_buf_windowed[i] = x->signal_buf[i] * x->hann[i];
	
		// take FT of window
		mayer_realfft(window, x->signal_buf_windowed);
	
		// unpack mayer_realfft results into R&I arrays
		for(i=0; i<=window_half; i++)
		{
			x->signal_R[i] = x->signal_buf_windowed[i];
			if(fabs(x->signal_R[i]) < 0.0001)
				x->signal_R[i] = 0.0;
		}
		
		x->signal_I[0]=0;  // DC
		
		for(i=(window-1), j=1; i>window_half; i--, j++)
		{
			x->signal_I[j] = x->signal_buf_windowed[i];
			if(fabs(x->signal_I[j]) < 0.0001)
				x->signal_I[j] = 0.0;
		}
		
		x->signal_I[window_half]=0;  // Nyquist
		
		
		for(i=0; i<=window_half; i++)
		{	
			x->real_chain[i] = (x->signal_R[i] * x->prev_re[i]) + (x->signal_I[i] * x->prev_im[i]);
			x->imaginary_chain[i] = (x->signal_I[i] * x->prev_re[i]) - (x->signal_R[i] * x->prev_im[i]);
		
			// copy into data structure for next block
			x->prev_re[i] = x->signal_R[i];
			x->prev_im[i] = x->signal_I[i];
		};

		// go through this precession loop only once if overlap is 2, twice if it's 4...
		for(j=0; j<(int)(overlap*0.5); j++)
		{
			// only update real_chain & imaginary_chain if we go through the loop more than once
			if(j>0)
			{
				for(i=0; i<=window_half; i++)
				{
					x->real_chain[i] = x->signal_R[i];
					x->imaginary_chain[i] = x->signal_I[i];
				}
			}
			
			// signal_R and I are free now, so put the squared versions in there.
			for(i=0; i<=window_half; i++)
			{
				x->signal_R[i] = (x->real_chain[i] * x->real_chain[i]) - (x->imaginary_chain[i] * x->imaginary_chain[i]);
				x->signal_I[i] = x->real_chain[i] * x->imaginary_chain[i] * 2;
			};
		}
			
		for(i=0; i<=window_half; i++)
		{
			x->normaled[i] = q8_rsqrt((x->signal_R[i] * x->signal_R[i]) + (x->signal_I[i] * x->signal_I[i]) + 0.00000000000000000001);
	
			x->signal_R[i] *= x->normaled[i];
			x->signal_I[i] *= x->normaled[i];

			// temp storage in real_chain since signal_R is needed  in the second calculation below
			x->real_chain[i] = (x->prev_re[i] * x->signal_R[i]) - (x->prev_im[i] * x->signal_I[i]);
			
			x->signal_I[i] = (x->prev_re[i] * x->signal_I[i]) + (x->prev_im[i] * x->signal_R[i]);
			
			// copy from temp storage to signal_R
			x->signal_R[i] = x->real_chain[i];
		};
		
		// pack real and imaginary parts in correct order for mayer_realifft
		for(i=0; i<=window_half; i++)
			x->signal_buf_windowed[i] = x->signal_R[i];
		
		for(i=(window_half+1), j=(window_half-1); i<window; i++, j--)
			x->signal_buf_windowed[i] = x->signal_I[j];
			
		// resynth
		mayer_realifft(window, x->signal_buf_windowed);
			
			
		// window
		for(i=0; i<window; i++)
		{
			x->signal_buf_windowed[i] *= x->hann[i];
			x->signal_buf_windowed[i] *= amp_scalar;
		}
			
		// shift nonoverlapped blocks back
		for(i=0; i<(overlap-1); i++)
			for(j=0; j<window; j++)
				x->nonoverlapped_output[(i*window)+j] = x->nonoverlapped_output[((i+1)*window)+j];
		
		// write the new block
		for(i=0; i<window; i++)
			x->nonoverlapped_output[((overlap-1)*window)+i] = x->signal_buf_windowed[i];
		
		// shift final_out back
		for(i=0; i<n; i++)
			x->final_output[i] = x->final_output[overlap_chunk+i];
			
		// init this chunk of the final output so it can be summed in the for() below
		for(i=0; i<overlap_chunk; i++)
			x->final_output[n+i] = 0.0;
			
		// do the overlap/add
		for(i=0; i<overlap; i++)
			for(j=0; j<overlap_chunk; j++)
				x->final_output[n+j] += x->nonoverlapped_output[(i*window)+((overlap-i-1)*overlap_chunk)+j];
	};
	
	
	
	// output
	for(i=0; i<n; i++, out++)
		*out = in[i] - (x->final_output[((x->dsp_tick%x->buffer_limit)*n)+i] * x->fader);

	x->dsp_tick++;



    return (w+5);
}


void dryUp_tilde_dsp(dryUp_tilde *x, t_signal **sp)
{
	dsp_add(
		dryUp_tilde_perform,
		4,
		x,
		sp[0]->s_vec,
		sp[1]->s_vec,
		sp[0]->s_n
	); 


// compare n to stored n and recalculate hann & filterbank if different
	if( sp[0]->s_n != x->n  ||  sp[0]->s_sr != x->sr )
	{		
		x->final_output = (t_sample *)t_resizebytes(x->final_output, (x->overlap_chunk + x->n) * sizeof(t_sample), (x->overlap_chunk + (int)sp[0]->s_n) * sizeof(t_sample));

		x->sr = sp[0]->s_sr;
		x->n = sp[0]->s_n;
		
		x->buffer_limit = x->window/x->n/x->overlap;
	};

};

static void dryUp_tilde_free(dryUp_tilde *x)
{	
    t_freebytes(x->hann, x->window*sizeof(t_float));
    t_freebytes(x->prev_re, x->window*sizeof(t_float));
    t_freebytes(x->prev_im, x->window*sizeof(t_float));
    t_freebytes(x->signal_buf, x->window*sizeof(t_sample));
    t_freebytes(x->signal_buf_windowed, x->window*sizeof(t_sample));
    t_freebytes(x->signal_R, (x->window_half+1)*sizeof(t_float));
    t_freebytes(x->signal_I, (x->window_half+1)*sizeof(t_float));
    t_freebytes(x->normaled, (x->window_half+1)*sizeof(t_float));
    t_freebytes(x->real_chain, (x->window_half+1)*sizeof(t_float));
    t_freebytes(x->imaginary_chain, (x->window_half+1)*sizeof(t_float));
    t_freebytes(x->nonoverlapped_output, (x->window*x->overlap)*sizeof(t_sample));
    t_freebytes(x->final_output, (x->overlap_chunk + x->n)*sizeof(t_sample));
};

void dryUp_tilde_setup(void)
{
    dryUp_tilde_class = 
    class_new(
    	gensym("dryUp~"),
    	(t_newmethod)dryUp_tilde_new,
    	(t_method)dryUp_tilde_free,
        sizeof(dryUp_tilde),
        CLASS_DEFAULT,
        A_DEFFLOAT,
        A_DEFFLOAT,
		0
    );

    CLASS_MAINSIGNALIN(dryUp_tilde_class, dryUp_tilde, x_f);

	class_addmethod(
		dryUp_tilde_class, 
        (t_method)dryUp_tilde_fader,
		gensym("fader"),
		A_DEFFLOAT,
		0
	);

	class_addmethod(
		dryUp_tilde_class, 
        (t_method)dryUp_tilde_window,
		gensym("window"),
		A_DEFFLOAT,
		0
	);
	
	class_addmethod(
		dryUp_tilde_class, 
        (t_method)dryUp_tilde_overlap,
		gensym("overlap"),
		A_DEFFLOAT,
		0
	);

	class_addmethod(
		dryUp_tilde_class, 
        (t_method)dryUp_tilde_amp_scalar,
		gensym("amp_scalar"),
		A_DEFFLOAT,
		0
	);

	class_addmethod(
		dryUp_tilde_class, 
        (t_method)dryUp_tilde_mystery_scalar,
		gensym("mystery_scalar"),
		A_DEFFLOAT,
		0
	);
	
    class_addmethod(
    	dryUp_tilde_class,
    	(t_method)dryUp_tilde_dsp,
    	gensym("dsp"),
    	0
    );
}