#include <stdio.h>
#include <math.h>
#include "mex.h"

/*******************************************************/
/* Puschel_stage                                       */
/* A nonlinear adaptation filter (Dau, 1996)           */
/*       By Aki Härmä 26. 11 1998  		       */
/*******************************************************/
void trans(double *sig, long int len, double b, double yy, 
	   double *RET, double *YY)
{
 long int q, w;
 double x, cc;

cc=1-b;

 for (q=0;q<len;q++)
   {
     x=sig[q]/yy; 
     RET[q]=x;
     yy=cc*x+b*yy;
   }
YY[0]=yy;
return;
}

/************************  MEXFUNCTION *****************************/
void mexFunction(
	int		nlhs,
	mxArray	*plhs[],
	int		nrhs,
	const mxArray    *prhs[]
	)
{
  long int mm,nn,len;
  double *sig, *RET, *YY;
  double b, yy;
 
        sig = mxGetPr(prhs[0]);
      	mm = (long int)mxGetM(prhs[0]);
	nn= (long int)mxGetN(prhs[0]);
	if(nn==1)len=mm; else len=nn;

	b = mxGetScalar(prhs[1]);
	yy= mxGetScalar(prhs[2]);

       plhs[0] = mxCreateDoubleMatrix(1, len, mxREAL);
	RET = mxGetPr(plhs[0]);
       plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
	YY = mxGetPr(plhs[1]);

	trans(sig,len,b,yy,RET,YY);
	
	return;
}