dawn@0: /* 
dawn@0: * minInterp.c
dawn@0: * determine position of minimum of parabolic fit to local minima
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: 
dawn@0: /* Output Arguments */
dawn@0: #define	MIN_OUT	plhs[0]
dawn@0: #define	IND_OUT	plhs[1]
dawn@0: 
dawn@0: 
dawn@0: static void mininterp(
dawn@0:             double *xp,
dawn@0:             unsigned int m,
dawn@0:             unsigned int n,
dawn@0:             double *minp,
dawn@0:             double *indp
dawn@0: 		   )
dawn@0: {
dawn@0:     double	x1,x2,x3,a,b,min,ind;
dawn@0:     int i,j,k;
dawn@0:   
dawn@0:   	/* each column is interpolated independently */
dawn@0:   	/* each local minimum is replaced by the value of the interpolated minimum, and the 
dawn@0:   	position of the interplolated minimum is recorded in indp */
dawn@0:   	
dawn@0:     for (j=0; j<n; j++) { 
dawn@0:         /* first row: don't touch */
dawn@0:         SET(indp[j*m])=0;         
dawn@0:         SET(minp[j*m])=GET(xp[j*m]);
dawn@0:         /* intermediate: refine */
dawn@0:         for (k=1; k<m-1; k++) {   
dawn@0:             i=k+j*m;
dawn@0: 			x1=GET(xp[i-1]); x2=GET(xp[i]); x3=GET(xp[i+1]);  
dawn@0: 	        SET(indp[i]) = 0;	/* default */
dawn@0: 	        SET(minp[i]) = x2; 
dawn@0: 		    if ((x2<=x1) && (x2<=x3)) { 
dawn@0: 			    a = (x1 + x3 - 2*x2)/2;
dawn@0: 			    b = (x3 - x1)/2;
dawn@0: 			    if (a) {                     
dawn@0: 	                SET(indp[i]) = - b / (2*a);
dawn@0: 	                SET(minp[i]) = x2 +  - (b*b) / (4*a);
dawn@0: 	            }
dawn@0: 	        }
dawn@0: 	   	}
dawn@0:         /* last row: don't touch */
dawn@0:         SET(indp[m-1 + j*m])=0;
dawn@0:         SET(minp[m-1 + j*m])=GET(xp[m-1 + j*m]);
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	*minp, *indp, *xp;
dawn@0:   unsigned int	m,n;
dawn@0:   
dawn@0:   /* Check for proper number of arguments */
dawn@0:   
dawn@0:   if (nrhs != 1) {
dawn@0:     mexErrMsgTxt("minParabolic requires one input argument.");
dawn@0:   } else if (nlhs > 2) {
dawn@0:     mexErrMsgTxt("minParabolic requires one or two output arguments");
dawn@0:   }
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("minParabolic requires that X be a 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<3) { mexErrMsgTxt("number of rows should be 3 or more"); }
dawn@0:   
dawn@0:   
dawn@0:   /* Create two matrices to return */
dawn@0:   
dawn@0:   MIN_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
dawn@0:   IND_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
dawn@0:   
dawn@0:   
dawn@0:   /* Assign pointers to the various parameters */
dawn@0:   
dawn@0:   minp = mxGetPr(MIN_OUT);
dawn@0:   indp = mxGetPr(IND_OUT);
dawn@0:   xp = mxGetPr(X_IN);  
dawn@0:   
dawn@0:   checkin_matrix((mxArray *) X_IN);
dawn@0:   checkin_matrix(MIN_OUT);
dawn@0:   checkin_matrix(IND_OUT);
dawn@0:   
dawn@0:   /* Do the actual computations in a subroutine */
dawn@0:   
dawn@0: 
dawn@0:   mininterp(xp,m,n,minp,indp);
dawn@0: 
dawn@0:   checkout_matrix((mxArray *) X_IN);
dawn@0:   checkout_matrix(MIN_OUT);
dawn@0:   checkout_matrix(IND_OUT);
dawn@0:   return;
dawn@0: }
dawn@0: 
dawn@0: