tomwalters@0: /*
tomwalters@0:     op.c        Ordered sequence of operations on input stream.
tomwalters@0:    ------
tomwalters@0: 
tomwalters@0:     Operations on each input item are done in the order they appear on the
tomwalters@0:     command line. Operations can be repeated as many times as needed on the
tomwalters@0:     command line.
tomwalters@0: 
tomwalters@0:     Input items may be binary numbers (short or float) or ascii lines
tomwalters@0:     (terminated by CR) according to the type option.
tomwalters@0: 
tomwalters@0:     The list of operations is as follows. Some operations take valued
tomwalters@0:     arguments (which are real numbers), others are flags, as indicated:
tomwalters@0: 
tomwalters@0:     operation           result for each input item x
tomwalters@0:     -----------------   --------------------------------------------
tomwalters@0:     add=<value>         x+<value>
tomwalters@0:     negate=on           -x
tomwalters@0:     multiply=<value>    x*<value>
tomwalters@0:     divide=<value>      x/<value>
tomwalters@0:     remainder=<value>   real remainder of x/<value>
tomwalters@0:     inverse=on          1/x
tomwalters@0:     round=on            x rounded to nearest integer
tomwalters@0:     absolute=on         |x|
tomwalters@0:     power=<value>       x^<value>
tomwalters@0:     exponent=<value>    <value>^x
tomwalters@0:     e-exponent=on       e^x
tomwalters@0:     log=<value>         log(x) to base <value>
tomwalters@0:     e-log=on            log(x) to base e.
tomwalters@0:     sin=on              sin(x)
tomwalters@0:     cos=on              cos(x)
tomwalters@0:     tan=on              tan(x)
tomwalters@0:     threshold=<value>   if ( x  <value> ) 0 else if ( x = <value> ) 1
tomwalters@0:     diff1=on            first difference
tomwalters@0:     diff2=on            second difference
tomwalters@0:     diff3=on            third difference
tomwalters@0:     cusum=on            cumulative sum
tomwalters@0:     cumean=on           cumulative (recursive) mean an
tomwalters@0:     cumin=on            cumulative min
tomwalters@0:     cumax=on            cumulative max
tomwalters@0: 
tomwalters@0:     The `type' option sets the data type for all input and output items.
tomwalters@0:     Recognised data types are: short, float, ascii.
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     Ascii data is read in line-by-line, and comment lines are allowed.
tomwalters@0:     The following operation apply only to ascii input items, enabling
tomwalters@0:     comment lines to be stripped or echoed directly to the output.
tomwalters@0: 
tomwalters@0:     strip=<value>       ignore lines beginning with <value> string
tomwalters@0:     echo=<value>        echo lines beginning with <value> string
tomwalters@0:     number=on           number output lines
tomwalters@0: 			If number=on then stripped lines are
tomwalters@0: 			not numbered; they simply dissappear.
tomwalters@0: 
tomwalters@0:     If there are no operations, (other than comment stripping/echoing or line
tomwalters@0:     numbering), then ascii lines are echoed.
tomwalters@0: 
tomwalters@0: 
tomwalters@0: Examples:
tomwalters@0: 
tomwalters@0: 1. To specify the following ordered sequence of operations on each binary
tomwalters@0: float in the input stream: invert (take reciprocal), scale up by 100,
tomwalters@0: square (raise to power 2), invert (take reciprocal), square root (raise to
tomwalters@0: power 0.5).
tomwalters@0: 
tomwalters@0: op type=float inverse=on multiply=100 power=2 inverse=on power=.5  file
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 2. To use as an interactive "desk calculator", eg to print out the result
tomwalters@0: of log2(17/16):
tomwalters@0: 
tomwalters@0: echo 17 | op type=ascii -div16 -log2
tomwalters@0: 
tomwalters@0: 3. To interactively operate on a list, eg add 1 to each of the list 1,2,3,
tomwalters@0: (note how to put newlines into echo'd lists):
tomwalters@0: 
tomwalters@0: echo "1\
tomwalters@0: 2\
tomwalters@0: 3" | op type=ascii -add1
tomwalters@0: 
tomwalters@0: 4. To read ascii data, echoing any comment lines beginning with the string
tomwalters@0: "**", and for each number x in infile, o/p the number 149-x, do:
tomwalters@0: 
tomwalters@0: op type=ascii echo="**" neg=on -add149  infile > outfile
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "options.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: 
tomwalters@0: char applic[] = "Ordered sequence of operations on each item in input stream.\n     (Input and output in binary shorts)." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr,
tomwalters@0: 		*addstr ,
tomwalters@0: 		*negstr ,
tomwalters@0: 		*mulstr ,
tomwalters@0: 		*divstr ,
tomwalters@0: 		*remstr ,
tomwalters@0: 		*invstr ,
tomwalters@0: 		*rndstr ,
tomwalters@0: 		*absstr ,
tomwalters@0: 		*powstr ,
tomwalters@0: 		*expstr ,
tomwalters@0: 		*eexpstr,
tomwalters@0: 		*logstr ,
tomwalters@0: 		*elogstr,
tomwalters@0: 		*sinstr ,
tomwalters@0: 		*cosstr ,
tomwalters@0: 		*tanstr ,
tomwalters@0: 		*threshstr,
tomwalters@0: 		*diff1str  ,
tomwalters@0: 		*diff2str  ,
tomwalters@0: 		*diff3str  ,
tomwalters@0: 		*cusumstr   ,
tomwalters@0: 		*cumeanstr  ,
tomwalters@0: 		*cuminstr   ,
tomwalters@0: 		*cumaxstr   ,
tomwalters@0: 		*typestr ,
tomwalters@0: 		*stripstr ,
tomwalters@0: 		*echostr  ,
tomwalters@0: 		*numstr    ;
tomwalters@0: 
tomwalters@0: static Options option[] = {
tomwalters@0:     {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                                  , DEBUG   },
tomwalters@0:     {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                                      , DEBUG   },
tomwalters@0:     {   "add"       ,   "off"       ,  &addstr      ,   "Add                    x = x+value  "                  , VAL     },
tomwalters@0:     {   "negate"    ,   "off"       ,  &negstr      ,   "Negate                 x = -x     "                    , SETFLAG },
tomwalters@0:     {   "multiply"  ,   "off"       ,  &mulstr      ,   "Multiply               x = x*value  "                  , VAL     },
tomwalters@0:     {   "divide"    ,   "off"       ,  &divstr      ,   "Divide                 x = x/value  "                  , VAL     },
tomwalters@0:     {   "remainder" ,   "off"       ,  &remstr      ,   "Remainder              x = real remainder of x/value"  , VAL     },
tomwalters@0:     {   "inverse"   ,   "off"       ,  &invstr      ,   "Inverse                x = 1/x    "                    , SETFLAG },
tomwalters@0:     {   "round"     ,   "off"       ,  &rndstr      ,   "Round                  x = x, to nearest integer"      , SETFLAG },
tomwalters@0:     {   "absolute"  ,   "off"       ,  &absstr      ,   "Absolute               x = |x|    "                    , SETFLAG },
tomwalters@0:     {   "power"     ,   "off"       ,  &powstr      ,   "Power                  x = x^value  "                  , VAL     },
tomwalters@0:     {   "exponent"  ,   "off"       ,  &expstr      ,   "Exponent               x = value^x  "                  , VAL     },
tomwalters@0:     {   "e-exponent",   "off"       ,  &eexpstr     ,   "e-exponent             x = e^x    "                    , SETFLAG },
tomwalters@0:     {   "log"       ,   "off"       ,  &logstr      ,   "Log                    x = log(x), base `value'."      , VAL     },
tomwalters@0:     {   "e-log"     ,   "off"       ,  &elogstr     ,   "elog                   x = log(x), base e.    "        , SETFLAG },
tomwalters@0:     {   "sin"       ,   "off"       ,  &sinstr      ,   "sin                    x = sin(x) "                    , SETFLAG },
tomwalters@0:     {   "cos"       ,   "off"       ,  &cosstr      ,   "cos                    x = cos(x) "                    , SETFLAG },
tomwalters@0:     {   "tan"       ,   "off"       ,  &tanstr      ,   "tan                    x = tan(x) "                    , SETFLAG },
tomwalters@0:     {   "threshold" ,   "off"       ,  &threshstr   ,   "Threshold              x = 0 (x<value) or 1 (x>=value)", VAL     },
tomwalters@0:     {   "diff1"     ,   "off"       ,  &diff1str    ,   "First difference "                                     , SETFLAG },
tomwalters@0:     {   "diff2"     ,   "off"       ,  &diff1str    ,   "Second difference "                                    , SETFLAG },
tomwalters@0:     {   "diff3"     ,   "off"       ,  &diff1str    ,   "Third difference "                                     , SETFLAG },
tomwalters@0:     {   "cusum"     ,   "off"       ,  &cusumstr    ,   "Cumulative sum         x = current sum"                , SETFLAG },
tomwalters@0:     {   "cumean"    ,   "off"       ,  &cumeanstr   ,   "Cumulative mean        x = current (recursive) mean"   , SETFLAG },
tomwalters@0:     {   "cumin"     ,   "off"       ,  &cuminstr    ,   "Cumulative min         x = current min  "              , SETFLAG },
tomwalters@0:     {   "cumax"     ,   "off"       ,  &cumaxstr    ,   "Cumulative max         x = current max  "              , SETFLAG },
tomwalters@0:     {   "type"      ,   "short"     ,  &typestr     ,   "Data type (short, float, ascii)"                       , VAL     },
tomwalters@0:     {   "strip"     ,   "off"       ,  &stripstr    ,   "Ignore lines beginning with value string"              , VAL     },
tomwalters@0:     {   "echo"      ,   "off"       ,  &echostr     ,   "Echo lines beginning with value string"                , VAL     },
tomwalters@0:     {   "number"    ,   "off"       ,  &numstr      ,   "Number output lines"                                   , SETFLAG },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #define SIZE        64  /* Max number of program command-line arguments */
tomwalters@0: 
tomwalters@0: #define Add          2  /* Functions */
tomwalters@0: #define Negate       3
tomwalters@0: #define Multiply     4
tomwalters@0: #define Divide       5
tomwalters@0: #define Remainder    6
tomwalters@0: #define Inverse      7
tomwalters@0: #define Round        8
tomwalters@0: #define Absolute     9
tomwalters@0: #define Power       10
tomwalters@0: #define Exp         11
tomwalters@0: #define Eexp        12
tomwalters@0: #define Log         13
tomwalters@0: #define Elog        14
tomwalters@0: #define Sin         15
tomwalters@0: #define Cos         16
tomwalters@0: #define Tan         17
tomwalters@0: #define Threshold   18
tomwalters@0: #define Diff1       19
tomwalters@0: #define Diff2       20
tomwalters@0: #define Diff3       21
tomwalters@0: #define Sum         22
tomwalters@0: #define Mean        23
tomwalters@0: #define Min         24
tomwalters@0: #define Max         25
tomwalters@0: #define Ascii       26
tomwalters@0: #define Strip       27
tomwalters@0: #define Echo        28
tomwalters@0: #define Number      29
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* data types */
tomwalters@0: 
tomwalters@0: #define SHORT        0
tomwalters@0: #define FLOAT        1
tomwalters@0: #define ASCII        2
tomwalters@0: 
tomwalters@0: int     Type ;
tomwalters@0: 
tomwalters@0: int     STRIP=0;
tomwalters@0: char    Stripstr[]="*";         /* Default strip comment string */
tomwalters@0: 
tomwalters@0: int     ECHO=0;
tomwalters@0: char    Echostr[]="**";         /* Default echo comment string */
tomwalters@0: 
tomwalters@0: int     NUMBER=0;               /* Flag for numbering ascii lines */
tomwalters@0: int     n=0;
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:  Rounding for both +ve and -ve numbers. Real numbers are truncated to ints in
tomwalters@0:  the direction of zero (EG both 0.9 and -0.9 are truncated to 0). Using round,
tomwalters@0:  0.9 becomes 1 and -0.9 becomes -1. (Note also 0.5 becomes 1, -0.5 becomes -1.
tomwalters@0: */
tomwalters@0: #define round(x)    ((x >= 0) ? ((int)(x+0.5)) : ((int)(x-0.5)))
tomwalters@0: 
tomwalters@0: /* Remainder of real division for +ve and -ve numbers */
tomwalters@0: #define rem(x,y)    ((x >= 0) ? (fabs((x) - (y)*(int)((x)/(y)))) : (fabs((x) - (y)*(int)((x)/(y)-0.5))))
tomwalters@0: 
tomwalters@0: /* Threshold and convert to binary (0,1) value */
tomwalters@0: #define threshold(x,T)  ((x >= T) ? 1 : 0)
tomwalters@0: 
tomwalters@0: float diff1();
tomwalters@0: float diff2();
tomwalters@0: float diff3();
tomwalters@0: float mean();
tomwalters@0: float min();
tomwalters@0: float max();
tomwalters@0: float sum();
tomwalters@0: 
tomwalters@0: struct opstr {      /* array of operations: functions and (optional) arg */
tomwalters@0:     char  func;
tomwalters@0:     float arg;
tomwalters@0: } op[SIZE];
tomwalters@0: 
tomwalters@0: int k=0;            /* number of operations stored in array */
tomwalters@0: 
tomwalters@0: main(argc, argv)
tomwalters@0: int   argc ;
tomwalters@0: char *argv[] ;
tomwalters@0: {
tomwalters@0:     FILE *fp, *fopen();
tomwalters@0:     float xf, ops();
tomwalters@0:     short xs ;
tomwalters@0:     char  s[256], stripcomment[32], echocomment[32];
tomwalters@0:     char *prog, *val ;
tomwalters@0:     int   i, index, span ;
tomwalters@0: 
tomwalters@0:     strcpy(stripcomment, Stripstr);
tomwalters@0:     strcpy(echocomment,  Echostr );
tomwalters@0: 
tomwalters@0:     for (i=0 ; option[i].name != (char *)0 ; i++)
tomwalters@0: 	*(option[i].val) = option[i].dflt ;
tomwalters@0: 
tomwalters@0:     prog = argv[0] ;
tomwalters@0:     while ( --argc > 0 && ++argv ) {
tomwalters@0: 	if ( ( i = whichopt( option, *argv, &index, &span ) ) == UNKNOWN )
tomwalters@0: 	    break ;
tomwalters@0: 	if ( ( val = checksyntax( *argv, span, option[index].type ) ) == (char *)0 )
tomwalters@0: 	    break ;
tomwalters@0: 	if ( i == AMBIGUOUS ) {
tomwalters@0: 	    fprintf(stderr,"%s: ambiguous option [%s]\n", prog, *argv ) ;
tomwalters@0: 	    exit ( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 	switch ( index ) {
tomwalters@0: 	    case  0 : operate( option, index,  val ) ;                      break;
tomwalters@0: 	    case  1 : operate( option, index,  val ) ;                      break;
tomwalters@0: 	    case  2 : op[k].func = Add;         op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case  3 : op[k].func = Negate;                             k++; break;
tomwalters@0: 	    case  4 : op[k].func = Multiply;    op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case  5 : op[k].func = Divide;      op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case  6 : op[k].func = Remainder;   op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case  7 : op[k].func = Inverse;                            k++; break;
tomwalters@0: 	    case  8 : op[k].func = Round;                              k++; break;
tomwalters@0: 	    case  9 : op[k].func = Absolute;                           k++; break;
tomwalters@0: 	    case 10 : op[k].func = Power;       op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case 11 : op[k].func = Exp;         op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case 12 : op[k].func = Eexp;                               k++; break;
tomwalters@0: 	    case 13 : op[k].func = Log;         op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case 14 : op[k].func = Elog;                               k++; break;
tomwalters@0: 	    case 15 : op[k].func = Sin;                                k++; break;
tomwalters@0: 	    case 16 : op[k].func = Cos;                                k++; break;
tomwalters@0: 	    case 17 : op[k].func = Tan;                                k++; break;
tomwalters@0: 	    case 18 : op[k].func = Threshold;   op[k].arg = atof(val); k++; break;
tomwalters@0: 	    case 19 : op[k].func = Diff1;                              k++; break;
tomwalters@0: 	    case 20 : op[k].func = Diff2;                              k++; break;
tomwalters@0: 	    case 21 : op[k].func = Diff3;                              k++; break;
tomwalters@0: 	    case 22 : op[k].func = Sum;                                k++; break;
tomwalters@0: 	    case 23 : op[k].func = Mean;                               k++; break;
tomwalters@0: 	    case 24 : op[k].func = Min;                                k++; break;
tomwalters@0: 	    case 25 : op[k].func = Max;                                k++; break;
tomwalters@0: 	    case 26 : Type = checktype( val ) ;                             break;
tomwalters@0: 	    case 27 : STRIP++; strcpy(stripcomment, val);                   break;
tomwalters@0: 	    case 28 : ECHO++;  strcpy(echocomment,  val);                   break;
tomwalters@0: 	    case 29 : NUMBER++;                                             break;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     }
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts( helpstr, prog, applic, option ) ;
tomwalters@0: 
tomwalters@0:     if (argc <= 0) fp = stdin;
tomwalters@0:     else if ((fp = fopen(*argv, "r")) == NULL) {
tomwalters@0: 	    fprintf(stderr,"can't open %s\n", *argv);
tomwalters@0: 	    exit(1);
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     switch ( Type ) {
tomwalters@0: 
tomwalters@0: 	case SHORT :
tomwalters@0: 	    while ( fread( &xs, sizeof(short), 1, fp ) ) {
tomwalters@0: 		xs = (short)ops( (float)xs ) ;
tomwalters@0: 		fwrite( &xs, sizeof(short), 1, stdout ) ;
tomwalters@0: 	    }
tomwalters@0: 	    break ;
tomwalters@0: 
tomwalters@0: 	case FLOAT :
tomwalters@0: 	    while ( fread( &xf, sizeof(float), 1, fp ) ) {
tomwalters@0: 		xf = ops( xf ) ;
tomwalters@0: 		fwrite( &xf, sizeof(float), 1, stdout ) ;
tomwalters@0: 	    }
tomwalters@0: 	    break ;
tomwalters@0: 
tomwalters@0: 	case ASCII :
tomwalters@0: 	    while (fgets(s, 256, fp)) {
tomwalters@0: 		if (STRIP && match(s,stripcomment)==0)
tomwalters@0: 		    ;                       /* strip comment (ie do nothing) */
tomwalters@0: 		else {
tomwalters@0: 		    if (NUMBER) printf("%d: ", ++n);
tomwalters@0: 		    if (ECHO && match(s,echocomment)==0)
tomwalters@0: 			fputs(s,stdout);        /* echo comment */
tomwalters@0: 		    else {
tomwalters@0: 			if ( k>0 ) {
tomwalters@0: 			    xf = ops( atof(s) ) ;
tomwalters@0: 			    printf("%.3f\n", xf ) ;
tomwalters@0: 			}
tomwalters@0: 			else fputs(s,stdout);   /* echo line (case of no ops) */
tomwalters@0: 		    }
tomwalters@0: 		}
tomwalters@0: 	    }
tomwalters@0: 	    break ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     fclose(fp);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: checktype( s )
tomwalters@0: char *s ;
tomwalters@0: {
tomwalters@0:     if ( iststr( s, "short" ) )  return SHORT ;
tomwalters@0:     if ( iststr( s, "float" ) )  return FLOAT ;
tomwalters@0:     if ( iststr( s, "ascii" ) )  return ASCII ;
tomwalters@0:     fprintf( stderr,"unknown datatype [%s]\n", s ) ;
tomwalters@0:     exit( 1 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float ops(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     int i;
tomwalters@0: 
tomwalters@0:     for (i=0 ; i<k ; i++)
tomwalters@0: 	switch (op[i].func) {
tomwalters@0: 	    case Add:           x = x+op[i].arg;            break;
tomwalters@0: 	    case Negate:        x = (-x);                   break;
tomwalters@0: 	    case Multiply:      x = x*op[i].arg;            break;
tomwalters@0: 	    case Divide:        x = x/op[i].arg;            break;
tomwalters@0: 	    case Remainder:     x = rem(x,op[i].arg);       break;
tomwalters@0: 	    case Inverse:       x = 1.0/x;                  break;
tomwalters@0: 	    case Round:         x = round(x);               break;
tomwalters@0: 	    case Absolute:      x = fabs(x);                break;
tomwalters@0: 	    case Power:         x = pow(x,op[i].arg);       break;
tomwalters@0: 	    case Exp:           x = pow(op[i].arg,x);       break;
tomwalters@0: 	    case Eexp:          x = exp(x);                 break;
tomwalters@0: 	    case Log:           x = log(x)/log(op[i].arg);  break;
tomwalters@0: 	    case Elog:          x = log(x);                 break;
tomwalters@0: 	    case Sin:           x = sin(x);                 break;
tomwalters@0: 	    case Cos:           x = cos(x);                 break;
tomwalters@0: 	    case Tan:           x = tan(x);                 break;
tomwalters@0: 	    case Threshold:     x = threshold(x,op[i].arg); break;
tomwalters@0: 	    case Diff1:         x = diff1(x);               break;
tomwalters@0: 	    case Diff2:         x = diff2(x);               break;
tomwalters@0: 	    case Diff3:         x = diff3(x);               break;
tomwalters@0: 	    case Sum:           x = sum(x);                 break;
tomwalters@0: 	    case Mean:          x = mean(x);                break;
tomwalters@0: 	    case Min:           x = min(x);                 break;
tomwalters@0: 	    case Max:           x = max(x);                 break;
tomwalters@0: 	}
tomwalters@0:     return x;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: match(s1,s2)
tomwalters@0: char *s1, *s2;
tomwalters@0: {
tomwalters@0:     int  n1, n2;
tomwalters@0: 
tomwalters@0:     if ((n1=strlen(s1)) < (n2=strlen(s2)))
tomwalters@0: 	return strncmp(s1,s2,n1);
tomwalters@0:     else
tomwalters@0: 	return strncmp(s1,s2,n2);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:  The past values used to compute differences are, (in the order they would
tomwalters@0:  appear in an array, where x is the current value):
tomwalters@0: 	1st diffs           ... x11 x
tomwalters@0: 	2nd diffs       ... x22 x21 x
tomwalters@0: 	3rd diffs   ... x33 x32 x31 x
tomwalters@0: *****************************************************************************/
tomwalters@0: float diff1(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float x11=0;
tomwalters@0:     float  d;
tomwalters@0: 
tomwalters@0:     d=fabs(x11-x);
tomwalters@0:     x11=x;
tomwalters@0:     return d;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: float diff2(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float x21=0, x22=0;
tomwalters@0:     float  d;
tomwalters@0: 
tomwalters@0:     d=fabs( fabs(x22-x21)-fabs(x21-x) );
tomwalters@0:     x22=x21;
tomwalters@0:     x21=x;
tomwalters@0:     return d;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: float diff3(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float x31=0, x32=0, x33=0;
tomwalters@0:     float  d;
tomwalters@0: 
tomwalters@0:     d=fabs( fabs( fabs(x33-x32)-fabs(x32-x31) ) - fabs( fabs(x32-x31)-fabs(x31-x) ) );
tomwalters@0:     x33=x32;
tomwalters@0:     x32=x31;
tomwalters@0:     x31=x;
tomwalters@0:     return d;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:  Recurisve estimation of the mean
tomwalters@0:  (Ref: Young, eqn 2.14 (p14), eqn VI(2) (p63), eqn 5.17 (p65))
tomwalters@0: ****************************************************************************/
tomwalters@0: float mean(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float a = 0;         /* initial mean value */
tomwalters@0:     static float p = 10000;     /* initial gain factor [Young p65] */
tomwalters@0: 
tomwalters@0:     p = p/(1+p);                /* update p (the gain factor)   */
tomwalters@0:     a = a - p*(a-x);            /* update a (the mean estimate) */
tomwalters@0:     return a;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:  Min, Max, and Sum
tomwalters@0: ****************************************************************************/
tomwalters@0: float min(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float MIN=1.0e30;
tomwalters@0: 
tomwalters@0:     if (x < MIN) MIN = x;
tomwalters@0:     return MIN;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: float max(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float MAX=0;
tomwalters@0: 
tomwalters@0:     if (x > MAX) MAX = x;
tomwalters@0:     return MAX;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: float sum(x)
tomwalters@0: float x;
tomwalters@0: {
tomwalters@0:     static float SUM=0;
tomwalters@0: 
tomwalters@0:     SUM += x;
tomwalters@0:     return SUM;
tomwalters@0: }