view 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 */
/*--------------------------------------------------------------------------*/