Hey Patrick,
So, I've probably screwed up some stuff here with GSL, but I've been learning a 
hell of a lot with regards to vector operations, and porting the vDSP 
accelerate stuff to GSL. Eventually I suppose it would be good to get rid of 
that, and only use standard libs, except where the FFT stuff is concerned when 
we'll use either d_fft_fftdg.c or d_fft_fftw.c from pd src.
Also, I was reading an article about dereferencing, and I think I might have 
used &x-> where just x-> is needed for gsl_vector_float operations, and I know 
the initialization functions gsl_vector_float_calloc() are wrong sized - I need 
to work out which need to be MAX_ORDER and which to be sized and initialized at 
size = nfft or size = nfft02 (It's 10pm and I'm also trying to master an album!)

Well, it's messy right now. Don't bust a gut on working all of it out, but if 
you can speed me up by spotting a few things it'd be appreciated.
I feel like I've absorbed a couple of megabytes in my head since the weekend.

Cheers,Ed
Enclosed - the original Max code and my incomplete hash of it - but I'm trying 
to port the hardest bit I think...Things will be renamed before a release to 
acknowledge Mark Cartwright, (as in mbc.lpc~ for the original MSP extern) but I 
haven't decided on the namespace options yet.


_-_-_-_-_-_-_-^-_-_-_-_-_-_-_

For Lone Shark releases, Pure Data software and published Research, go to 
http://sharktracks.co.uk  

    On Tuesday, 6 March 2018, 19:41:16 GMT, Pagano, Patrick 
<patrick.pag...@uconn.edu> wrote:  
 
 #yiv6390742159 #yiv6390742159 -- P 
{margin-top:0;margin-bottom:0;}#yiv6390742159 
Let me know how i can help ed




pp




Patrick Pagano B.S, M.F.A

Assistant Professor in Residence

Digital Media & Design


Web & Interactive Technologies

UCONN ECE Faculty Coordinator

University of Connecticut, Stamford

(352)-226-2016
From: Pd-list <pd-list-boun...@lists.iem.at> on behalf of Simon Iten 
<itensi...@gmail.com>
Sent: Tuesday, March 6, 2018 2:22:54 PM
To: Ed Kelly
Cc: pd list
Subject: Re: [PD] Porting Max MSP externals to Pure Data cool, thanks!

On 6 Mar 2018, at 17:25, Ed Kelly <morph_2...@yahoo.co.uk> wrote:
Encoding and decoding of LPC streams.I know you eagerly anticipate this! I have 
a lot of work to do...Ed


_-_-_-_-_-_-_-^-_-_-_-_-_-_-_

For Lone Shark releases, Pure Data software and publishedResearch, go to 
http://sharktracks.co.uk 

On Sunday, 4 March 2018, 22:11:46 GMT, Simon Iten <itensi...@gmail.com> wrote:

you are my hero. what will they do exactly? lcp encoding, decoding? lcp-streams 
to audio?
cheers

On 4 Mar 2018, at 15:30, Ed Kelly via Pd-list <pd-list@lists.iem.at> wrote:
Hi List,
I'm porting a library of LPC externals from Max/MSP to Pd.I wonder if someone 
could point me towards an example of another MSP external that has successfully 
been ported to Pd - preferably a DSP external that has creation arguments, with 
code for both so I can identify the differences between coding for the two 
platforms.
Cheers,Ed

_-_-_-_-_-_-_-^-_-_-_-_-_-_-_

For Lone Shark releases, Pure Data software and published Research, go 
tohttp://sharktracks.co.uk _______________________________________________
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management -> 
https://lists.puredata.info/listinfo/pd-list




  
/**
	@file
	mbc.lpc~ - an MSP object shell
	mark cartwright - mcartwri...@gmail.com	

	@ingroup	lpcToolkit	
*/

#include "ext.h"							// standard Max include, always required (except in Jitter)
#include "ext_obex.h"						// required for new style objects
#include "z_dsp.h"							// required for MSP objects
#include <math.h>
#include <Accelerate/Accelerate.h>

#define MAX_ORDER 200
#define NLOG2(x) (ceil(log2(x)))
#define POW2(x) (1 << x);
#define DEFAULT_FS 44100
#define DEFAULT_FRAMERATE 100
#define DEFAULT_V_SIZE 64

////////////////////////// object struct
typedef struct _lpc 
{
	t_pxobject					ob;					// the object itself (t_pxobject in MSP)
	t_float*					l_frame_buff;		//input frame buffer
	t_float*					l_winframe_buff;	//windowed input frame buffer
	t_float*					l_outCoeff_buff;	//coefficient signal
	t_float*					l_outParcor_buff; 	//PARCOR coeffs
	t_float*					l_outError_buff;	//error signal
	t_float*					l_win;				//analysis window
	t_float*					l_R;
	double*						l_W;
	double*						l_E;
	double*						l_K;
	double 						l_G;
	double**					l_A;
	double 						l_x1;				//last samples of pre-emphasis filter
	float 						l_b1;				//pre-emph coefficient
	int 						l_order;			//predictor order
	int 						l_order_max;		//max order according to fs = order * frame_rate
	int 						l_preemph;			//pre-epmhasis filter on/off
	int 						l_frame_rate;		//analysis frame rate
	int 						l_frame_size;		//analysis frame size, where fs = frame_rate * frame_size * 2
	int 						l_hop_size;			//hop_size = frame_size * 2 (b/c of overlap)
	int 						l_inframe_idx;		//current inframe buffer index
	int 						l_outframe_idx;		//current outframe buffer index
	long 						l_v_size;			//vector size
	float 						l_fs;				//sampling rate
	int 						l_maxnfft;			//fft length
	int 						l_log2nfft;			//log2(fft length)
	FFTSetup 					l_fftsetup;			//FFTSetup for vDSP FFT functions
	DSPSplitComplex 			l_fftSplitTemp;		//temp split complex vector structure
} t_lpc;

///////////////////////// function prototypes
//// standard set
void *lpc_new(long order, long framerate, long preemph);
void lpc_free(t_lpc *x);
void lpc_free_arrays(t_lpc *x);
void lpc_assist(t_lpc *x, void *b, long m, long a, char *s);

void lpc_dsp(t_lpc *x, t_signal **sp, short *count);
t_int *lpc_perform(t_int *w);

void lpc_order(t_lpc *x, int n);
void lpc_preemph(t_lpc *x, int n);
void lpc_init(t_lpc *x);
void lpc_hanning(t_lpc *x);
void lpc_hamming(t_lpc *x);
void lpc_bartlett(t_lpc *x);

//////////////////////// global class pointer variable
void *lpc_class;


int main(void)
{	
	// object initialization, note the use of dsp_free for the freemethod, which is required
	// unless you need to free allocated memory, in which case you should call dsp_free from
	// your custom free function.

	// OLD METHOD
	// setup((t_messlist **)&lpc_class, (method)lpc_new, (method)dsp_free, (short)sizeof(t_lpc), 0L, A_GIMME, 0);
	// addfloat((method)lpc_float);
	// you need this
	// addmess((method)lpc_dsp,				"dsp",			A_CANT, 0);
    // addmess((method)lpc_assist,			"assist",		A_CANT, 0);  
	// you need this
    // dsp_initclass();
	
	// NEW METHOD
	t_class *c;
	
	c = class_new("mbc.lpc~", (method)lpc_new, (method)lpc_free, (long)sizeof(t_lpc), 0L, A_DEFLONG, A_DEFLONG, A_DEFLONG, 0); //arglist: order, framerate, preemph
	
	class_addmethod(c, (method)lpc_dsp, "dsp", A_CANT, 0);
	class_addmethod(c, (method)lpc_order,"order",A_DEFLONG,0);
	class_addmethod(c, (method)lpc_preemph,"preemph",A_DEFLONG,0);
	class_addmethod(c, (method)lpc_assist,"assist",A_CANT,0);
	
	class_dspinit(c);				// new style object version of dsp_initclass();
	class_register(CLASS_BOX, c);	// register class as a box class
	lpc_class = c;
	
	return 0;
}

t_int *lpc_perform(t_int *w) 
{
	t_lpc *x = (t_lpc *) (w[1]);
	t_float *in = (t_float *)(w[2]);
	t_float *out_error = (t_float *)(w[3]);
	t_float *out_gain = (t_float *)(w[4]);
	t_float *out_coeff = (t_float *)(w[5]);
	t_float *out_parcor = (t_float *)(w[6]);
	t_float *out_index = (t_float *)(w[7]);
	int n = (int)(w[8]);
	int p = x->l_order;
	int length = x->l_frame_size;
	int log2nfft = NLOG2(length+p+1);
	int nfft = POW2(log2nfft);
	int nfftO2 = nfft/2;
	float recipn = 1.0/nfft;
	int inframeidx = x->l_inframe_idx;
	int outframeidx = x->l_outframe_idx;
	
	int i, j, i1, ji;
	int in_idx = 0, out_idx = 0;
	double val;
	float scale;
	
	if (x->l_preemph) 
	{
		while (n--) 
		{
			val = in[in_idx];
			in[in_idx] = val + x->l_b1 * x->l_x1;
			x->l_x1 = val;
			in_idx++;
		}
		n = (int)(w[8]);
		in_idx = 0;
	}
	
	while (n--) 
	{
		if (inframeidx < length) {
			//copy input into frame buff
			x->l_frame_buff[inframeidx] = in[in_idx];
			
			out_gain[out_idx] = x->l_G;
			out_error[out_idx] = x->l_outError_buff[outframeidx]; //for now
			if (outframeidx < x->l_order) {
				out_coeff[out_idx] = x->l_outCoeff_buff[outframeidx];
				out_parcor[out_idx] = x->l_outParcor_buff[outframeidx];
				out_index[out_idx] = outframeidx + 1;
			} else {
				out_coeff[out_idx] = 0.0;
				out_parcor[out_idx] = 0.0;
				out_index[out_idx] = 0;
			}
			
			inframeidx++;
			in_idx++;
			outframeidx++;
			out_idx++;
		} else {
			//perform durbin-levinson - for right now, just count to the order---------------
			//clear memory, is this necessary?
			for (i=0; i < p+1; i++){
				x->l_R[i] = 0.0;
				x->l_W[i] = 0.0;
				x->l_E[i] = 0.0;
				x->l_K[i] = 0.0;			
			}
			for(i=0; i<=p; i++) {
				for(j=0; j < p; j++) x->l_A[i][j] = 0.0;
			}
			//window frame buff
			vDSP_vmul(x->l_frame_buff, 1, x->l_win, 1, x->l_winframe_buff, 1, length);
			#ifdef DEBUG
				for(i=0;i<length;i++)cpost("\nwinframe(%d) = %g;",i+1,x->l_winframe_buff[i]);
			#endif
			
			//create r from auto correlation
			if ((2*nfft*log2(nfft)+nfft) > length*p) { //NOTE: change this to update only when order is changed!
			//time domain method
				for(i=0; i < p+1; i++) vDSP_dotpr(x->l_winframe_buff,1,x->l_winframe_buff+i,1,x->l_R+i,length-i);
			} else {
			//frequency domain method
				// zero pad
				vDSP_vclr(x->l_winframe_buff+length,1,nfft-length);
				//convert to split complex vector
				vDSP_ctoz( ( DSPComplex * ) x->l_winframe_buff, 2, &x->l_fftSplitTemp, 1, nfftO2);
				//perform forward in place fft
				vDSP_fft_zrip(x->l_fftsetup, &x->l_fftSplitTemp, 1, log2nfft, FFT_FORWARD);
				//scaling
				scale = 0.5;
				vDSP_vsmul(x->l_fftSplitTemp.realp,1,&scale,x->l_fftSplitTemp.realp,1,nfftO2);
				vDSP_vsmul(x->l_fftSplitTemp.imagp,1,&scale,x->l_fftSplitTemp.imagp,1,nfftO2);
				//compute PSD
				vDSP_zvmags(&x->l_fftSplitTemp,1,x->l_fftSplitTemp.realp,1,nfftO2);
				//clear imaginary part
				vDSP_vclr(x->l_fftSplitTemp.imagp,1,nfftO2);
				//perform inverse in place fft
				vDSP_fft_zrip(x->l_fftsetup, &x->l_fftSplitTemp, 1, log2nfft, FFT_INVERSE);
				//scaling
				vDSP_vsmul(x->l_fftSplitTemp.realp,1,&recipn,x->l_fftSplitTemp.realp,1,nfftO2);
				vDSP_vsmul(x->l_fftSplitTemp.imagp,1,&recipn,x->l_fftSplitTemp.imagp,1,nfftO2);
				//convert back to real number vector
				vDSP_ztoc(&x->l_fftSplitTemp, 1, (DSPComplex *)x->l_R, 2, nfftO2);	
			}
			
			
			/*for(i=0; i < p+1; i++) {
				x->l_R[i] = 0.0;
				for(j=0; j < length - i; j++) {
					x->l_R[i] += x->l_winframe_buff[j] * x->l_winframe_buff[i+j];
					//x->l_R[i] += x->l_win[j] * x->l_win[i+j];
				}
			}*/
			#ifdef DEBUG
					for(i=0;i< p+1; i++) cpost("\nR(%d) = %f;",i+1,x->l_R[i]);
			#endif
			
			x->l_W[0] = x->l_R[1];
			x->l_E[0] = x->l_R[0];
			
			for (i = 1; i <= p; i++) {
				x->l_K[i] = x->l_W[i-1] / x->l_E[i-1];
				
				x->l_A[i][i] = x->l_K[i];
				i1 = i - 1;
				if (i1 >= 1) {
					for (j = 1; j <=i1; j++) {
						ji = i - j;
						x->l_A[j][i] = x->l_A[j][i1] - x->l_K[i] * x->l_A[ji][i1];
					}				
				}
				
				x->l_E[i] = x->l_E[i-1] * (1.0 - x->l_K[i] * x->l_K[i]);
				
				if (i != p) {
					x->l_W[i] = x->l_R[i+1];
					for (j = 1; j <= i; j++) {
						x->l_W[i] -= x->l_A[j][i] * x->l_R[i-j+1];
					}
				}
			}
			
			x->l_G = sqrt(x->l_E[p]);
			for (i=0; i < p; i++) {
				x->l_outCoeff_buff[i] = (float)(x->l_A[i+1][p]);
				x->l_outParcor_buff[i] = (float)(x->l_K[i+1]);
				#ifdef DEBUG
					//cpost("\nParcor(%d) = %g;",i+1,x->l_K[i+1]);
					//cpost("\nCoeff(%d) = %g;",i+1,x->l_A[i+1][p]);
				#endif
			}
			
			//--------------------------------------------------------------------------------
			
			//copy right side to left side, move delayed input to output
			for (i=0; i < x->l_hop_size; i++) {
				x->l_outError_buff[i] = x->l_frame_buff[i];
				x->l_frame_buff[i] = x->l_frame_buff[i + x->l_hop_size];
			}
			
			inframeidx = x->l_hop_size;
			outframeidx = 0;
			n++; //to avoid skipping a sample (since we already decremented
			while (n--) {
				x->l_frame_buff[inframeidx] = in[in_idx];
				
				out_gain[out_idx] = (float)(x->l_G);
				out_error[out_idx] = x->l_outError_buff[outframeidx]; //for now
				if (outframeidx < x->l_order) {
					out_coeff[out_idx] = x->l_outCoeff_buff[outframeidx];
					out_parcor[out_idx] = x->l_outParcor_buff[outframeidx];
					out_index[out_idx] = (float)(outframeidx + 1);
				} else {
					out_coeff[out_idx] = 0.0;
					out_parcor[out_idx] = 0.0;
					out_index[out_idx] = 0;
				}
				
				inframeidx++;
				in_idx++;
				outframeidx++;
				out_idx++;
			}
			break;
		}
	}
	
	x->l_inframe_idx = inframeidx;
	x->l_outframe_idx = outframeidx;

	return (w+9);
}

void lpc_dsp(t_lpc *x, t_signal **sp, short *count)
{
	if (sp[0]->s_n != x->l_v_size || sp[0]->s_sr != x->l_fs) {
		x->l_v_size = sp[0]->s_n;
		x->l_fs = sp[0]->s_sr;
		lpc_init(x);
	}
	
	dsp_add(lpc_perform, 8, x, sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[3]->s_vec, sp[4]->s_vec, sp[5]->s_vec, sp[0]->s_n);
}

void lpc_assist(t_lpc *x, void *b, long msg, long arg, char *dst)
{
	if (msg==ASSIST_INLET) {
		switch (arg) {
			case 0: sprintf(dst,"(signal) LPC Input"); break;
		}
	}
	else {
		switch (arg) {
			case 0: sprintf(dst,"(signal) Error Signal"); break;
			case 1: sprintf(dst,"(signal) Filter Gain"); break;
			case 2: sprintf(dst,"(signal) Filter Coefficients"); break;
			case 3: sprintf(dst,"(signal) PARCOR Coefficients"); break;
			case 4: sprintf(dst,"(signal) Coefficient Index"); break;
		}
	}
}

void lpc_free(t_lpc *x) {
	dsp_free((t_pxobject *) x);
	lpc_free_arrays(x);
}

void lpc_preemph(t_lpc *x, int n) {
	x->l_preemph = n;
}

void lpc_order(t_lpc *x, int n) {
	if (x->l_order_max < n) {
		error("mbc.lpc~: framerate * order must be less than or equal to sampling rate.  At the current setting maxorder is %d.  For a higher order, lower the framerate.", x->l_order_max);
		x->l_order = x->l_order_max;
	} else {
		x->l_order = n;
	}
	//lpc_init(x) ?
}

void lpc_init(t_lpc *x) {
	int i;
	
	if (x->l_fs < x->l_frame_rate * x->l_order) {
		x->l_frame_rate = (long)(x->l_fs / x->l_order);
		error("mbc.lpc~: framerate * order must be less than or equal to sampling rate.  framerate has been changed to %d", x->l_frame_rate);
	}
	
	x->l_hop_size = (int)(x->l_fs / x->l_frame_rate);
	
	if (x->l_v_size > x->l_hop_size) {
		error("mbc.lpc~: warning: frame_size is less than vector size.  this may cause errors, please reduce the frame rate or increase the vector size");
	}
	
	x->l_order_max = x->l_fs / x->l_frame_rate;
	x->l_frame_size = x->l_hop_size * 2;
	x->l_inframe_idx = x->l_hop_size;
	x->l_outframe_idx = 0;
	x->l_log2nfft = (int) NLOG2(x->l_frame_size+MAX_ORDER+1);
	x->l_maxnfft = POW2(x->l_log2nfft);
	
	// free memory if needed
	lpc_free_arrays(x);
	
	//allocate memory
	x->l_frame_buff = (t_float *) sysmem_newptrclear( x->l_frame_size * sizeof(t_float));
	x->l_winframe_buff = (t_float *) sysmem_newptrclear( x->l_maxnfft * sizeof(t_float));
	x->l_outCoeff_buff = (t_float *) sysmem_newptrclear(MAX_ORDER * sizeof(t_float));
	x->l_outParcor_buff = (t_float *) sysmem_newptrclear(MAX_ORDER * sizeof(t_float));
	x->l_outError_buff = (t_float *) sysmem_newptrclear( x->l_hop_size * sizeof(t_float));
	x->l_win = (t_float *) sysmem_newptrclear( x->l_frame_size * sizeof(t_float));
	x->l_R = (t_float *) sysmem_newptrclear( x->l_maxnfft * sizeof(t_float));
	x->l_W = (double *) sysmem_newptrclear( (MAX_ORDER + 1) * sizeof(double));
	x->l_E = (double *) sysmem_newptrclear( (MAX_ORDER + 1) * sizeof(double));
	x->l_K = (double *) sysmem_newptrclear( (MAX_ORDER + 1) * sizeof(double));
	x->l_A = (double **) sysmem_newptrclear( (MAX_ORDER + 1) * sizeof(double*));
	for(i=0; i<MAX_ORDER; i++) {
		x->l_A[i] = (double *)sysmem_newptrclear( (MAX_ORDER + 1) * sizeof(double));
	}
	x->l_fftSplitTemp.realp = (float *)sysmem_newptrclear( (x->l_maxnfft / 2) * sizeof(float));
	x->l_fftSplitTemp.imagp = (float *)sysmem_newptrclear( (x->l_maxnfft / 2) * sizeof(float));
		
	x->l_x1 = 0;
	x->l_G = 0.0;
	
	x->l_b1 = -0.98;
	
	//calculate window
	lpc_hamming(x);
	//create FFTsetups
	x->l_fftsetup = vDSP_create_fftsetup(x->l_log2nfft,FFT_RADIX2);
	if (!x->l_fftsetup) error("mbc.lpc~: not enough available memory");
	
}

void lpc_free_arrays(t_lpc *x)
{
	int i;
	
	if (x->l_frame_buff) {
		sysmem_freeptr(x->l_frame_buff);
		x->l_frame_buff = NULL;
	}
	if (x->l_winframe_buff) {
		sysmem_freeptr(x->l_winframe_buff);
		x->l_winframe_buff = NULL;
	}
	if (x->l_outCoeff_buff) {
		sysmem_freeptr(x->l_outCoeff_buff);
		x->l_outCoeff_buff = NULL;
	}
	if (x->l_outError_buff) {
		sysmem_freeptr(x->l_outError_buff);
		x->l_outError_buff = NULL;
	}
	if (x->l_outParcor_buff) {
		sysmem_freeptr(x->l_outParcor_buff);
		x->l_outParcor_buff = NULL;
	}
	if (x->l_win) {
		sysmem_freeptr(x->l_win);
		x->l_win = NULL;
	}
	if (x->l_R) {
		sysmem_freeptr(x->l_R);
		x->l_R = NULL;
	}
	if (x->l_W) {
		sysmem_freeptr(x->l_W);
		x->l_W = NULL;
	}
	if (x->l_E) {
		sysmem_freeptr(x->l_E);
		x->l_E = NULL;
	}
	if (x->l_K) {
		sysmem_freeptr(x->l_K);
		x->l_K = NULL;
	}
	if (x->l_A) {
		for(i=0; i<MAX_ORDER; i++) {
			if (x->l_A[i]) {
				sysmem_freeptr(x->l_A[i]);
				x->l_A[i] = NULL;
			}
		}
		sysmem_freeptr(x->l_A);
		x->l_A = NULL;
	}
	if (x->l_fftSplitTemp.realp) {
		sysmem_freeptr(x->l_fftSplitTemp.realp);
		x->l_fftSplitTemp.realp = NULL;
	}
	if (x->l_fftSplitTemp.imagp) {
		sysmem_freeptr(x->l_fftSplitTemp.imagp);
		x->l_fftSplitTemp.imagp = NULL;
	}

	if (x->l_fftsetup) {
		vDSP_destroy_fftsetup(x->l_fftsetup);
		x->l_fftsetup = NULL;
	}
}

void lpc_hanning(t_lpc *x) {
	int N = x->l_frame_size;
	int n;
	
	for (n=0; n<N; n++) {
		x->l_win[n] = 0.5 * (1.0 - cos((TWOPI*(n+1))/(N+1)));
	}
}

void lpc_hamming(t_lpc *x) {
	int N = x->l_frame_size;
	int n;
	
	for (n=0; n<N; n++) {
		x->l_win[n] = (0.54 - 0.46*cos((TWOPI*n)/(N-1)));
	}
}

void lpc_bartlett(t_lpc *x) {
	int N = x->l_frame_size;
	int n;
	
	for (n=0; n<N; n++) {
		x->l_win[n] = 2.0/(N - 1) * ((N-1)/2 - fabs(n - (N-1)/2.0));
	}
}

void *lpc_new(long order, long framerate, long preemph)
{
	t_lpc *x = NULL;
	
	if (order == 0)
	{
		error("mbc.lpc~: must specify the lpc order");
		return NULL;
	}
	
	if (framerate == 0)
	{
		framerate = DEFAULT_FRAMERATE;
	}
	
	if (x = (t_lpc *)object_alloc(lpc_class)) {
		dsp_setup((t_pxobject *)x,1);
		outlet_new((t_pxobject *)x, "signal");
		outlet_new((t_pxobject *)x, "signal");
		outlet_new((t_pxobject *)x, "signal");
		outlet_new((t_pxobject *)x, "signal");
		outlet_new((t_pxobject *)x, "signal");
		
		// from arguments:
		x->l_frame_rate = framerate;
		x->l_order = order;
		x->l_preemph = preemph;
		
		// system variables (note: these need to be checked again in lpc_dsp since they may differ if in a poly~ for instance)
		x->l_fs = sys_getsr();
		if (x->l_fs == 0)
		{
			x->l_fs = DEFAULT_FS;
		}
		
		x->l_v_size = sys_getblksize();
		if (x->l_v_size == 0)
		{
			x->l_v_size = DEFAULT_V_SIZE;
		}
		
		lpc_init(x);
	}
	return (x);
}
/*
 *   LPC Toolkit
 *   By Mark Cartwright
 *
 *
 *
 *
 */

#include "m_pd.h"
#include <math.h>
//#include <gsl/gsl_math.h>
#include <gsl/gsl_vector_float.h>
#include <gsl/gsl_blas.h>

static t_class *mbc_lpc_class;

#define MAX_ORDER 200
#define NLOG2(x) (ceil(log2(x)))
#define POW2(x) (1 << x);
//#define DEFAULT_FS 44100
#define DEFAULT_FRAMERATE 100
#define DEFAULT_V_SIZE 64
#define DEFAULT_ORDER 32

typedef struct _DSPcomplex
{
  gsl_vector_float* real;
  gsl_vector_float* imag;
} t_DSPcomplex;

//#define REAL(z,i) gsl_vector_float_set(z,2*(i))
//#define IMAG(z,i) gsl_vector_float_setz(z,2*(i)+1)

typedef struct _lpc 
{
  t_object          x_obj;	           // the object itself (t_pxobject in MSP)
  //t_float*          l_frame_buff;	   //input frame buffer
  gsl_vector_float*          l_frame_buff;	   //input frame buffer
  //t_float*          l_winframe_buff;       //windowed input frame buffer
  gsl_vector_float*          l_winframe_buff;       //windowed input frame buffer
  //t_float*	    l_outCoeff_buff;       //coefficient signal
  gsl_vector_float*	    l_outCoeff_buff;       //coefficient signal
  //t_float*	    l_outParcor_buff; 	   //PARCOR coeffs
  gsl_vector_float*	    l_outParcor_buff; 	   //PARCOR coeffs
  //t_float*          l_outError_buff;	   //error signal
  gsl_vector_float*          l_outError_buff;	   //error signal
  //t_float*          l_win;	           //analysis window
  gsl_vector_float*          l_win;	           //analysis window
  //t_float*	    l_R;
  gsl_vector_float*	    l_R;
  double*	    l_W;
  double*	    l_E;
  double*	    l_K;
  double 	    l_G;
  double**	    l_A;
  double 	    l_x1;                  //last samples of pre-emphasis filter
  float 	    l_b1;	           //pre-emph coefficient
  int 	            l_order;	           //predictor order
  int 	            l_order_max;           //max order according to fs = order * frame_rate
  int 	            l_preemph;             //pre-epmhasis filter on/off
  int 	            l_frame_rate;          //analysis frame rate
  int 	            l_frame_size;          //analysis frame size, where fs = frame_rate * frame_size * 2
  int 	            l_hop_size;            //hop_size = frame_size * 2 (b/c of overlap)
  int 	            l_inframe_idx;         //current inframe buffer index
  int 	            l_outframe_idx;	   //current outframe buffer index
  long 	            l_v_size;	           //vector size
  float 	    l_fs;                  //sampling rate
  int 	            l_maxnfft;             //fft length
  int 	            l_log2nfft;	           //log2(fft length)
  int initOrder, initFramerate, initPreemph; // initialization values
  //  FFTSetup          l_fftsetup;            //FFTSetup for vDSP FFT functions
  t_DSPcomplex      split;
  //  DSPSplitComplex   l_fftSplitTemp;        //temp split complex vector structure
} t_lpc;

t_int *lpc_perform(t_int *w) 
{
  t_lpc   *x          = (t_lpc *)  (w[1]);
  t_float *in         = (t_float *)(w[2]);
  t_float *out_error  = (t_float *)(w[3]);
  t_float *out_gain   = (t_float *)(w[4]);
  t_float *out_coeff  = (t_float *)(w[5]);
  t_float *out_parcor = (t_float *)(w[6]);
  t_float *out_index  = (t_float *)(w[7]);
  int n               = (int)      (w[8]);
  int p = x->l_order;
  int length = x->l_frame_size;
  int log2nfft = NLOG2(length+p+1);
  int nfft = POW2(log2nfft);
  int nfftO2 = nfft/2;
  t_float recipn = 1.0/nfft;
  int inframeidx = x->l_inframe_idx;
  int outframeidx = x->l_outframe_idx;
  
  int i, j, i1, ji;
  int in_idx = 0, out_idx = 0;
  double val;
  t_float scale, rCorr, re, im;
  gsl_vector_float *vectorBuffer;
  vectorBuffer = gsl_vector_float_calloc(MAX_ORDER);	   //input frame buffer
  
  
  if (x->l_preemph)
    {
      while (n--) 
	{
	  val = in[in_idx];
	  in[in_idx] = val + x->l_b1 * x->l_x1;
	  x->l_x1 = val;
	  in_idx++;
	}
      n = (int)(w[8]);
      in_idx = 0;
    }
  
  while (n--) 
    {
      if (inframeidx < length) {
	//copy input into frame buff
	//also copy into winframe buff since GSL puts the result into one of the arg vectors
	gsl_vector_float_set(&x->l_frame_buff,inframeidx,in[in_idx]);// = x->l_frame_buff[inframeidx] = in[in_idx];
	
	out_gain[out_idx] = x->l_G;
	out_error[out_idx] = gsl_vector_float_get(&x->l_outError_buff,outframeidx); //for now
	if (outframeidx < x->l_order) {
	  out_coeff[out_idx] = gsl_vector_float_get(&x->l_outCoeff_buff,outframeidx);
	  out_parcor[out_idx] = gsl_vector_float_get(&x->l_outParcor_buff,outframeidx);
	  out_index[out_idx] = outframeidx + 1;
	} else {
	  out_coeff[out_idx] = 0.0;
	  out_parcor[out_idx] = 0.0;
	  out_index[out_idx] = 0;
	}
	
	inframeidx++;
	in_idx++;
	outframeidx++;
	out_idx++;
      } else {
	gsl_vector_float_memcpy(&x->l_winframe_buff,&x->l_frame_buff);
	 
	//perform durbin-levinson - for right now, just count to the order---------------
	//clear memory, is this necessary?
	for (i=0; i < p+1; i++){
	  x->l_R[i] = 0.0;
	  x->l_W[i] = 0.0;
	  x->l_E[i] = 0.0;
	  x->l_K[i] = 0.0;			
	}
	for(i=0; i<=p; i++) {
	  for(j=0; j < p; j++) x->l_A[i][j] = 0.0;
	}
	//window frame buff
	//vDSP_vmul(x->l_frame_buff, 1, x->l_win, 1, x->l_winframe_buff, 1, length);
	gsl_vector_float_mul(&x->l_winframe_buff,&x->l_win);
#ifdef DEBUG
	for(i=0;i<length;i++)post("\nwinframe(%d) = %g;",i+1,gsl_vector_float_get(&x->l_winframe_buff,i));
#endif
	
	//create r from auto correlation
	if ((2*nfft*log2(nfft)+nfft) > length*p) { //NOTE: change this to update only when order is changed!
	  //time domain method
	  //for(i=0; i < p+1; i++) vDSP_dotpr(x->l_winframe_buff,1,x->l_winframe_buff+i,1,x->l_R+i,length-i);
	  //gsl_vector_float_memcpy(vectorBuffer,&x->l_winframe_buff);
	  for(i=0; i < p+1; i++)
	    {
	      gsl_blas_sdot(&x->l_winframe_buff,&x->l_winframe_buff+i,&rCorr);
	      gsl_vector_float_set(&x->l_R,i,rCorr);
	    }
	} else {
	  //frequency domain method
	  // zero pad
	  //vDSP_vclr(x->l_winframe_buff+length,1,nfft-length);
	  for(i=length;i<nfft;i++)
	    {
	      gsl_vector_float_set(&x->l_winframe_buff,i,0);
	    }
	  //for(i=0;i<nfft02;i++)
	  //  {
	  //    //convert to split complex vector
	  //    gsl_vector_float_set(&x->split.real,i,gsl_vector_float_get(x->l_winframe_buff,i*2));
	  //    gsl_vector_float_set(&x->split.imag,i,gsl_vector_float_get(x->l_winframe_buff,i*2+1));
	  //  }
	  //vDSP_ctoz( ( DSPComplex * ) x->l_winframe_buff, 2, &x->l_fftSplitTemp, 1, nfftO2);
	  //perform forward in place fft
	  gsl_fft_complex_radix2_forward(&x->l_winframe_buff,1,nfft);
	  
	  //vDSP_fft_zrip(x->l_fftsetup, &x->l_fftSplitTemp, 1, log2nfft, FFT_FORWARD);
	  //scaling
	  scale = 0.5;

	  // Ed Kelly: scale real and imaginary parts and calculate magnitudes squared
	  for(i=0;i<nfft02;i++)
	    {
	      j = i*2;
	      re = gsl_vector_float_get(&x->l_winframe_buff,j);
	      im = gsl_vector_float_get(&x->l_winframe_buff,j+1);
	      re *= scale;
	      re = re * re;
	      im *= scale;
	      im = im * im;
	      //I can't see how we need the split struct now!!!
	      //gsl_vector_float_set(&x->split.real,i,re+im);
	      //gsl_vector_float_set(&x->split.imag,i,0.0);
	      gsl_vector_float_set(vectorBuffer,j,re+im);
	    }

	  //vDSP_vsmul(x->l_fftSplitTemp.realp,1,&scale,x->l_fftSplitTemp.realp,1,nfftO2);
	  //vDSP_vsmul(x->l_fftSplitTemp.imagp,1,&scale,x->l_fftSplitTemp.imagp,1,nfftO2);
	  //compute PSD
	  //vDSP_zvmags(&x->l_fftSplitTemp,1,x->l_fftSplitTemp.realp,1,nfftO2);

	  //clear imaginary part
	  //for(i=0;i<nfft02;i++)
	  //{
	      //IMAG(&x->l_winframe_buff,i) = 0.0;
	      //gsl_vector_float_set(&x->l_winframe_buff,i*2+1,0.0);
	  //}
	  //vDSP_vclr(x->l_fftSplitTemp.imagp,1,nfftO2);
	  //perform inverse in place fft
	  gsl_fft_complex_radix2_backward(vectorBuffer,1,nfft);
	  //vDSP_fft_zrip(x->l_fftsetup, &x->l_fftSplitTemp, 1, log2nfft, FFT_INVERSE);
	  //scaling
	  //vDSP_vsmul(x->l_fftSplitTemp.realp,1,&recipn,x->l_fftSplitTemp.realp,1,nfftO2);
	  //vDSP_vsmul(x->l_fftSplitTemp.imagp,1,&recipn,x->l_fftSplitTemp.imagp,1,nfftO2);
	  gsl_vector_float_scale(vectorBuffer,&recipn);
	  //convert back to real number vector
	  //vDSP_ztoc(&x->l_fftSplitTemp, 1, (DSPComplex *)x->l_R, 2, nfftO2);
	  gsl_vector_float_memcpy(x->l_R,vectorBuffer);
	}
	
	
	/*for(i=0; i < p+1; i++) {
	  x->l_R[i] = 0.0;
	  for(j=0; j < length - i; j++) {
	  x->l_R[i] += x->l_winframe_buff[j] * x->l_winframe_buff[i+j];
	  //x->l_R[i] += x->l_win[j] * x->l_win[i+j];
	  }
	  }*/
#ifdef DEBUG
	for(i=0;i< p+1; i++) post("\nR(%d) = %f;",i+1,x->l_R[i]);
#endif
	
	x->l_W[0] = gsl_vector_float_get(x->l_R,1);
	x->l_E[0] = gsl_vector_float_get(x->l_R,0);
	
	for (i = 1; i <= p; i++) {
	  x->l_K[i] = x->l_W[i-1] / x->l_E[i-1];
	  
	  x->l_A[i][i] = x->l_K[i];
	  i1 = i - 1;
	  if (i1 >= 1) {
	    for (j = 1; j <=i1; j++) {
	      ji = i - j;
	      x->l_A[j][i] = x->l_A[j][i1] - x->l_K[i] * x->l_A[ji][i1];
	    }				
	  }
	  
	  x->l_E[i] = x->l_E[i-1] * (1.0 - x->l_K[i] * x->l_K[i]);
	  
	  if (i != p) {
	    x->l_W[i] = gsl_vector_float_get(x->l_R,i+1);
	    for (j = 1; j <= i; j++) {
	      x->l_W[i] -= x->l_A[j][i] * gsl_vector_float_get(x->l_R,i-j+1);
	    }
	  }
	}
	
	x->l_G = sqrt(x->l_E[p]);
	for (i=0; i < p; i++) {
	  x->l_outCoeff_buff[i] = (float)(x->l_A[i+1][p]);
	  x->l_outParcor_buff[i] = (float)(x->l_K[i+1]);
#ifdef DEBUG
	  //post("\nParcor(%d) = %g;",i+1,x->l_K[i+1]);
	  //post("\nCoeff(%d) = %g;",i+1,x->l_A[i+1][p]);
#endif
	}
	
	//--------------------------------------------------------------------------------
	
	//copy right side to left side, move delayed input to output
	for (i=0; i < x->l_hop_size; i++) {
	  x->l_outError_buff[i] = x->l_frame_buff[i];
	  x->l_frame_buff[i] = x->l_frame_buff[i + x->l_hop_size];
	}
	
	inframeidx = x->l_hop_size;
	outframeidx = 0;
	n++; //to avoid skipping a sample (since we already decremented
	while (n--) {
	  x->l_frame_buff[inframeidx] = in[in_idx];
				
	  out_gain[out_idx] = (t_float)(x->l_G);
	  out_error[out_idx] = x->l_outError_buff[outframeidx]; //for now
	  if (outframeidx < x->l_order) {
	    out_coeff[out_idx] = x->l_outCoeff_buff[outframeidx];
	    out_parcor[out_idx] = x->l_outParcor_buff[outframeidx];
	    out_index[out_idx] = (t_float)(outframeidx + 1);
	  } else {
	    out_coeff[out_idx] = 0.0;
	    out_parcor[out_idx] = 0.0;
	    out_index[out_idx] = 0;
	  }
	  
	  inframeidx++;
	  in_idx++;
	  outframeidx++;
	  out_idx++;
	}
	break;
      }
    }
  
  x->l_inframe_idx = inframeidx;
  x->l_outframe_idx = outframeidx;
  
  return (w+9);
}





void lpc_hanning(t_lpc *x)
{
  int N = x->l_frame_size;
  int n;
	
  for (n=0; n<N; n++)
    {
      x->l_win[n] = 0.5 * (1.0 - cos((TWOPI*(n+1))/(N+1)));
    }
}

void lpc_hamming(t_lpc *x)
{
  int N = x->l_frame_size;
  int n;
	
  for (n=0; n<N; n++)
    {
      x->l_win[n] = (0.54 - 0.46*cos((TWOPI*n)/(N-1)));
    }
}

void lpc_bartlett(t_lpc *x)
{
  int N = x->l_frame_size;
  int n;

  for (n=0; n<N; n++)
    {
      x->l_win[n] = 2.0/(N - 1) * ((N-1)/2 - fabs(n - (N-1)/2.0));
    }
}

void lpc_preemph(t_lpc *x, int n)
{
  x->l_preemph = n;
}

void lpc_order(t_lpc *x, int n)
{
  if (x->l_order_max < n)
    {
      error("mbc.lpc~: framerate * order must be less than or equal to sampling rate.  At the current setting maxorder is %d.  For a higher order, lower the framerate.", x->l_order_max);
      x->l_order = x->l_order_max;
    }
  else
    {
      x->l_order = n;
    }
    //lpc_init(x) ?
}

void lpc_free_arrays(t_lpc *x)
{
  //int i;
	
  gsl_vector_float_free(x->l_frame_buff);
  x->l_frame_buff = NULL;
  gsl_vector_float_free(x->l_winframe_buff);
  x->l_winframe_buff = NULL;
  gsl_vector_float_free(x->l_outCoeff_buff);
  x->l_outCoeff_buff = NULL;
  gsl_vector_float_free(x->l_outError_buff);
  x->l_outError_buff = NULL;
  gsl_vector_float_free(x->l_outParcor_buff);
  x->l_outParcor_buff = NULL;
  gsl_vector_float_free(x->l_win);
  x->l_win = NULL;
  gsl_vector_float_free(x->l_R);
  x->l_R = NULL;
  if (x->l_W) {
    freebytes(x->l_W);
    x->l_W = NULL;
  }
  if (x->l_E) {
    freebytes(x->l_E);
    x->l_E = NULL;
  }
  if (x->l_K) {
    freebytes(x->l_K);
    x->l_K = NULL;
  }
  if (x->l_A) {		
    freebytes(x->l_A);
    x->l_A = NULL;
  }
  /*	if (x->l_fftSplitTemp.realp) {
	sysmem_freeptr(x->l_fftSplitTemp.realp);
	x->l_fftSplitTemp.realp = NULL;
	}
	if (x->l_fftSplitTemp.imagp) {
	sysmem_freeptr(x->l_fftSplitTemp.imagp);
	x->l_fftSplitTemp.imagp = NULL;
	}
	
	if (x->l_fftsetup) {
	vDSP_destroy_fftsetup(x->l_fftsetup);
	x->l_fftsetup = NULL;
	}*/
}

void *lpc_new(t_symbol *s, int argc, t_atom *argv)
{
  t_lpc *x = (t_lpc *)pd_new(lpc_class);

  //#define DEFAULT_FS 44100
  //#define DEFAULT_FRAMERATE 100
  //#define DEFAULT_V_SIZE 64
  //#define DEFAULT_ORDER 32
  /*if(argc >= 3)
    {
      x->initOrder = atom_getfloat(argv);
      x->initFramerate = atom_getfloat(argv+1);
      x->initPreemph = atom_getfloat(argv+2);
      }
      else
      if(argc == 2)*/
  if(argc >= 2)
    {
      x->initOrder = atom_getfloat(argv);
      //x->initFramerate = atom_getfloat(argv+1);
      x->initPreemph = atom_getfloat(argv+1);
    }
  else if(argc == 1)
    {
      x->initOrder = atom_getfloat(argv);
      //x->initFramerate = DEFAULT_FRAMERATE;
      x->initPreemph = 0;
    }
  else
    {
      x->initOrder = DEFAULT_ORDER;
      //x->initFramerate = DEFAULT_FRAMERATE;
      x->initPreemph = 0;
    }
      
  if(x->initOrder < 1)
    {
      post("Order must be from 1 to 200! Setting to %d",DEFAULT_ORDER);
      x->l_order = DEFAULT_ORDER;
    }
  else if(x->initOrder > MAX_ORDER)
    {
      post("Order cannot be greater than 200!");
      x->l_order = MAX_ORDER;
    }
  else
    {
      post("Setting order to %d",x->initOrder);
      x->l_order = x->initOrder;
    }
  x->l_fs = sys_getsr();
  x->l_frame_rate = x->l_fs / x->l_order;
  
  //int length = x->l_frame_size;
  //int log2nfft = NLOG2(length+p+1);
  //int nfft = POW2(log2nfft);
  //int nfftO2 = nfft/2;
    //MAX_ORDER is wrong - check fft size calculation!
  x->l_frame_buff = gsl_vector_float_calloc(MAX_ORDER);	   //input frame buffer
  x->l_winframe_buff = gsl_vector_float_calloc(MAX_ORDER);       //windowed input frame buffer
  x->l_outCoeff_buff = gsl_vector_float_calloc(MAX_ORDER);       //coefficient signal
  x->l_outParcor_buff = gsl_vector_float_calloc(MAX_ORDER); 	   //PARCOR coeffs
  x->l_outError_buff = gsl_vector_float_calloc(MAX_ORDER);	   //error signal
  x->l_win = gsl_vector_float_calloc(MAX_ORDER);	           //analysis window
  



}
_______________________________________________
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management -> 
https://lists.puredata.info/listinfo/pd-list

Reply via email to