dawn@0: /* 
dawn@0: * interp_inplace.c
dawn@0: * linear interpolation
dawn@0: * 
dawn@0: * Alain de Cheveigné, CNRS/Ircam	
dawn@0: * (c) 2003 CNRS
dawn@0: */
dawn@0: 
dawn@0: #include <math.h>
dawn@0: #include "mex.h"
dawn@0: 
dawn@0: /* #define MWRAP */
dawn@0: #include "mwrap_check.h"
dawn@0: 
dawn@0: /* Input Arguments */
dawn@0: #define X_IN	prhs[0]
dawn@0: #define IDX_IN	prhs[1]
dawn@0: 
dawn@0: /* Output Arguments */
dawn@0: 
dawn@0: static void interp_inplace(
dawn@0: 				double *xp,		/* input vector */
dawn@0: 				double *idxp,	/* index vector */
dawn@0: 				int m, 			/* rows input vector */
dawn@0: 				int midx		/* rows index vector */
dawn@0: 			) 
dawn@0: {
dawn@0: 	int j,idxi;
dawn@0: 	double a,b,idx,idxf;
dawn@0: 		
dawn@0: 	
dawn@0: 	for (j=0;j<midx;j++) {								/* for each index */
dawn@0: 		idx=GET(idxp[j]);
dawn@0: 		idxi=floor(idx);
dawn@0: 		idxf=idx-idxi;
dawn@0: 		if (idxi<0) {
dawn@0: 			a=GET(xp[0]);
dawn@0: 			b=GET(xp[1]);
dawn@0: 			SET(idxp[j]) = a + idx*(b-a);
dawn@0: 		} else if (idxi>=m) {
dawn@0: 			a=GET(xp[m-2]);
dawn@0: 			b=GET(xp[m-1]);
dawn@0: 			SET(idxp[j]) = b + (idx-m+1)*(b-a);
dawn@0: 		} else {
dawn@0: 			a=GET(xp[idxi]);
dawn@0: 			b=GET(xp[idxi+1]);
dawn@0: 			SET(idxp[j]) = a + idxf*(b-a);
dawn@0: 		}
dawn@0: 	}
dawn@0: 	
dawn@0: 	return;						
dawn@0: }
dawn@0: 
dawn@0: void mexFunction(
dawn@0: 				int nlhs, mxArray *plhs[],
dawn@0: 				int nrhs, const mxArray *prhs[]
dawn@0: 				)
dawn@0: {
dawn@0: 	double	*xp, *idxp; 
dawn@0: 	int nx, mx, nidx, midx;
dawn@0: 	
dawn@0: 	/* Check for proper number of arguments */
dawn@0: 	if (nrhs !=2 ) {
dawn@0: 		mexErrMsgTxt("INTERP_INPLACE takes 2 input arguments");
dawn@0: 	}
dawn@0: 	
dawn@0: 	/* Check type of input */
dawn@0: 	if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || 
dawn@0: 		mxIsSparse(X_IN)	|| !mxIsDouble(X_IN) ) {
dawn@0: 		mexErrMsgTxt("INTERP_INPLACE: X should be doubles");
dawn@0: 	}
dawn@0: 	if (!mxIsNumeric(IDX_IN) || mxIsComplex(IDX_IN) || 
dawn@0: 		mxIsSparse(IDX_IN)	|| !mxIsDouble(IDX_IN) ) {
dawn@0: 		mexErrMsgTxt("INTERP_INPLACE: Y should be doubles");
dawn@0: 	}
dawn@0: 	mx=mxGetM(X_IN);		/* rows */
dawn@0: 	nx=mxGetN(X_IN);		/* columns */
dawn@0: 	midx=mxGetM(IDX_IN);		
dawn@0: 	nidx=mxGetN(IDX_IN);		
dawn@0: 
dawn@0: 	if (nx>1 || nidx>1) {
dawn@0: 		mexErrMsgTxt("INTERP_INPLACE: X and IDX should be column vectors");
dawn@0: 	}
dawn@0: 	if (mx<1) {
dawn@0: 		mexErrMsgTxt("INTERP_INPLACE: X should have at least two samples");
dawn@0: 	}
dawn@0: 		
dawn@0: 	xp = mxGetPr(X_IN);
dawn@0: 	idxp = mxGetPr(IDX_IN);
dawn@0: 	
dawn@0: 	checkin_matrix((mxArray *) IDX_IN);
dawn@0: 	checkin_matrix((mxArray *) X_IN);
dawn@0: 
dawn@0: 	/* Do the actual computations in a subroutine */
dawn@0: 	interp_inplace(xp,idxp,mx,midx);
dawn@0: 
dawn@0: 	checkout_matrix((mxArray *) X_IN);
dawn@0: 	checkout_matrix((mxArray *) IDX_IN);
dawn@0: 	return;
dawn@0: }
dawn@0: