/* digitar.c
 * Fri Jun 18 17:37:07 PDT 2010 Kevin Karplus
 *
 * Crude implementation of the digitar algorithm
 * using the version in Alex Strong's patent, rather than the original paper.
 */
 
#include <stdlib.h>	/* for random() */
#include <assert.h>
#include <stdio.h>

#include "digitar.h"

#define RAND_SCALE 	(1./2147483648.)

inline double drandom(void)
{   return random()* RAND_SCALE;
}

struct wavetable* alloc_table(int max_length)
{
    struct wavetable* wav;
    wav = calloc(1, sizeof(struct wavetable));
    assert(wav);
    wav->table = calloc(max_length, sizeof(short int));
    assert(wav->table);
    wav->alloc_length = max_length;
    /* note: length and sum_table start at 0 */
    return wav;
}

/* free_table(&table) frees the table and sets table=NULL */
void free_table(struct wavetable** tab_p)
{    free((*tab_p)->table);
    free(*tab_p);
    *tab_p=0;
}



/* return a sample from a wave table of length tab_length
 * phase should be in the range [0..tab_length) and is updated
 * frequency is not in Hz, but in steps/sample.
 * freq = f_Hz / SAMPLE_RATE * tab_length
 */
inline short int wave_osc(struct wavetable* wav,
	fixed_point freq)
{
    register fixed_point p = wav->phase;
    p += freq;
    if (p >= (wav->length)<<FIXED_POINT_SHIFT) p -= (wav->length)<<FIXED_POINT_SHIFT;
    wav->phase = p;
    return wav->table[p>>FIXED_POINT_SHIFT];
}

inline void fill_table(struct wavetable *wav, int amplitude)
{
    int i;	/* index of element to fill */
    int x;	/* value of element */
    int sum=0;	/* sum of table */
    
    int num_hi_needed = wav->length/2;
    
    for (i=wav->length-1; i>=0; i--)
    {	x = amplitude;
        if (drandom() * i < num_hi_needed) 	{    num_hi_needed--;	}
	else {   x = 0- x; }
    	wav->table[i] = x;
	sum+=x;
    }
    wav->sum_table = sum;
    fprintf(stderr,"DEBUG: num_hi_needed=%d, sum=%d\n", num_hi_needed, sum);
}

inline short int update_table(struct wavetable *wav, int i, short int x)
{
    wav->sum_table -= wav->table[i];
    wav->sum_table += x;
    wav->table[i] = x;
    return x;
}

/* decay_one is done as a macro, rather than an in-line
 * because the 3 subscripts all need to be updated, and
 * I'd like them to remain in registers if possible.
 * I don't trust gcc's optimization to be able to do this correctly
 * with int * arguments.
 */
#define decay_one(wav, at_minus,at,at_plus)	\
	{   update_table(wav,at, 		\
		((wav->table[at_minus] +wav->table[at_plus])>>1)	\
		+(wav->sum_table>0? \
		   (wav->table[at]*wav->length>wav->sum_table? -4: 0): \
		   (wav->table[at]*wav->length<wav->sum_table? 4: 0)) ); \
	    at_plus = at;    at = at_minus;    at_minus--;	\
	    if (at_minus<0) at_minus=wav->length-1;		\
	}	    


/* Fill wav with +-amplitude  (in range 0..32767)
 * then decay for decay_passes full passes
 */
void pluck(struct wavetable *wav, int amplitude, int decay_passes)
{    
    int i;	/* counter for initial decay */
    int at_minus=0, at=1, at_plus=2;
    assert(wav->length >2);
    fill_table(wav,amplitude);
    for (i=wav->length*decay_passes; i>0 ; i--)
    {   decay_one(wav, at_minus,at,at_plus);
    }
}


/* play num_samples long note at frequency f_Hz
 * adding to array out, using wavetable  wav
 * Decay rate is given in decays/sample (0 <= decay_rate )
 */
short int* decay_note(struct wavetable*wav,
	float f_Hz,
	float decay_rate,
	int num_samples,
	short int *out)
{
    int i;
    
    /* phase oscillator for playing note */
    fixed_point freq_fp = (int)(f_Hz*wav->length/SAMPLE_RATE*FIXED_POINT_SCALING);
    
    /* decay phase for timing and decay pointers */
    fixed_point decay_phase=0;
    fixed_point decay_fp = (int)(decay_rate*FIXED_POINT_SCALING);
    int at_minus=0, at = 1, at_plus=2;
    assert(wav->length >2);
    
    fprintf(stderr,"DEBUG: f_Hz=%f freq_fp=%d length=%d sum_table=%d\n",
		    f_Hz, freq_fp, wav->length, wav->sum_table);
    
    for (i=0; i<num_samples; i++)
    {   out[i] += wave_osc(wav, freq_fp);
        decay_phase += decay_fp;
	while (decay_phase>=FIXED_POINT_SCALING)
	{   decay_one(wav, at_minus,at,at_plus);
	    decay_phase-=FIXED_POINT_SCALING;
	}
    }
    fprintf(stderr,"DEBUG: sum_table=%d\n", wav->sum_table);
    return out+num_samples;
}


/* play num_samples long note at frequency f_Hz
 * adding into array out, using wavetable  wav
 * Decay rate is given in decays/sample (0 <= decay_rate )
 *
 * probability of bow slipping in (delta_v/v)^2/force
 *	(clipped at 1, of course)
 * bow_velocity should be fairly small (100? 1000?)
 * bow_force should be fairly small (10.?)
 */
void bow_note(struct wavetable*wav,
	float f_Hz,
	float decay_rate,
	float bow_force,
	int bow_velocity,
	int num_samples,
	short int **out_p)
{
    int i;
    short int *out = *out_p;

    /* phase oscillator for playing note */
    fixed_point freq_fp = (int)(f_Hz*wav->length/SAMPLE_RATE*FIXED_POINT_SCALING);
    
    /* decay phase for timing and decay pointers */
    fixed_point decay_phase=0;
    fixed_point decay_fp = (int)(decay_rate*FIXED_POINT_SCALING);
    int at_minus=0, at = 1, at_plus=2;
    
    /* bowing variables */
    short int x;	/* current position of string */
    short int old_x;	/* previous position of string */
    int bow_x;	/* modified position of string (bowing) */
    int new_x;	/* modified position of string (decaying) */
    int delta_v;	/* velocity(string)-velocity(bow) */
    float prob_scale_was_slip = 1./(bow_force*bow_velocity*bow_velocity);
    float prob_scale_was_stick = 0.1*prob_scale_was_slip;
    char was_sticking=1;
    
    float prob_slip;	/* probability of bow slipping */
    
    assert(wav->length >2);
    
    fprintf(stderr,"DEBUG: f_Hz=%f freq_fp=%d\n", f_Hz, freq_fp);
    
    for (i=0; i<num_samples; i++)
    {   out[i] += wave_osc(wav, freq_fp);

        decay_phase += decay_fp;
	while (decay_phase>=FIXED_POINT_SCALING)
	{   x = wav->table[at];			
	    old_x = wav->table[at_minus];	
	    new_x = ((old_x + wav->table[at_plus])>>1)
	    	+ (wav->sum_table>0? -4: 4);
	    delta_v = (x-old_x - bow_velocity);	
	    prob_slip = (was_sticking?prob_scale_was_stick:prob_scale_was_slip)*delta_v*delta_v ;
	    was_sticking=0;
//	    fprintf(stderr,"DEBUG: prob_slip=%f\n", prob_slip);	
	    if (prob_slip<1.0 && drandom() > prob_slip)	
	    {   /* string sticks to bow  */	
	        bow_x = old_x +bow_velocity;	
//		fprintf(stderr,"DEBUG: old_x=%d, x=%d, new_x=%d\n", old_x,x, new_x);	
	        if (-32768 <= bow_x && bow_x <= 32767)
	    	{   /* no overflow */		
	    	    new_x = bow_x;
		    was_sticking=1;
		}				
	    }
	    update_table(wav,at,new_x);
	    at_plus = at;    at = at_minus;    at_minus--;	
	    if (at_minus<0) at_minus=wav->length-1;		
	    decay_phase-=FIXED_POINT_SCALING;
	}
    }
    *out_p = out+num_samples;

}