/*#include <string.h>*/
#include "mex.h"
#include <math.h>

void locatePeak(double *stream, mwSize length_stream, mwSize *p_start, mwSize *p_extr, mwSize *p_end, double *p_ampl);
/*Locates next peak in stream after index p_start */
/*Input: */
/*stream: the stream in which peaks are to be found */
/*length_stream: the length of the stream (in datapoints) */
/*p_start: from that index on the next peak will be detected in stream */
/*Output: */
/*p_start: the beginning of the peak i.e. the index of the closest local minimum left to it */
/*p_extr: the index of the peak */
/*p_end: the end of the peak i.e. the index of the closest local minimum right to it */
/*p_ampl: the amplitude of the peak */

void add_to_array_mwSize(mwSize **array, mwSize *size_of_array, mwSize nbr_of_elemnts, mwSize element);
/*Adds the element "element" to the mwSize-array "array". */
/*Input: */
/*array: the array  */
/*size_of_array: Number of elements that can fit in the array. If this number is exceeded that array size */
/*is automatically doubled and the array and size_of_array values are updated. */
/*nbr_of_elements: the current number of used elements in the array */
/*element: the element to be added */

void add_to_array_double(double **array, mwSize *size_of_array, mwSize nbr_of_elements, double element);
/*Adds the element "element" to the double-array "array". */
/*Input: */
/*array: the array  */
/*size_of_array: Number of elements that can fit in the array. If this number is exceeded that array size */
/*is automatically doubled and the array and size_of_array values are updated. */
/*nbr_of_elements: the current number of used elements in the array */
/*element: the element to be added */

void rm_from_array_mwSize(mwSize *array, mwSize nbr_of_elements, mwSize element_idx);
/*Removes the element_idx-th element from the mwSize-array "array" */
/*Input: */
/*array: the array from which an element should be removed */
/*nbr_of_elements: the number of elements the array consists of */
/*element_idx: the index of the element to be removed */

void rm_from_array_double(double *array, mwSize nbr_of_elements, mwSize element_idx);
/*Removes the element_idx-th element from the double-array "array" */
/*Input: */
/*array: the array from which an element should be removed */
/*nbr_of_elements: the number of elements the array consists of */
/*element_idx: the index of the element to be removed */

void rm_multiplets(mwSize **start_array, mwSize **extr_array, mwSize **end_array, double **ampl_array, mwSize *nbr_of_peaks, mwSize min_dist, int idx_policy, int ampl_policy);
/*Removes peaks which are close together ("multiplets") */
/*Input: */
/*start_array, extr_array, end_array, ampl_array: arrays describing the peaks (start_array: indices of the closest left local minima of the peaks, */
/*extr_array: indices of the peaks in the stream, end_array: indices of the closest right local minima of the peaks, ampl_array: the amplitudes of the peaks) */
/*nbr_of_peaks: the number of peaks detected so far */
/*min_dist: two peaks are seen as a multiplet and coalesced if they are separated by less than min_dist data points */
/*idx_policy: defines at what position a new peak representing several coalesced peaks should appear */
/*idx_policy==0: new peak at position of leftmost peak;  */
/*idx_policy==1: new peak at position of rightmost peak;  */
/*idx_policy==2: new peak in the middle of interval spanned by the peaks to be coalesced */
/*idx_policy==3; new peak at the position of the highest peak */
/*ampl_policy: defines what amplitude a new peak which is representing several coalesced peaks should have */
/*ampl_policy==0: new peak with same amplitude as leftmost peak */
/*ampl_policy==1: new peak with same amplitude as rightmost peak */
/*ampl_policy==2: new peak with average amplitude of the peaks to be replaced */
/*ampl_policy==3: new peak with the amplitude of the highest peak */

mwSize max_idx_double(double* ampl_array, mwSize nbr_of_elements);
/*Returns the index of the sample with maximal amplitude of an array */
/*Input: */
/*ampl_array: the array containging the amplitude information */
/*nbr_of_elements: the length of the array */
/*Output: */
/*the index of the sample with maximal amplitude */

mwSize max_idx_mwSize(mwSize* ampl_array, mwSize nbr_of_elements);
/*Returns the index of the sample with maximal amplitude of an array */
/*Input: */
/*ampl_array: the array containging the amplitude information */
/*nbr_of_elements: the length of the array */
/*Output: */
/*the index of the sample with maximal amplitude */

double average(double* ampl_array, mwSize nbr_of_elements);
/*Computes the average of the elements of an array */
/*Input: */
/*ampl_array: array containing the samples to be averaged */
/*nbr_of_elements: the number of elements ampl_array consists of */
/*Output: */
/*The average of the values of ampl_array */

void height_wrt_minima(double* stream, double *ampl_array, mwSize *start_array, mwSize *end_array, double *left_height_array, double *right_height_array, mwSize nbr_of_peaks);
/*Computes the heights of peaks w.r.t to their closest left and right local minimum, respectively. */
/*Input: */
/*stream: the data stream */
/*ampl_array: the array containing the amplitudes of the peaks */
/*start_array: the array containing the indices of the closes left local minima of the individual peaks */
/*end_array: the array containing the indices of the closest right local minima of the individual peaks */
/*mwSize nbr_of_peaks: the number of peaks */
/*Output: */
/*left_height_array: array containing the heights of the peaks with respect to the closest left local minima */
/*right_height_array: array containing the heights of the peaks with respect to the closest right local minima */

void fwhm(double* stream, mwSize *start_array, mwSize *end_array, double* left_height_array, double* right_height_array, double *fwhm_array, mwSize nbr_of_peaks);
/*Computes a kind of pseude full width at half maximum of the peaks */
/*The half maximum of the k-th peak is computed as hm=0.5*min(rigth_height_array[k],left_height_array[k]); */
/*The fwhm is then computed with respect to hm. */
/*Input: */
/*stream: the recorded data stream */
/*left_height_array: the heights of the peaks with respect to the closest left minimum */
/*right_height_array: the heights of the peaks with respect to the closest right minimum */
/*nbr_of_peaks: the number of peaks */
/*Output: */
/*fwhm_array: array containing the pseude fwhms of the peaks */


void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
/*Matlab call: [start_array,ext_array,end_array,ampl_array,left_height_array,right_height_array,fwhm_array]=peakdetector(stream,min_dist,idx_policy,ampl_policy) */
{
	/*Stream realted: */
	mwSize length_stream, k;	/*Length of provided stream */
	double *stream;		/*array containing stream */
	
	/*Primary peak detection */
	mwSize p_start,p_extr,p_end; 	/*define peak */
									/*p_start: index of closest minima left of peak */
									/*p_extr: index of maximum */
									/*p_end: index of closest minima right of peak */
	double p_ampl;					/*p_ampl: the amplitude of the peak */
	mwSize *start_array = NULL;	    /*arrays collecting the individual peak informations */
    mwSize *extr_array = NULL;
    mwSize *end_array = NULL;
	double *ampl_array = NULL;		/*array containing the amplitudes of teh individual peaks */
    mwSize start_array_size = 0;    /*the number of elements allocated for the arraya above */
    mwSize extr_array_size = 0;
    mwSize end_array_size = 0;  
    mwSize ampl_array_size = 0;     
	mwSize nbr_of_peaks=0;			/*The number of peaks found in the stream */
	
    /*Amplitude filter */
    double min_ampl;
    
	/*Multiplet reduction related */
	int min_dist;	/*if two peaks are closer than this value they are coalesced */
	int idx_policy; /*policy indicating at which point a new peak representing some coalesced peaks should appear */
	int ampl_policy; /*policy indicating what amplitude a new peak representign some coalesced peaks should have */

	/*Peak analysis[start_array,ext_array,end_array,ampl_array,left_height_array,right_height_array,fwhm_array]=peak_detector(stream<,min_dist,idx_policy,ampl_policy>) */
	double *left_height_array, *right_height_array; /*arrays containing the heights of the peaks w.r.t. to the closest left, respectively right local minimum */
	double *fwhm_array; /*array containing the full width at half maximum of the peaks */

	/*Return data */
	double *z; /*auxiliary pointer needed to write double arrays to arrays of type mxArray */


	/*Check if right number of arguments was provided */
	if ( nrhs!=1 && nrhs!=2 && nrhs!=4 && nrhs!=5 )
	{
		mexErrMsgTxt("Invalid number of input arguments");
	}

	/*Get stream */
	length_stream=mxGetN(prhs[0]);
	stream=mxGetPr(prhs[0]);

	/*Get input parameters */
	if (nrhs==2)
	{
       min_ampl=(double)(*(mxGetPr(prhs[1])));
    }
    else if(nrhs==4)
    {
		min_dist=(int)(*(mxGetPr(prhs[1])));
		idx_policy=(int)(*(mxGetPr(prhs[2])));
		ampl_policy=(int)(*(mxGetPr(prhs[3])));
	}
    else if(nrhs==5)
    {
        min_ampl=(double)(*(mxGetPr(prhs[1])));
      	min_dist=(int)(*(mxGetPr(prhs[2])));
		idx_policy=(int)(*(mxGetPr(prhs[3])));
		ampl_policy=(int)(*(mxGetPr(prhs[4])));
    }

	p_start=0;	/*Start to search for peaks from the beginning of stream */
	while (p_start+1 < length_stream)
	{
			locatePeak(stream,length_stream,&p_start,&p_extr,&p_end,&p_ampl);
			/*Add peak info to corresponding arrays if it satisfies the min_map criterion */
            if( (nrhs==1 || nrhs==4) || ((nrhs==2 || nrhs==5) && p_ampl > min_ampl) )
            {
                add_to_array_mwSize(&start_array,&start_array_size,nbr_of_peaks,p_start);
                add_to_array_mwSize(&extr_array,&extr_array_size,nbr_of_peaks,p_extr);
                add_to_array_mwSize(&end_array,&end_array_size,nbr_of_peaks,p_end);
                add_to_array_double(&ampl_array,&ampl_array_size,nbr_of_peaks,p_ampl);
                nbr_of_peaks++;
            }
            p_start=p_end;		
	}

	/*Check if last local maximum was indeed a peak (and not only the last sample at the end of the stream) */
	if (nbr_of_peaks > 0 && p_extr==p_end) 
/* 	if (nbr_of_peaks > 0 && p_start==p_end)*/
	{
		nbr_of_peaks--;
	}

	/*Remove multiplets */
	if (nrhs>=4)
	{
		rm_multiplets(&start_array,&extr_array,&end_array,&ampl_array,&nbr_of_peaks,min_dist,idx_policy,ampl_policy);
	}
    
	/*Determine left and right heights */
    if (nlhs >= 5)
    {
        /*Compute heights of peaks w.r.t. closest left and right local minima */
    	/*Allocate memory */
        left_height_array=mxMalloc(nbr_of_peaks*sizeof(double));
    	right_height_array=mxMalloc(nbr_of_peaks*sizeof(double));
    	height_wrt_minima(stream, ampl_array, start_array, end_array, left_height_array, right_height_array, nbr_of_peaks);
    }

	/*Compute fwhm of the peaks */
    if (nlhs >= 7)
    {
    	fwhm_array=mxMalloc(nbr_of_peaks*sizeof(double));
        fwhm(stream, start_array, end_array, left_height_array, right_height_array, fwhm_array, nbr_of_peaks);
    }

	/*Return results to MATLAB */
	if (nlhs>=1)
	{
		plhs[0]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[0]);
		for (k=0; k<nbr_of_peaks; k++)
		{
			z[k]=(double)start_array[k]+1;
		}
	}
	
	if (nlhs>=2)
	{
		plhs[1]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[1]);
		for (k=0; k<nbr_of_peaks; k++)
		{
			z[k]=(double)extr_array[k]+1;
		}
	}
	
	if (nlhs>=3)
	{
		plhs[2]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[2]);
		for (k=0; k<nbr_of_peaks; k++)
		{
			z[k]=(double)end_array[k]+1;
		}
	}

	if (nlhs>=4)
	{
		plhs[3]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[3]);
        memcpy(z, ampl_array, nbr_of_peaks*sizeof(double));
	}

	if (nlhs>=5)
	{
		plhs[4]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[4]);
        memcpy(z, left_height_array, nbr_of_peaks*sizeof(double));
	}

	if (nlhs>=6)
	{
		plhs[5]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[5]);
        memcpy(z, right_height_array, nbr_of_peaks*sizeof(double));	}

	if (nlhs>=7)
	{
		plhs[6]=mxCreateDoubleMatrix(1,nbr_of_peaks,mxREAL);
		z=mxGetPr(plhs[6]);
        memcpy(z, fwhm_array, nbr_of_peaks*sizeof(double));
	}

/*	mexCallMATLAB(0,NULL,1,prhs,"plot"); */
	
/*	mexPrintf("free memory\n"); */
/*	mexCallMATLAB(0,NULL,0,NULL,"pause"); */

	/*Free memory */
    if(nbr_of_peaks>0) {
        mxFree(start_array);
        start_array=NULL;
        mxFree(extr_array);
        extr_array=NULL;
        mxFree(end_array);
        end_array=NULL;
        mxFree(ampl_array);
        ampl_array=NULL;
        
        if (nlhs >=5)
        {
            mxFree(left_height_array);
            left_height_array=NULL;
            mxFree(right_height_array);
            right_height_array=NULL;
        }
        
        if (nlhs >= 7)
        {
            mxFree(fwhm_array);
            fwhm_array=NULL;
        }
    }
}


void locatePeak(double *stream, mwSize length_stream, mwSize *p_start, mwSize *p_extr, mwSize *p_end, double *p_ampl)
{
    /*If we discover a plateau while searching for the end of the peak, we store the  */
    /*first idx of the plateau in this variable, so we can later decide if this */
    /*point marks the end of the peak or not. If the signal starts to climb at */
    /*the other end of the plateau, this was the end of the current peak. If the */
    /*signal continues to fall after the end of the plateau, we keep searching  */
    /*for the end of the peak. */
	mwSize tmp_end = 0;
    
	/*Search for first local minimum right to *p_start */
	while ((*p_start)+1 < length_stream)
	{
		if (stream[(*p_start)+1]-stream[*p_start] > 0)
		{		
			break;
		}
		else
		{
			(*p_start)++;
		}
	}
	
	/*Search for extremum (actual peak) */
	*p_extr=*p_start;
	while((*p_extr)+1 < length_stream)
	{
		if (stream[(*p_extr)+1]-stream[*p_extr]<0)
		{
			break;
		}
		else
		{
			(*p_extr)++;
		}
	}

	/*Set amplitude */
	*p_ampl=stream[*p_extr];

	/*Search for end of peak */
	*p_end=*p_extr;
	while ((*p_end)+1 < length_stream)
	{
        double diff = stream[(*p_end)+1]-stream[*p_end];
		if (diff > 0.0)
		{
            if (tmp_end != 0)
            {
                *p_end = tmp_end;
            }
			break;
		}
        else
		{
            if (diff == 0.0 && tmp_end == 0)
            {
                tmp_end = *p_end;
            }
            else if(diff < 0.0 && tmp_end != 0)
            {
                tmp_end = 0;
            }
			(*p_end)++;
		}
	}
}

void add_to_array_mwSize(mwSize **array, mwSize *array_size, mwSize nbr_of_elements, mwSize element)
{
	/*mwSize k; */
        
    /*If there is no more room in the allocated array: allocate a */
    /*new array which is twice as long */
    if(*array_size < nbr_of_elements+1)
    {
        mwSize *new_array;

        /*If this is the first call allocate an inital array with room for 100 elements */
        if(*array_size == 0)
        {
            *array_size = 128;
        }
        else
        {
            *array_size *= 2;
        }
        new_array = mxMalloc(*array_size * sizeof(mwSize));
        
        /*Copy old values to new array */
        memcpy(new_array, *array, nbr_of_elements*sizeof(mwSize));
        
        /*Free old array */
        if (nbr_of_elements!=0)
        {
            mxFree(*array);
        }
        
        /*Use the new array from now on */
        *array = new_array;
    }
    
    /*Add new element */
	(*array)[nbr_of_elements] = element;
}

void add_to_array_double(double **array, mwSize *array_size, mwSize nbr_of_elements, double element)
{
	/*mwSize k; */
    
    /*If there is no more room in the allocated array: allocate a */
    /*new array which is twice as long */
    if(*array_size < nbr_of_elements+1)
    {
        double *new_array;
        
        /*If this is the first call allocate an inital array with room for 100 elements */
        if(*array_size == 0)
        {
            *array_size = 128;
        }
        else
        {
            *array_size *= 2;
        }
        new_array = mxMalloc(*array_size * sizeof(double));
        
        /*Copy old values to new array */
        memcpy(new_array, *array, nbr_of_elements*sizeof(double));        
        
        /*Free old array */
        if (nbr_of_elements!=0)
        {
            mxFree(*array);
        }
        
        /*Use the new array from now on */
        *array = new_array;
    }
    
    /*Add new element */
	(*array)[nbr_of_elements] = element;
}

void rm_from_array_mwSize(mwSize *array, mwSize nbr_of_elements, mwSize element_idx)
{
	mwSize k;
	
	/*Check if element_idx is in [0 nbr_of_elements-1] */
	if (element_idx < 0 || element_idx >= nbr_of_elements)
	{
		mexErrMsgTxt("Element to be removed out of array range!");
	}
	
	/*move elements to their new positions */
	for (k=element_idx+1; k<nbr_of_elements; k++)
	{
		array[k-1]=array[k];
	}
}

void rm_from_array_double(double *array, mwSize nbr_of_elements, mwSize element_idx)
{
	mwSize k;
	
	/*Check if element_idx is in [0 nbr_of_elements-1] */
	if (element_idx < 0 || element_idx >= nbr_of_elements)
	{
		mexErrMsgTxt("Element to be removed out of array range!");
	}
	
	/*move elements to their new positions */
	for (k=element_idx+1; k<nbr_of_elements; k++)
	{
		array[k-1]=array[k];
	}
}

void rm_multiplets(mwSize **start_array, mwSize **extr_array, mwSize **end_array, double **ampl_array, mwSize *nbr_of_peaks, mwSize min_dist, int idx_policy, int ampl_policy)
{
	mwSize l;
	mwSize k=0;
	/*Auxiliary counter */

	mwSize distance;
	/*distance between two consecutive peaks */

	mwSize peakseg_start, peakseg_end;
	/*Indicating the indices of the peaks to be coalesced */

	/*Loop over peak list */
	while (k+1<(*nbr_of_peaks))
	{
		peakseg_start=k;
		peakseg_end=k;
		distance=(*extr_array)[peakseg_start+1]-(*extr_array)[peakseg_start];

        /*Check if consecutive peaks are closer than min_dist. */
		while (distance < min_dist)
		{
			k++;
			peakseg_end=k;
			if (k+1 == *nbr_of_peaks)
			{
				break;
			}
			distance=(*extr_array)[k+1]-(*extr_array)[k];
        }
         
        
		/*Check if a multiplet was found. If so coalesce it. Info about new peak will be stored in element peakseg_start of the */
		/*corresponding arrays */
		if (peakseg_end != peakseg_start)
		{
				/*Set new peak and peak range	 */
				/*set new peak range */
				(*end_array)[peakseg_start]=(*end_array)[peakseg_end];

				switch (idx_policy)
				{
					/*case 0: Set new peak to leftmost position. Nothing to do. New peak position equals position of the leftmost peak.  */
					case 1:	/*Set new peak to rightmost position */
						(*extr_array)[peakseg_start]=(*extr_array)[peakseg_end];
						break;
					case 2:	/*Set new peak in middle of all peaks */
						(*extr_array)[peakseg_start]= (mwSize)((((double)(*extr_array)[peakseg_start]+(double)(*extr_array)[peakseg_end])/2)+0.5);
						break;
					case 3: /*Set new peak at the position of the highest peak */
                        (*extr_array)[peakseg_start]=(*extr_array)[peakseg_start+max_idx_double(*ampl_array+peakseg_start,peakseg_end-peakseg_start+1)];
						break;
				}

				/*Set peak amplitude */
				switch (ampl_policy)
				{
					/*case 0: set peak amplitude to that of leftmost peak. Nothing to do. New peak amplitude equals amplitude of leftmost peak */
					case 1:	/*Set new peak's amplitude to that one of the rightmost peak */
						(*ampl_array)[peakseg_start]=(*ampl_array)[peakseg_end];
						break;
					case 2: /*Set new peak's amplitude to average of detected peaks */
						(*ampl_array)[peakseg_start]=average(*ampl_array+peakseg_start,peakseg_end-peakseg_start+1);
						break;
					case 3: /*Set new peaks's amplitude to the one of the largest peak */
						(*ampl_array)[peakseg_start]=(*ampl_array)[peakseg_start+max_idx_double(*ampl_array+peakseg_start,peakseg_end-peakseg_start+1)];
						break;
				}

				/*Remove superfluous peaks */
				for (l=peakseg_end; l>peakseg_start; l--)
				{
					rm_from_array_mwSize(*start_array,*nbr_of_peaks,l);
					rm_from_array_mwSize(*extr_array,*nbr_of_peaks,l);
					rm_from_array_mwSize(*end_array,*nbr_of_peaks,l);
					rm_from_array_double(*ampl_array,*nbr_of_peaks,l);
					/*Adjust peaks counter */
					(*nbr_of_peaks)--;
                    k--;
				}
              
		}
		k++;
	}
}

mwSize max_idx_double(double* ampl_array, mwSize nbr_of_elements)
{
	mwSize max_idx=0;	/*contains the index of the sample with maximal amplitude */
	double max_val;		/*auxiliary variable used to compare the amplitudes of the samples */
	mwSize k;			/*auxiliary counter */

	/*Determine sample with maximal amplitude */
	max_val=ampl_array[max_idx];
	for (k=0; k<nbr_of_elements; k++)
	{
		if(ampl_array[k]>max_val)
		{
			max_val=ampl_array[k];
			max_idx=k;
		}
	}
	return(max_idx);
}

mwSize max_idx_mwSize(mwSize* ampl_array, mwSize nbr_of_elements)
{
	mwSize max_idx=0;	/*contains the index of the sample with maximal amplitude */
	mwSize max_val;		/*auxiliary variable used to compare the amplitudes of the samples */
	mwSize k;			/*auxiliary counter */

	/*Determine sample with maximal amplitude */
	max_val=ampl_array[max_idx];
	for (k=0; k<nbr_of_elements; k++)
	{
		if(ampl_array[k]>max_val)
		{
			max_val=ampl_array[k];
			max_idx=k;
		}
	}
	return(max_idx);
}

double average(double* ampl_array, mwSize nbr_of_elements)
{
	mwSize k; /*auxiliary counter */
	double average=0;	/*average of samples */

	/*Compute average */
	for (k=0; k<nbr_of_elements; k++)
	{
		average+=ampl_array[k];
	}
	average/=(double)(nbr_of_elements);
	
	return(average);
}

void height_wrt_minima(double *stream, double *ampl_array, mwSize *start_array, mwSize *end_array, double *left_height_array, double *right_height_array, mwSize nbr_of_peaks)
{
	mwSize k; /*auxiliary counter */
	
	for (k=0; k<nbr_of_peaks; k++)
	{
		left_height_array[k]=ampl_array[k]-stream[start_array[k]];
		right_height_array[k]=ampl_array[k]-stream[end_array[k]];
	}
}

void fwhm(double* stream, mwSize *start_array, mwSize *end_array, double* left_height_array, double* right_height_array, double *fwhm_array, mwSize nbr_of_peaks)
{
	mwSize k; /*auxiliary counter */
	double hm; /* the half maximum of the individual peaks */
	mwSize p_left, p_right; /*pointers indicating at which point left and right of the peak, respectively, the half maximum is reached */

	double debug;

	/*Loop over peaks */
	for (k=0; k<nbr_of_peaks; k++)
	{	
		/*Determine half maximum */
		if (left_height_array[k]<right_height_array[k])
		{
			hm=0.5*left_height_array[k]+stream[start_array[k]];
		}
		else
		{
			hm=0.5*right_height_array[k]+stream[end_array[k]];
		}

		/*Search for point left of peak at which half maximum is reached */
		p_left=start_array[k];
		while (stream[p_left]<hm)
		{
			debug=stream[p_left];
			p_left++;
		}
		/*Search for point left of peak at which half maximum is reached */
		p_right=end_array[k];
		while (stream[p_right]<hm)
		{
			p_right--;
		}

		/*Compute fwhm */
		fwhm_array[k]=p_right-p_left;
	}
}
