dawn@0: /* 
dawn@0: * cumnorm_inplace.c
dawn@0: * cumulative mean-normalization of diff function
dawn@0: * 
dawn@0: * Alain de Cheveigné, CNRS/Ircam	Dec 2001
dawn@0: * (c) 2001 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: 
dawn@0: /* Output Arguments */
dawn@0: 
dawn@0: static void cumsum_inplace(
dawn@0: 				double *xp,		/* matrix to cumulative-mean-normalize (along columns) */
dawn@0: 				int m, 			/* rows */
dawn@0: 				int n	 		/* columns */ 
dawn@0: 			) 
dawn@0: {
dawn@0: 	int j,k;
dawn@0: 	double z, sum, mean;
dawn@0: 		
dawn@0: 	
dawn@0: 	for (k=0; k<n; k++) {							/* for each column */
dawn@0: 		SET(xp[k*m+0])=1;
dawn@0: 		sum = 0.0;
dawn@0: 		for (j=1;j<m;j++) {							/* for each row */
dawn@0: 			z = GET(xp[k*m+j]);
dawn@0: 			z = (z > 0) ? z:0;						/* clip to remove numerical artifacts */
dawn@0: 			sum += z;
dawn@0: 			z = (sum>0) ? (z / (sum/j)) : 1;		/* cumulative-mean-normalize */
dawn@0: 			SET(xp[k*m+j])=z;
dawn@0: 		}
dawn@0: 	}
dawn@0: /*
dawn@0: 	for (j=0;j<m;j++) {									// for each row (time) 
dawn@0: 		sum=0;
dawn@0: 		SET(xp[j])=1;
dawn@0: 		for (k=1; k<n; k++) {							// for each column (lag)
dawn@0: 			z=GET(xp[m*k+j]);
dawn@0: 			sum+=z;
dawn@0: 			mean=sum/k;
dawn@0: 			SET(xp[m*k+j]) = (mean > 0) ? z/mean : 1;
dawn@0: 		}
dawn@0: 	}
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; 
dawn@0: 	int nx, mx;
dawn@0: 	
dawn@0: 	/* Check for proper number of arguments */
dawn@0: 	if (nrhs != 1) {
dawn@0: 		mexErrMsgTxt("CUMSUM_INPLACE takes 1 input argument");
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("CUMSUM_INPLACE: X should be doubles");
dawn@0: 	}
dawn@0: 	mx=mxGetM(X_IN);		/* rows */
dawn@0: 	nx=mxGetN(X_IN);		/* columns */
dawn@0: 	if (nx*mx == 0) {
dawn@0: 		mexErrMsgTxt("CUMSUM_INPLACE: input matrix is empty");
dawn@0: 	}	
dawn@0: 		
dawn@0: 	xp = mxGetPr(X_IN);	
dawn@0: 	checkin_matrix((mxArray *) X_IN);
dawn@0: 			
dawn@0: 	
dawn@0: 	/* Do the actual computations in a subroutine */
dawn@0: 	cumsum_inplace(xp,mx,nx);
dawn@0: 
dawn@0: 	checkout_matrix((mxArray *) X_IN);
dawn@0: 	return;
dawn@0: }
dawn@0: