comparison saitools/findpeaks.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

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