tomwalters@0: /*
tomwalters@0:     Copyright (c) Applied Psychology Unit, Medical Research Council. 1993
tomwalters@0:     ===========================================================================
tomwalters@0: 
tomwalters@0:     Permission to use, copy, modify, and distribute this software without fee 
tomwalters@0:     is hereby granted for research purposes, provided that this copyright 
tomwalters@0:     notice appears in all copies and in all supporting documentation, and that 
tomwalters@0:     the software is not redistributed for any fee (except for a nominal 
tomwalters@0:     shipping charge). Anyone wanting to incorporate all or part of this 
tomwalters@0:     software in a commercial product must obtain a license from the Medical 
tomwalters@0:     Research Council.
tomwalters@0: 
tomwalters@0:     The MRC makes no representations about the suitability of this 
tomwalters@0:     software for any purpose.  It is provided "as is" without express or 
tomwalters@0:     implied warranty.
tomwalters@0:  
tomwalters@0:     THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
tomwalters@0:     ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL 
tomwalters@0:     THE A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES 
tomwalters@0:     OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, 
tomwalters@0:     WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 
tomwalters@0:     ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS 
tomwalters@0:     SOFTWARE.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*--------------------------------------------------------------------------*/
tomwalters@0: 
tomwalters@0: /*   findpeaks.c 
tomwalters@0: *   -------------
tomwalters@0: *
tomwalters@0: * Part of the temporal regularity programs
tomwalters@0: *
tomwalters@0: * M.Akeroyd. Summer 1993. Revised Winter 1994.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <string.h>
tomwalters@0: #include <stdlib.h>
tomwalters@0: 
tomwalters@0: #include "tip.h"
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*--------------------------------------------------------------------------*/
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int findpeaksreverse()
tomwalters@0: {
tomwalters@0:   /* What this code does: 
tomwalters@0:    * Its primary job is to find out where the tops of local maxima ("peaks") 
tomwalters@0:    * are. Its secondary job is to define the edges of those peaks.
tomwalters@0:    * The locations and heights of the peaks are held in the peak[n] structure,
tomwalters@0:    * where n is the number of the peak. The edges are held in the 
tomwalters@0:    * nextevent[n] structure. 
tomwalters@0:    *
tomwalters@0:    * Its called 'findpeaksreverse' because it traverses the array from 'right'
tomwalters@0:    * (ie, +ve Image Time) to 'left' (-ve Image Time): peaks get numbered
tomwalters@0:    * accordingly. eg, if nwidth = 0, then the TriggerPeak is peak 1.
tomwalters@0:    *
tomwalters@0:    * Edge defintions:
tomwalters@0:    *    for a peak n
tomwalters@0:    *        nextevent[n] is the LEFT edge (because we are going right-to-left
tomwalters@0:    *                                       on samples);
tomwalters@0:    *        nextevent[n-1] the RIGHT edge.
tomwalters@0:    * 
tomwalters@0:    * Peak definitions:
tomwalters@0:    *   simple peak:   f(sample - 1 ) <  f(sample)  > f(sample + 1)
tomwalters@0:    *   complex peak is a plateau: hopefully the code finds the middle, or 
tomwalters@0:    *                the leftedge of the middle (if the middle is an 
tomwalters@0:    *                even number of samples)
tomwalters@0:    *
tomwalters@0:    * Nextevents definitions: 
tomwalters@0:    *   minima:      f(sample - 1 ) > f(sample) < f(sample + 1)
tomwalters@0:    *   rightzero:   the first f(sample) = 0 to the Right of a peak
tomwalters@0:    *   leftzero:       ...       ...       ...     Left ...
tomwalters@0:    *                (so, because we're going backwards, the RIGHTZERO of a 
tomwalters@0:    *                 peak is actually held in nextevent[n-1], the LEFTZERO 
tomwalters@0:    *                 in nextevent[n]).
tomwalters@0:    *
tomwalters@0:    *
tomwalters@0:    * The first and last samples of a channel are dealt with differently, 
tomwalters@0:    * if only because the inequalities will crash if the (sample-1) or 
tomwalters@0:    * (sample+1) (respectively) don't exist. But I'm not sure if they work 
tomwalters@0:    * properly. I'm pretty certain they don't: so, they are just ignored.
tomwalters@0:    *
tomwalters@0:    * Finally, the number of peaks found is RETURNed.
tomwalters@0:    *
tomwalters@0:    * 
tomwalters@0:    * Version of 21 May 1993. Revised slighlty Winter 1994.
tomwalters@0:    */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*---------------------------*/
tomwalters@0: 
tomwalters@0:   extern short inputdata[MAX_DATA];                       
tomwalters@0:   extern struct Peak peak[MAX_PEAKS];               
tomwalters@0:   extern struct NextEvent nextevent[MAX_NEXTEVENTS];
tomwalters@0:   extern char progname[MAX_STRING_LENGTH];
tomwalters@0:   extern int framewidth_samples;           /* pwidth + nwidth * samplerate */
tomwalters@0: 
tomwalters@0:   int sample;                              /* left to right */
tomwalters@0:   int n = 0;                               /* peak counter: right to left */
tomwalters@0:   int zeroflag = OFF;                      /* used within nextevents */   
tomwalters@0:   int localsample, flat_top, truesample;   /* All 3 used in finding the 
tomwalters@0: 					      centre of a plateau */
tomwalters@0: 
tomwalters@0:   /* WARNING: goes from right to left ... */
tomwalters@0: 
tomwalters@0: /*---------------------------*/
tomwalters@0: 
tomwalters@0:   /* do peak check for first sample: its framewidth - 1, because the channel 
tomwalters@0:    * starts at 0 */
tomwalters@0:   /* NB: removed 18th June 1993*/
tomwalters@0:  
tomwalters@0: /*  if (inputdata[framewidth_samples - 1] > inputdata[framewidth_samples -2]){
tomwalters@0: *    peak[++n].sample = framewidth_samples - 1;
tomwalters@0: *    peak[n].mag = inputdata[framewidth_samples-1];}
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*---------------------------*/
tomwalters@0: 
tomwalters@0:   /* for the rest of the samples */
tomwalters@0:   for (sample=framewidth_samples - 2; sample >= 1; sample--) {
tomwalters@0: 
tomwalters@0:     /* simple peak ?  */
tomwalters@0:     if ((inputdata[sample] > inputdata[sample - 1]) && 
tomwalters@0: 	(inputdata[sample] > inputdata[sample + 1])){
tomwalters@0:       peak[++n].sample = sample;
tomwalters@0:       peak[n].mag = inputdata[sample];
tomwalters@0:       continue;}
tomwalters@0: 
tomwalters@0:     /* complicated peak? ie, a plateau to it  */
tomwalters@0:     /* but plateaus on the sides of big peaks don't count ... */
tomwalters@0:     if ((inputdata[sample] >= inputdata[sample - 1]) && 
tomwalters@0: 	(inputdata[sample] > inputdata[sample + 1])){
tomwalters@0:       localsample = sample;      /* counter to leftedge of the plateau */
tomwalters@0:       flat_top = 1;              /* width, in samples, of the top */
tomwalters@0:       truesample = sample;
tomwalters@0:       while (inputdata[--localsample] == inputdata[sample]) 
tomwalters@0: 	flat_top ++;            /* find out the width of the plateau */
tomwalters@0: 
tomwalters@0:       /* is this plateau a mere bump on the edge of a real peak? 
tomwalters@0:        * If so, its not a peak.
tomwalters@0:        * localsample is just off the leftedge of the plateau */
tomwalters@0:       if (inputdata[localsample] > inputdata[localsample+1])
tomwalters@0: 	continue;
tomwalters@0: 
tomwalters@0:       truesample = sample - (int) (flat_top / 2);  
tomwalters@0:             /* I think this ensures it is a leftward as possible */
tomwalters@0:       peak[++n].sample = truesample;
tomwalters@0:       peak[n].mag = inputdata[truesample];
tomwalters@0:       continue;}
tomwalters@0:     
tomwalters@0:     /*---------------------------*/
tomwalters@0: 
tomwalters@0:     /* nextevents ...  */
tomwalters@0: 
tomwalters@0:     /* This is an errorcheck ...
tomwalters@0:      *  If we've reached a LEFTZERO, but the current nextevent is a 
tomwalters@0:      *  RIGHTZERO, then something has gone very wrong ...  
tomwalters@0:      *  we should have hit a peak on the way, so incrementing the counter */
tomwalters@0:     if ((inputdata[sample] == 0) && nextevent[n].type == RIGHTEDGE) {
tomwalters@0:       fprintf(stderr, " %s : something's gone wrong: didn't find a peak between two zeros.\n", progname);
tomwalters@0:       exit(-1); }
tomwalters@0: 
tomwalters@0:     if ((inputdata[sample] == 0) && (zeroflag == OFF)) {
tomwalters@0:       nextevent[n].type = LEFTZERO;
tomwalters@0:       nextevent[n].sample = sample;
tomwalters@0:       nextevent[n].mag = inputdata[sample]; 
tomwalters@0:       zeroflag = ON;
tomwalters@0:       continue;}
tomwalters@0: 
tomwalters@0:     if ((inputdata[sample] != 0) && (zeroflag == ON)) {
tomwalters@0:       nextevent[n].type = RIGHTZERO;
tomwalters@0:       nextevent[n].sample = sample+1;    /* +1 because you can only notice 
tomwalters@0: 					  * such an event when you are off it,
tomwalters@0: 					  * and because we are moving 
tomwalters@0: 					  * in reverse  */
tomwalters@0:       nextevent[n].mag = inputdata[sample+1]; 
tomwalters@0:       zeroflag = OFF;
tomwalters@0:       continue;}
tomwalters@0:     
tomwalters@0:     /* This next bit removed at some stage */
tomwalters@0:     /*	
tomwalters@0:      *   if ((inputdata[sample] <= inputdata[sample +1]) && 
tomwalters@0:      *       (inputdata[sample] < inputdata[sample -1]) && 
tomwalters@0:      *       (nextevent[n].type == NULL)) { 
tomwalters@0:      */
tomwalters@0: 
tomwalters@0:     if ((inputdata[sample] <= inputdata[sample +1]) && 
tomwalters@0: 	(inputdata[sample] < inputdata[sample -1])) {
tomwalters@0:       nextevent[n].type = MINIMUM;
tomwalters@0:       nextevent[n].sample = sample;
tomwalters@0:       nextevent[n].mag = inputdata[sample];
tomwalters@0:       continue;}
tomwalters@0:  
tomwalters@0:   } /* sample */
tomwalters@0: 
tomwalters@0: /*-----------------------------*/
tomwalters@0:       
tomwalters@0:   /* repeat for last sample*/
tomwalters@0: 
tomwalters@0:   /*  peak ?  */
tomwalters@0:   /*NB: Removed 18th June 1993 */
tomwalters@0:   
tomwalters@0:   /*  if (inputdata[0] > inputdata[1] ){
tomwalters@0:    *    peak[++n].sample = 0;
tomwalters@0:    *    peak[n].mag = inputdata[0];}
tomwalters@0:    */
tomwalters@0:  
tomwalters@0: 
tomwalters@0:   /* nextevents ... */
tomwalters@0: 
tomwalters@0:   if (inputdata[0] == 0 && zeroflag == OFF) {
tomwalters@0:     nextevent[n].type = LEFTZERO;
tomwalters@0:     nextevent[n].sample = 0;
tomwalters@0:     nextevent[n].mag = inputdata[0]; 
tomwalters@0:     zeroflag = ON;}
tomwalters@0:   
tomwalters@0:   
tomwalters@0:   if (inputdata[0] != 0 && zeroflag == ON) {
tomwalters@0:     nextevent[n].type = RIGHTZERO;
tomwalters@0:     nextevent[n].sample = 0;  /* -1 because you only notice such an event 
tomwalters@0: 			       * when you've gone past*/
tomwalters@0:     nextevent[n].mag = inputdata[0]; 
tomwalters@0:     zeroflag = OFF;}
tomwalters@0:   
tomwalters@0:   if (inputdata[0] < inputdata[1]){
tomwalters@0:     nextevent[n].type = MINIMUM;
tomwalters@0:     nextevent[n].sample = 0;
tomwalters@0:     nextevent[n].mag = inputdata[0]; }
tomwalters@0:   
tomwalters@0:   if (nextevent[n].type == UNSET) {
tomwalters@0:     nextevent[n].type = END_OF_WIDTH;
tomwalters@0:     nextevent[n].sample = 0;
tomwalters@0:     nextevent[n].mag = inputdata[0]; }
tomwalters@0: 
tomwalters@0:   return n;  /* total number of peaks found, starting at 1 */
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* The End */
tomwalters@0: /*--------------------------------------------------------------------------*/
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: