/*
Copyright (C) 2004 Antoine Rousseau 
all material Copyright (c) 1997-1999 Miller Puckette.

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

This library 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
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA  

*/

/* a "sweet" saw generator, as described in pddoc/3.audio/K05.bandlimited .*/

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

#define UNITBIT32 1572864.  /* 3*2^19; bit 32 has place value 1 */

    /* machine-dependent definitions.  These ifdefs really
    should have been by CPU type and not by operating system! */
#ifdef IRIX
    /* big-endian.  Most significant byte is at low address in memory */
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#define int32 long  /* a data type that has 32 bits */
#else
#ifdef _WIN32
    /* little-endian; most significant byte is at highest address */
#define HIOFFSET 1
#define LOWOFFSET 0
#define int32 long
#else
#ifdef __FreeBSD__
#include <machine/endian.h>
#if BYTE_ORDER == LITTLE_ENDIAN
#define HIOFFSET 1
#define LOWOFFSET 0
#else
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#endif /* BYTE_ORDER */
#include <sys/types.h>
#define int32 int32_t
#endif
#ifdef __linux__

#include <endian.h>

#if !defined(__BYTE_ORDER) || !defined(__LITTLE_ENDIAN)                         
#error No byte order defined                                                    
#endif                                                                          
                                                                                
#if __BYTE_ORDER == __LITTLE_ENDIAN                                             
#define HIOFFSET 1                                                              
#define LOWOFFSET 0                                                             
#else                                                                           
#define HIOFFSET 0    /* word offset to find MSB */                             
#define LOWOFFSET 1    /* word offset to find LSB */                            
#endif /* __BYTE_ORDER */                                                       

#include <sys/types.h>
#define int32 int32_t

#else
#ifdef __APPLE__
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#define int32 int  /* a data type that has 32 bits */

#endif /* __APPLE__ */
#endif /* __linux__ */
#endif /* _WIN32 */
#endif /* SGI */

union tabfudge
{
    double tf_d;
    int32 tf_i[2];
};


/* -------------------------- ssaw~ ------------------------------ */
static t_class *ssaw_class, *scalarssaw_class;
static float ssaw_array[1002];
#define SAW_ARRAY_LEN 1002

typedef struct _ssaw
{
	t_object x_obj;
    //from phasor~:
	double x_phase;
    float x_conv;
    float x_f;	    /* scalar frequency */
    float x_band;	/* band limit (Hertz)*/
} t_ssaw;

static void *ssaw_new(t_floatarg f)
{
    t_ssaw *x = (t_ssaw *)pd_new(ssaw_class);
    x->x_f = f;
    inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_float, gensym("ft1"));
    x->x_phase = 0;
    x->x_conv = 0;
    x->x_band = 22000;
    outlet_new(&x->x_obj, gensym("signal"));
    return (x);
}

static t_int *ssaw_perform(t_int *w)
{
    t_ssaw *x = (t_ssaw *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int i,n = (int)(w[4]);
    double dphase = x->x_phase + UNITBIT32;
    union tabfudge tf;
    int normhipart;
    float conv = x->x_conv;
	float band=x->x_band*.33;
    float *buf = ssaw_array;

    tf.tf_d = UNITBIT32;
    normhipart = tf.tf_i[HIOFFSET];
    tf.tf_d = dphase;

    for (i = 0; i < n; i++)
    //while (n--)
    {
		float phase,band2,findex /*= *in++*/;
		int index /*= findex*/;
		float frac,  a,  b,  c,  d, cminusb, *fp;
		
    	tf.tf_i[HIOFFSET] = normhipart;
    	band2=abs(*in);
		if(band2>999999) band2=999999; else if(band2<1) band2=1;
		band2=band/band2;
		dphase += *in++ * conv;
		/**out++*/phase = (tf.tf_d - UNITBIT32)-0.5;
		tf.tf_d = dphase;

		findex=phase*band2;
		if(findex>0.5) findex=0.5; else if(findex<-0.5) findex=-0.5;

		/*findex=findex*1000+501;
		index=findex;*/
		/*if (index < 1)
		    index = 1, frac = 0;
		else if (index > maxindex)
		    index = maxindex, frac = 1;
		else*/ frac = findex - index;
		/*fp = buf + index;
		a = fp[-1];
		b = fp[0];
		c = fp[1];
		d = fp[2];
		cminusb = c-b;
		*out++ = 0.5+ phase - (
				 b + frac * ( cminusb - 0.1666667f * (1.-frac) * (
				(d - a - 3.0f * cminusb) * frac + (d + 2.0f*a - 3.0f*b))) 
				);*/
		*out++ = 0.5+ phase - buf[(int)(findex*1000+501)];
    }

    tf.tf_i[HIOFFSET] = normhipart;
    x->x_phase = tf.tf_d - UNITBIT32;
    return (w+5);
}

static void ssaw_dsp(t_ssaw *x, t_signal **sp)
{
    x->x_conv = 1./sp[0]->s_sr;
    dsp_add(ssaw_perform, 4, x, sp[0]->s_vec, sp[1]->s_vec, sp[0]->s_n);
}

static void ssaw_ft1(t_ssaw *x, t_float f)
{
    x->x_phase = f;
}

static void ssaw_initarray(void)
{
	int i;
	float j;

	for(i=0;i<1002;i++){
		j=(i-1)*M_PI/1000.0; //period 2000 sample, 1 sample back phase
		ssaw_array[i]= 0.57692* 
			(-1*cos(j) + 0.333333*cos(j*3.0) -0.2* cos(j*5.0));
	}
}

void ssaw_tilde_setup(void)
{
    ssaw_class = class_new(gensym("ssaw~"), (t_newmethod)ssaw_new, 0,
    	sizeof(t_ssaw), 0, A_DEFFLOAT, 0);
    CLASS_MAINSIGNALIN(ssaw_class, t_ssaw, x_f);
    class_addmethod(ssaw_class, (t_method)ssaw_dsp, gensym("dsp"), 0);
    class_addmethod(ssaw_class, (t_method)ssaw_ft1,
    	gensym("ft1"), A_FLOAT, 0);
	ssaw_initarray();

}