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

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

/*******************************************************/
/* AN_refractory91                                     */
/* A multi-channel implementation of a model of        */
/* refractory effects in Auditory Nerve from           */
/* (Meddis, 1991)                                      */
/*       By Aki Härmä 16. 12 1998  		       */
/*******************************************************/
void trans(double *sig, double fs, long int len, long int dim, 
	   double *OUT)
{
 long int t, j, ind;

 double *psum;
 double y;
 int dt;

 psum=NFArray(dim);

 dt=(int)fabs(0.001*fs); /* 1 ms refractory period */
 
 for (t=dt;t<len;t++)
   {
   for (j=0;j<dim;j++)
     {
     ind=t+j*len;
     y=sig[ind]*(1-psum[j]);
     psum[j]+=y-OUT[ind-dt];
     OUT[ind]=y;
     }
   }
return;
}

/************************  MEXFUNCTION *****************************/
void mexFunction(
	int		nlhs,
	mxArray	*plhs[],
	int		nrhs,
	const mxArray    *prhs[]
	)
{
  long int len, dim, m, n;
  double *sig, *OUT;
  double yy, fs;
  
        sig = mxGetPr(prhs[0]);
      	len = (long int)mxGetM(prhs[0]);
	dim = (long int)mxGetN(prhs[0]);

        fs = mxGetScalar(prhs[1]); /* Sampling rate */
	
	/*	mexPrintf("Dim: %d   Len: %d",dim,len);*/

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