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

double  *NFArray(int size)
{
	double	*p;
	p = (double *)mxCalloc(sizeof(*p),size);
	return p;
}

/*******************************************************/
/* Puschel_stage                                       */
/* A multi-channel implementation of
/* a nonlinear adaptation filter (Dau, 1996)           */
/*       By Aki Härmä 30. 11 1998  		       */
/*******************************************************/
void trans(double *sig, long int len, long int dim, double b,
	   double *RET)
{
 long int q, w, ind;
 int m,n;
 double x, cc;
 double *yy,*y;
 yy=NFArray(dim); y=NFArray(dim);
 for (q=0;q<dim;q++){y[q]=1.0;yy[q]=1.1;}
 cc=1-b;

 for (q=0;q<len;q++)
   {
   for (w=0;w<dim;w++)
     {
     ind=q+w*len;
     x=sig[ind]/y[w]; 
     RET[ind]=x;
     y[w]=cc*x+b*y[w];
     }
   }
return;
}

/************************  MEXFUNCTION *****************************/
void mexFunction(
	int		nlhs,
	mxArray	*plhs[],
	int		nrhs,
	const mxArray    *prhs[]
	)
{
  long int len, dim, m, n;
  double *sig, *RET;
  double b, yy;
 
        sig = mxGetPr(prhs[0]);
      	len = (long int)mxGetM(prhs[0]);
	dim = (long int)mxGetN(prhs[0]);
	/*mexPrintf("Dim: %d   Len: %d",dim,len);*/
	if(len<dim)mexPrintf("More channels than samples?");

	b = mxGetScalar(prhs[1]);

       plhs[0] = mxCreateDoubleMatrix(len, dim, mxREAL);
	RET = mxGetPr(plhs[0]);
	
	if(len<dim)mexPrintf("Puschel: More channels than samples?");
	else trans(sig,len,dim,b,RET);
	
	return;
}