diff saitools/saigraph.c @ 0:5242703e91d3 tip

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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/saitools/saigraph.c	Fri May 20 15:19:45 2011 +0100
@@ -0,0 +1,465 @@
+/*
+    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.
+*/
+
+/*------------------------------------------------------------------------*/
+
+/*  saigraph.c
+*  --------------
+*  
+* Works out the (first-order) temporal intervals between peaks in a .sai file,
+* and saves them in a SINGLE-FRAME .sai
+* In essence, a histogram.
+*
+* Uses ALL peaks: if you just want the locallly maxima peaks, then use 
+* trbigpeak first.
+*
+* Output is another .sai file: but only a single frame. The header is
+* modified accordingly.
+* 
+*
+*
+* M. Akeroyd.   May 1993.   Version 2.0
+*
+*               Autumn 1993: added grayscaleflag
+*               Winter 1994: allowed weighting, and other general revisions.
+*/
+
+
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include "tip.h"
+
+
+/* Function declarations */
+/*------------------------------------------------------------*/
+
+void parsecommandline (int argc, char *argv[], char inputfn[], char outputfn[]);
+
+int readheader(FILE *inputfp);
+void checkheader();
+void writeheader(FILE *outputfp);
+
+void changeheader(int no_frames, int new_pwidth, int new_nwidth, int framewidth_samples_output);
+
+int findpeaksreverse (); 
+void find_intervals(int total_peaks, int channel, int max_interval);
+
+void writedata_output (FILE *outputfp, int samples_to_write);
+void writedata_output_fg (FILE *figurefp, FILE *groundfp, int samples_to_write);
+FILE *open_file (char filefn[], FILE *dir_default, int streamtype);
+
+void copytooutput_fg_single (int total_peaks, int frame, int channel, int trigger_peak);
+
+void write_matrix_hori_time(int framewidth, FILE *outputfigurefp);
+void write_matrix_hori_freq(int framewidth, FILE *outputfigurefp);
+
+
+
+
+/* Data arrays: global */
+/*-----------------------------------------------------------------*/
+
+
+short inputdata[MAX_DATA];        /* input .sai: a single channel of a frame */
+short outputdata[MAX_DATA];       /* output .sai: a single channel */
+short outputfiguredata[MAX_DATA]; /* output figure.sai: a single channel */
+short outputgrounddata[MAX_DATA]; /* output ground.sai: a single channel */
+short clip[MAX_CHANNELS];
+
+short interval[MAX_CHANNELS][MAX_DATA];      /* histogram */ 
+short summation[MAX_CHANNELS];
+
+
+
+/* other variables */
+/*----------------------------------------------------------------*/
+
+struct Peak peak[MAX_PEAKS];                    /* used in findpeaksreverse():
+						 * single channel */
+struct NextEvent nextevent[MAX_NEXTEVENTS];     /* used in findpeaksreverse():
+						 * single channel */
+
+
+
+/* variables read from header */
+/*---------------------------------------------------------------*/
+
+char header[MAX_LINES_HEADER][MAX_LINE_LENGTH];
+int header_lines;
+
+int no_frames;
+int frameheight;             /* number of channels */
+int framewidth_samples;      /* pwidth + nwidth * samplerate */
+int frameshift_samples;      /* frstep_aid * samplerate */
+int pwidth;                  /* in msecs */
+int nwidth;                  /* in msecs: NEGATIVE */
+int width_win;               /* pixels */
+int height_win;              /* pixels */
+long samplerate;             /* samples per sec */
+int mincf;                   /* Hz */
+int maxcf;                   /* Hz */
+
+
+
+
+/* misc */
+/*-----------------------------------------------------------*/
+
+char progname[MAX_STRING_LENGTH];
+char outputfigurefn[MAX_STRING_LENGTH];
+
+
+int verboseflag = OFF;               /* -v */
+int removeflag = OFF;                /* -r */
+double removevalue = 0.1;            
+int widthflag = OFF;                 /* -width : if OFF, copy pwdith/nwidth 
+				      *      off Input .sai */
+int total_width;                     /* -width : total width of new .sai */
+int asciiflag = OFF;                 /* -a */
+int asciireverseflag = OFF;          /* -ar */
+int summationflag = OFF;             /* -s */
+int greyscaleflag = OFF;             /* -g */
+int weightflag = OFF;             /* -weight */
+int oppositearchflag = OFF;          /* -oparch: if ON, then the header-reading
+				      * code will load one (1) extra byte. This
+				      * is so that Sun .sai/.nap can be read
+				      * on a DEC, and equally the reverse.
+				      * It is NOT required if reading Sun on
+				      * Sun, or DEC on DEC. */
+
+
+
+/*...............           Main           ................................*/
+/*.........................................................................*/
+/*.........................................................................*/
+
+
+
+void main (int argc, char *argv[])
+{
+  int sample;
+  int n;
+
+  int frame, channel;
+  int trigger_peak = 0;
+  int no_peak = 0, total_peaks =0, total_peaks_2;
+  int framewidth_samples_input,
+      framewidth_samples_output;
+  int new_pwidth, 
+      new_nwidth;
+  int header_bytes = 0;
+ 
+  char inputfn[MAX_STRING_LENGTH], 
+       outputbasefn[MAX_STRING_LENGTH];
+
+  FILE *inputfp, *outputfigurefp;
+
+  int clipflag = OFF;
+
+/*-----------------------------------------*/
+
+  strcpy(progname, argv[0]);
+  strcpy(inputfn, "");
+  strcpy(outputbasefn, "");
+  strcpy(outputfigurefn, "");
+
+  /* parse command line */
+  parsecommandline(argc, argv, inputfn, outputbasefn);
+
+  if (strcmp(outputbasefn, "") == 0) 
+    strcpy(outputfigurefn, "");
+  else 
+    strcpy(outputfigurefn, outputbasefn);
+
+
+  /* open files */
+  /* default directions are:
+   * input = stdin 
+   * figure = stdout
+   */
+
+  inputfp = open_file(inputfn, stdin, READ);
+  outputfigurefp = open_file(outputfigurefn, stdout, WRITE);
+
+  /* read Header */
+  header_bytes = readheader(inputfp);
+  checkheader();
+
+
+  /* reset the variables */
+  framewidth_samples_input = framewidth_samples;
+  if (widthflag == OFF) {
+    new_pwidth = + pwidth;
+    new_nwidth = -nwidth;
+    framewidth_samples_output = framewidth_samples_input;}
+  else {
+    new_pwidth = + NEW_PWIDTH;
+    new_nwidth = (new_pwidth + total_width);
+    framewidth_samples_output = (int) ((int) (abs(new_pwidth) + abs (new_nwidth)) * samplerate ) / 1000;}
+
+  if (framewidth_samples_output >= MAX_DATA) {
+    fprintf(stderr, "%s: new frame is  too wide at %i: only allowed %i samples\n", progname, framewidth_samples_output, MAX_DATA);
+    exit(-1); }
+
+
+  /* change header */
+  changeheader(1, new_pwidth, new_nwidth, framewidth_samples_output);
+
+  /* write header */
+  if (asciiflag == OFF) 
+    writeheader(outputfigurefp);
+  else {
+    /* begining of ascii header. The rest depends on the matrix format ... */
+    fprintf(outputfigurefp, "# saigraph output.  %s\n", inputfn);
+    fprintf(outputfigurefp, "# \n");
+    fprintf(outputfigurefp, "# samplerate=%i  pwidth=%d  nwidth=%d  frames=%i ", samplerate, pwidth, nwidth, no_frames);
+    fprintf(outputfigurefp, " mincf=%i  maxcf=%i  channels=%i \n", mincf, maxcf,  frameheight);
+    fprintf(outputfigurefp, "# \n");
+  }
+
+/*------------------------------------------*/
+
+  /* Main loop */
+
+  if  (verboseflag == ON) { 
+    fprintf(stderr, " frames "); 
+    fflush(stderr);}
+
+  
+  /* clear arrays */
+  
+  for (channel=0; channel<MAX_CHANNELS; channel++) {
+    for(sample=0; sample<MAX_DATA; sample++)
+      interval[channel][sample] = 0;
+    clip[channel]=OFF;
+  }
+
+
+  /* reset this, so everything still works */
+  framewidth_samples = framewidth_samples_input;
+
+
+  for (frame=1; frame<=no_frames; frame++) {
+
+    if  (verboseflag == ON) { 
+      fprintf(stderr, " %i ", frame); 
+      fflush(stderr);}
+
+    /*--------------------------------------------*/
+
+    /* load the whole frame's worth of data first */
+
+    for (channel=1; channel <=frameheight; channel ++) {
+
+      /* clear arrays */
+      for(sample=0; sample<framewidth_samples; sample++) 
+	inputdata[sample] = outputdata[sample] = 0;
+      for(n=0; n < MAX_PEAKS; n++){
+	peak[n].tent = UNSET;
+	peak[n].sample=0;}
+
+      /* Nextevents are unnecessary */
+
+      /* load a single channel's worth of data */
+      fread (inputdata, 2, (size_t) framewidth_samples, inputfp);
+
+      /* This next is a simple input check: if any numbers < 0, then say so */
+      for(sample=0; sample<framewidth_samples; sample++)
+	if (inputdata[sample] < 0 )  {
+	  fprintf(stderr, "%s: something's gone wrong: the data is negative.\n", progname);
+	  exit(-1); }
+
+      /* find peaks. WARNING: in reverse */
+      total_peaks = findpeaksreverse();
+
+      /* bin any small peaks ... */
+      if ((removeflag == ON) && (total_peaks !=0 ) ) {
+        total_peaks_2 = removepeaks(total_peaks);
+        for(n=0; n < MAX_PEAKS; n++) {
+          peak[n].tent = UNSET;
+          peak[n].sample = 0;}
+        total_peaks = findpeaksreverse();}
+
+      /* get intervals */
+      find_intervals(total_peaks, channel, framewidth_samples_output);
+
+    } /* channel */
+
+  } /* frame */
+
+
+/*------------------------------------------*/
+/*------------------------------------------*/
+
+
+  fflush(stdout);
+
+
+  /* If required, sum the values */
+  if (summationflag==ON) {
+
+    for (channel = 1; channel <= frameheight; channel ++)
+	summation[channel] = 0;
+
+    for (channel = 1; channel <= frameheight; channel ++) {
+      clipflag = OFF;
+      for (sample = 0; sample < framewidth_samples_output; sample++)
+	summation[channel] += interval[channel][sample];
+
+      /* clip check */
+      if ((summation[channel] > 32765) || 
+	  (summation[channel] < 0)){
+	    if (clipflag == OFF) {
+	      fprintf(stderr, " \nclipping: channel %i, summations.\n ", channel);
+	      clipflag = ON;}
+	    summation[channel] = 32765;
+	  }
+    }
+  }
+
+/*-------------------------------------------*/
+
+  /* write output */
+
+  if (asciiflag == OFF) {
+    for (channel = 1; channel <= frameheight; channel ++){
+      for (sample = 0; sample < framewidth_samples_output; sample++)
+	outputfiguredata[sample] = interval[channel][sample];
+      writedata_output(outputfigurefp, framewidth_samples_output);
+    }
+  }  
+  else {
+    if (asciireverseflag == OFF) 
+      write_matrix_hori_freq(framewidth_samples_output, outputfigurefp);
+    else
+      write_matrix_hori_time(framewidth_samples_output, outputfigurefp);
+  }    
+  
+  /*-------------------------------------------*/
+
+  /* Tidy up and exit */
+    
+  if (verboseflag == ON) 
+    fprintf(stderr, "\n");
+  fclose(inputfp); 
+  fclose(outputfigurefp); 
+ 
+  if ( ferror(inputfp) != 0) {
+    fprintf(stderr, " %s : error closing input file.\n", progname);
+    exit(-1);}
+
+  if ( ferror(outputfigurefp) != 0) {
+    fprintf(stderr, " %s : error closing figure file.\n", progname);
+    exit(-1);}
+
+  exit(0);
+
+} /* Main */
+
+
+
+
+
+
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
+
+
+
+
+void parsecommandline(int argc, char *argv[], char inputfn[], char outputbasefn[])
+{
+  int x=1, helpflag = OFF;
+
+  while (x < argc){
+    if      (!strcmp(argv[x], "-o")) {strcpy(outputbasefn, argv[x+1]); x+=2;}
+    else if (!strcmp(argv[x], "-output")) {strcpy(outputbasefn, argv[x+1]); x+=2;}
+    else if (!strcmp(argv[x], "-i")) {strcpy(inputfn, argv[x+1]); x+=2;}
+    else if (!strcmp(argv[x], "-input")) {strcpy(inputfn, argv[x+1]); x+=2;}
+    else if (!strcmp(argv[x], "-help")) {helpflag = ON; x+=1;}
+    else if (!strcmp(argv[x], "-h")) {helpflag = ON; x+=1;}
+    else if (!strcmp(argv[x], "-verbose")) {verboseflag = ON; x+=1;}
+    else if (!strcmp(argv[x], "-v")) {verboseflag = ON; x+=1;}
+    else if (!strcmp(argv[x], "-weight")) {weightflag = NORMAL_WEIGHTING; x+=1;}
+    else if (!strcmp(argv[x], "-weightlog")) {weightflag = LOG_WEIGHTING; x+=1;}
+    else if (!strcmp(argv[x], "-width")) {widthflag = ON; total_width=atoi(argv[x+1]); x+=2;}
+    else if (!strcmp(argv[x], "-r")) { removeflag = ON; removevalue = atof(argv[x+1]); x+=2;}
+    else if (!strcmp(argv[x], "-a")) { asciiflag = ON;x+=1;}
+    else if (!strcmp(argv[x], "-ar")) { asciiflag = ON; asciireverseflag=ON; x+=1;}
+    else if (!strcmp(argv[x], "-s")) { summationflag=ON; x+=1;}
+    else if (!strcmp(argv[x], "-g")) { greyscaleflag=ON; x+=1;}
+    else if (!strcmp(argv[x], "-grey")) { greyscaleflag=ON; x+=1;}
+    else if (!strcmp(argv[x], "-gray")) { greyscaleflag=ON; x+=1;}
+    else if (!strcmp(argv[x], "-oparch")) { oppositearchflag=ON; x+=1;}
+    else {fprintf(stderr, "%s: unknown option %s\n", progname, argv[x]);
+	  exit(-1);}
+  }
+  
+  if (helpflag == ON)
+    {
+      fprintf(stderr, "\n");
+      fprintf(stderr, "                         %s\n", progname);
+      fprintf(stderr, " -----------------------------------------------------------------------------\n");
+      fprintf(stderr, " Works out a histogram of peak intervals, from a .sai\n");
+      fprintf(stderr, " Default output is a .sai file: if -a or -ar specified,\n");
+
+      fprintf(stderr, " then output is a ASCII matrix.\n");
+      fprintf(stderr, "\n");
+      fprintf(stderr, " The -oparch option is REQUIRED if reading a Sun .sai files \n");
+      fprintf(stderr, " on a DEC cpu or vice-versa.\n");
+      fprintf(stderr, " Add .sai to both input and output filenames.\n");
+
+      fprintf(stderr, " -----------------------------------------------------------------------------\n");
+      fprintf(stderr, "\n");
+      fprintf(stderr, " -i <.sai file>      input file             default = stdin\n");
+      fprintf(stderr, " -o <.sai file>      output file            default = stdout\n");
+      fprintf(stderr, " -oparch             assume input was generated on an opposing architecture\n");
+      fprintf(stderr, "\n");
+      fprintf(stderr, " -a                  ASCII output: rows = channels\n");
+      fprintf(stderr, " -ar                 ASCII output: rows = time intervals\n");
+      fprintf(stderr, " -s                  include summations over time intervals (ASCII options only)\n");
+      fprintf(stderr, "\n");
+      fprintf(stderr, " -width <int>        width of graph.sai (msecs).\n");
+      fprintf(stderr, " -r <fraction>       remove any peaks F(n) < F(n+1) & F(n-1) * fract \n");
+      fprintf(stderr, "                                            default = %.3f\n", removevalue);
+      fprintf(stderr, " -weight             weight histogram entries by mean height of peaks\n");
+      fprintf(stderr, " -weightlog          as per -weight, but use 10log(mean height)\n");
+      fprintf(stderr, "\n");
+      fprintf(stderr, " -v                  verbose output (to stderr)\n");
+      fprintf(stderr, " -g                  change .sai 'view' format to 'grayscale'.\n");
+      fprintf(stderr, "\n\n");
+      exit(-1);}
+
+}  
+
+
+
+
+/* The End */
+/*---------------------------------------------------------------------------*/