dawn@0: /* 
dawn@0: * minInRange.c
dawn@0: * for each sample, return index of the smallest value within an interval 
dawn@0: * of that sample
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: 
dawn@0: #define	X_IN	       prhs[0]
dawn@0: #define	INTERVAL_IN	   prhs[1]
dawn@0: 
dawn@0: /* Output Arguments */
dawn@0: 
dawn@0: #define	IDX_OUT	plhs[0]
dawn@0: 
dawn@0: static void mininrange(
dawn@0:             double *xp,
dawn@0:             double *intp,
dawn@0:             double *idxp,
dawn@0:             unsigned int n
dawn@0:             )
dawn@0: {
dawn@0:     int h,k, idx, left, right, interval;
dawn@0:     double min;
dawn@0: 	double max;
dawn@0:     
dawn@0: 	
dawn@0:     for (k=0; k<n; k++) {
dawn@0:     	interval=(int) GET(intp[k]);
dawn@0: 
dawn@0:         left = k - interval/2;
dawn@0:         right = left+interval;  
dawn@0:         
dawn@0: 		if (left<0) left=0; 	/* clip */
dawn@0: 		if (right>n) right=n;
dawn@0: 			
dawn@0:         min = GET(xp[k]);
dawn@0:         idx = k;
dawn@0:         for (h=left;h<right; h++) {
dawn@0:             if (GET(xp[h]) < min) {  /* update min */
dawn@0:                 min = GET(xp[h]);
dawn@0:                 idx = h;
dawn@0:             }
dawn@0:         }
dawn@0:         SET(idxp[k])=idx+1;
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, *intp;
dawn@0:   unsigned int	m,n;
dawn@0:   
dawn@0:   /* Check for proper number of arguments */
dawn@0:   
dawn@0:   if (nrhs != 2) {
dawn@0:     mexErrMsgTxt("MININRANGE requires two input arguments");
dawn@0:   } else if (nlhs != 1) {
dawn@0:     mexErrMsgTxt("MININRANGE requires one output argument");
dawn@0:   }
dawn@0:   
dawn@0:   /* Check type of input */
dawn@0:   
dawn@0:   if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || 
dawn@0:       mxIsSparse(X_IN)  || !mxIsDouble(X_IN)) {
dawn@0:     mexErrMsgTxt("MININRANGE: X should be matrix of doubles");
dawn@0:   }
dawn@0:   if (!mxIsNumeric(INTERVAL_IN) || mxIsComplex(INTERVAL_IN) || 
dawn@0:       mxIsSparse(INTERVAL_IN)  || !mxIsDouble(INTERVAL_IN)) {
dawn@0:     mexErrMsgTxt("MININRANGE X should be matrix of doubles");
dawn@0:   }
dawn@0: 
dawn@0:   m = mxGetM(X_IN); /* rows */
dawn@0:   n = mxGetN(X_IN); /* columns */ 
dawn@0:   if (m>1 || n <=1) {
dawn@0:       mexErrMsgTxt("MININRANGE X should be row vector");
dawn@0:   }
dawn@0:   if (m != mxGetM(INTERVAL_IN) || n != mxGetN(INTERVAL_IN)) {
dawn@0: 	  mexErrMsgTxt("MININRANGE: INTERVAL should be of same size as X");
dawn@0:   }
dawn@0:   
dawn@0:   /* Create matrix to return */
dawn@0:   
dawn@0:   IDX_OUT = mxCreateDoubleMatrix(1, n, mxREAL);  
dawn@0:   
dawn@0:   /* Assign pointers to the various parameters */
dawn@0:   
dawn@0:   xp = mxGetPr(X_IN);
dawn@0:   intp = mxGetPr(INTERVAL_IN);
dawn@0:   idxp = mxGetPr(IDX_OUT);
dawn@0:   
dawn@0:   checkin_matrix((mxArray *) X_IN);
dawn@0:   checkin_matrix((mxArray *) INTERVAL_IN);
dawn@0:   checkin_matrix(IDX_OUT);
dawn@0:   
dawn@0:   /* Do the actual computations in a subroutine */
dawn@0:   
dawn@0:   mininrange(xp,intp,idxp,n);
dawn@0: 
dawn@0:   checkout_matrix((mxArray *) X_IN);
dawn@0:   checkout_matrix((mxArray *) INTERVAL_IN);
dawn@0:   checkout_matrix(IDX_OUT);
dawn@0:   return;
dawn@0: }
dawn@0: 
dawn@0: