/* localmax.c
 *
 * Written by Tracy Harton 11/95
 *
 * This function finds the local maximums of a set of data
 *
 */


#include <malloc.h>
#include <memory.h>
#include <stdio.h>

#include "mex.h"

void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
  int mb, nb;
  double *ap, *bp, *cp, *dp;
  int length;
  int i;
  double pos[100], max[100];
  
  if ((nrhs != 2)||(nlhs != 2))
    {
      mexErrMsgTxt("Usage: [pos,max]=peak(width, data)");
    }
  
  mb = mxGetM(prhs[1]);
  nb = mxGetN(prhs[1]);

  ap = mxGetPr(prhs[0]);
  bp = mxGetPr(prhs[1]);

  length = localmax(nb, *ap, bp, &pos, &max);
 
  plhs[0] = mxCreateFull(length, 1, REAL);
  plhs[1] = mxCreateFull(length, 1, REAL);
  cp = mxGetPr(plhs[0]);
  dp = mxGetPr(plhs[1]);

  for(i = 0; i < length; i++)
  {
	cp[i] = pos[i];
	dp[i] = max[i];
  }
  
}

int localmax(int num, double thresh, double *data, double *pos, double *max)
{
  double lmax = 0.0;
  int lpos = 0;
  int flag = 0;
  int i;
  int width_loop;
  int peak_num = 0;
  int above = 0;

  printf("Looking for peaks of thresh %f in %d datapoints\n", thresh, num);

  for(i = 0; i < num; i++)
  {
	flag = 0;

	if((data[i] > thresh) && !above)
	{
		flag = 1;
		above = 1;
	}

	if((data[i] < thresh) && above)
	{
		flag = 1;
		above = 0;
	}


/*
	if((data[i-1] > data[i]) || (data[i-2] > data[i-1]))
		flag = 0;

	if((data[i+1] > data[i]) || (data[i+2] > data[i+1]))
		flag = 0;

	if((data[i] - data[i-20]) < thresh)
		flag = 0;
*/

	if(flag==1)
	{
		printf("Peak: %d %f\n", i, data[i]);
		pos[peak_num] = (double) i;
		max[peak_num] = data[i];
		peak_num++;
	}

  }
  return(peak_num);
}