tomwalters@0: /*
tomwalters@0:   stats.c      calculate some statistics (min, max, mean etc.)
tomwalters@0:   -------
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "options.h"
tomwalters@0: #include "units.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: 
tomwalters@0: char applic[]     = "print statistics of input files" ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *typestr, *statstr ;
tomwalters@0: static char *sampstr, *frstr,    *widstr,  *sizestr ;
tomwalters@0: static char *linestr, *precstr,  *fieldstr          ;
tomwalters@0: 
tomwalters@0: static Options option[] = {
tomwalters@0:     {   "help"      ,   "off"       ,  &helpstr     ,   "help"                          , DEBUG   },
tomwalters@0:     {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"              , DEBUG   },
tomwalters@0:     {   "samplerate",   "20kHz"     ,  &sampstr     ,   "samplerate "                   , VAL     },
tomwalters@0:     {   "frame"     ,   "1"         ,  &frstr       ,   "select frames inclusively"     , VAL     },
tomwalters@0:     {   "width"     ,   "max"       ,  &widstr      ,   "frame width"                   , VAL     },
tomwalters@0:     {   "type"      ,   "short"     ,  &typestr     ,   "data type"                     , VAL     },
tomwalters@0:     {   "stat"      ,   "mean"      ,  &statstr     ,   "statistic list (comma separated)"      , VAL     },
tomwalters@0:     {   "line"      ,   "off"       ,  &linestr     ,   "force ascii output"            , SETFLAG },
tomwalters@0:     {   "fieldwidth",   "10"        ,  &fieldstr    ,   "field width for ascii printing", SVAL    },
tomwalters@0:     {   "precision" ,   "3"         ,  &precstr     ,   "precision for ascii printing"  , SVAL    },
tomwalters@0:     {   "SIZE"      ,   "262144p"   ,  &sizestr     ,   "buffer size (s, ms, or p)"     , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: typedef float (*PF)() ; /* pointer-to-function type returning float */
tomwalters@0: 
tomwalters@0: typedef struct {        /* a struct to hold both function name and pointer */
tomwalters@0:     char  *name ;       /* string name of function */
tomwalters@0:     PF     func ;       /* pointer to function     */
tomwalters@0:     char  *help ;       /* explanation of use      */
tomwalters@0: } Func ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float   num(),      sum(),      sumabs(),   sumsq(),    mean()      ;
tomwalters@0: float   rms(),      variance(), stddev(),   min(),      max()       ;
tomwalters@0: float   absmax(),   absrange(), pc()                                            ;
tomwalters@0: 
tomwalters@0: float   call_by_name() ;
tomwalters@0: 
tomwalters@0: Func  statistic[] = {
tomwalters@0:     { "n"         ,     num         ,    "sample size (ie. frame width)"                        } ,
tomwalters@0:     { "sum"       ,     sum         ,    "area"                                                 } ,
tomwalters@0:     { "sa"        ,     sumabs      ,    "sum of absolute values"                               } ,
tomwalters@0:     { "ss"        ,     sumsq       ,    "sum of squares"                                       } ,
tomwalters@0:     { "mean"      ,     mean        ,    "sum/n"                                                } ,
tomwalters@0:     { "rms"       ,     rms         ,    "root mean square: square root of ss/n"                } ,
tomwalters@0:     { "variance"  ,     variance    ,    "computed with n-1 degrees of freedom"                 } ,
tomwalters@0:     { "stddev"    ,     stddev      ,    "standard deviation: square root of variance"          } ,
tomwalters@0:     { "min"       ,     min         ,    "minimum value"                                        } ,
tomwalters@0:     { "max"       ,     max         ,    "maximum value"                                        } ,
tomwalters@0:     { "absmax"    ,     absmax      ,    "maximum absolute value, ie. max( |max|, |min| )"      } ,
tomwalters@0:     { "range"     ,     absrange    ,    "max-min"                                              } ,
tomwalters@0:     { "pc_1of2"   ,     pc          ,    "percentage correct, ie. greater first of each pair"   } ,
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int    samplerate ;
tomwalters@0: int    bytes      ;
tomwalters@0: int    type       ;     /* datatype index */
tomwalters@0: int    nstat      ;     /* number of statistics required */
tomwalters@0: int    width      ;
tomwalters@0: int    SIZE       ;
tomwalters@0: int    lineout    ;     /* flag for forcing ascii output */
tomwalters@0: 
tomwalters@0: float *buf ;            /* args to every stat function */
tomwalters@0: int    n   ;
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int    argc;
tomwalters@0: char **argv;
tomwalters@0: {
tomwalters@0:     FILE   *fp   ;
tomwalters@0:     char  **list ;
tomwalters@0:     float   y    ;
tomwalters@0:     int     i, j, k, a, b, helpstats() ;
tomwalters@0: 
tomwalters@0:     i = getopts( option, argc, argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts1( helpstr, argv[0], applic, option, helpstats ) ;
tomwalters@0: 
tomwalters@0:     if ( ( type = typeindex( typestr ) ) < 0 ) {
tomwalters@0: 	fprintf( stderr, "stats: bad type [%s]\n", typestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     bytes = typebytes( type ) ;
tomwalters@0: 
tomwalters@0:     if ( selector( frstr, &a, &b ) == 0 ) {
tomwalters@0: 	fprintf(stderr,"stats: bad frame selector [%s]\n", frstr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     PRECISION  = atoi( precstr ) ;
tomwalters@0:     FIELDWIDTH = atoi( fieldstr ) ;
tomwalters@0:     lineout    = ison( linestr ) ;
tomwalters@0:     samplerate = to_Hz( sampstr, 0) ;
tomwalters@0:     SIZE       = to_p( sizestr, samplerate ) ;
tomwalters@0:     if ( ismax( widstr ) ) width = (-1) ;
tomwalters@0:     else                   width = to_p( widstr, samplerate ) ;
tomwalters@0: 
tomwalters@0:     nstat  = flistsize( statistic ) ;   /* max number of statistics */
tomwalters@0:     list   = (char **)malloc( nstat * sizeof( char *) ) ;
tomwalters@0:     if ( ( nstat  = tokens( statstr, list, nstat, ',' ) ) == 0 ) {
tomwalters@0: 	fprintf(stderr,"stats: incorrect stat list\n" ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     do {
tomwalters@0: 
tomwalters@0: 	if ( i == 0 ) fp = stdin ;
tomwalters@0: 	else if ( ( fp = fopen( argv[argc-i], "r" ) ) == (FILE *)0 ) {
tomwalters@0: 	    fprintf( stderr,"stats: can't open %s\n", argv[argc-i] ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	if ( width > 0 ) buf = (float *)malloc( width * sizeof(float) ) ;
tomwalters@0: 	else             buf = (float *)malloc( SIZE  * sizeof(float) ) ;
tomwalters@0: 
tomwalters@0: 	if ( width > 0 && a > 1 )
tomwalters@0: 	    if ( seekstart( (a-1)*width, bytes, fp ) < (a-1)*width )
tomwalters@0: 		fprintf( stderr, "stats warning: insufficient input\n" ) ;
tomwalters@0: 
tomwalters@0: 	for ( j = a ; ( j <= b || b == 0 ) && ( n = Readitem( buf, type, width, fp, SIZE ) ) > 0 ; j++ ) {
tomwalters@0: 
tomwalters@0: 	    for ( k = 0 ; k < nstat ; k++ ) {
tomwalters@0: 
tomwalters@0: 		y = call_by_name( list[k], statistic ) ;
tomwalters@0: 
tomwalters@0: 		if ( lineout )
tomwalters@0: 		    printf( "%*.*f  ", FIELDWIDTH, PRECISION, y ) ;
tomwalters@0: 		else
tomwalters@0: 		    writeitem( &y, type, 1, stdout ) ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	free( buf ) ;
tomwalters@0: 	fclose( fp ) ;
tomwalters@0: 	if ( lineout ) printf( "\n" ) ;
tomwalters@0: 
tomwalters@0:     } while ( --i > 0 ) ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Return the number of structs in the null-terminated array of Func structs.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: flistsize( flist )
tomwalters@0: Func *flist ;
tomwalters@0: {
tomwalters@0:     int  i ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; flist[i].name != (char *)0 ; i++ )
tomwalters@0: 	;
tomwalters@0:     return i ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Call a function given its string name. The function's address is found against
tomwalters@0: the matching name in the given `flist'.
tomwalters@0: The function returns an float.
tomwalters@0: Rather than return at the first matching name, search all the list for names
tomwalters@0: which match the (possibly abbreviated) string name. Print a warning if
tomwalters@0: ambiguity found.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: float call_by_name( str, flist )
tomwalters@0: char  *str   ;
tomwalters@0: Func  *flist ;
tomwalters@0: {
tomwalters@0:     int  i, j = (-1) ;
tomwalters@0: 
tomwalters@0:     for ( i=0 ; flist[i].name != (char *)0 ; i++ ) {
tomwalters@0: 	if ( iststr( str, flist[i].name ) ) {
tomwalters@0: 	    if ( j >= 0 )
tomwalters@0: 		fprintf( stderr,"warning: ambiguous stat names (%s and %s)\n", flist[j].name, flist[i].name ) ;
tomwalters@0: 	    j = i ;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0:     return ( (*flist[j].func)() ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Read items of given type up to max of SIZE items (when n=-1).
tomwalters@0: Return n or 0 if eof.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: Readitem( y, type, n, fp, SIZE )
tomwalters@0: float  *y    ;
tomwalters@0: int     type ;
tomwalters@0: int     n    ;
tomwalters@0: FILE   *fp   ;
tomwalters@0: int     SIZE ;
tomwalters@0: {
tomwalters@0:     if ( n == (-1) ) {
tomwalters@0: 	for ( n = 0 ; n < SIZE && readitem( &y[n], type, 1, fp ) ; n++ )
tomwalters@0: 	    ;
tomwalters@0: 	if ( n == SIZE )
tomwalters@0: 	    fprintf( stderr, "stats warning: buffer full\n" ) ;
tomwalters@0:     }
tomwalters@0:     else if ( readitem( y, type, n, fp ) == 0 ) return 0 ;
tomwalters@0: 
tomwalters@0:     return n ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Print help on request
tomwalters@0: */
tomwalters@0: 
tomwalters@0: helpstats()
tomwalters@0: {
tomwalters@0:     int  i ;
tomwalters@0: 
tomwalters@0:     fprintf(stderr,"\nstatistics:       \n");
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; statistic[i].name != (char *)0 ; i++ )
tomwalters@0: 	fprintf(stderr,"            %-10s  %s\n", statistic[i].name, statistic[i].help ) ;
tomwalters@0: 
tomwalters@0:     exit( 1 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*************************** functions ************************************/
tomwalters@0: /* The functions all take two args: char *buf; int n;                     */
tomwalters@0: /* These are declared externally.                                         */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float num()
tomwalters@0: {
tomwalters@0:     return n ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float sum()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float sum = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	sum += buf[i] ;
tomwalters@0: 
tomwalters@0:     return ( sum ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float sumabs()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float sum = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	sum += fabs( (double)buf[i] ) ;
tomwalters@0: 
tomwalters@0:     return ( sum ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float sumsq()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float sumsq = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	sumsq += ( buf[i] * buf[i] ) ;
tomwalters@0: 
tomwalters@0:     return ( sumsq ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float mean()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float sum = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	sum += buf[i] ;
tomwalters@0: 
tomwalters@0:     return ( sum / n ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float rms()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float sumsq = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	sumsq += ( buf[i] * buf[i] ) ;
tomwalters@0: 
tomwalters@0:     return ( sqrt( sumsq / n ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float variance()    /* computed with n-1 degrees of freedom */
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float sum = 0, sumsq = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 	sum   += buf[i]   ;
tomwalters@0: 	sumsq += ( buf[i] * buf[i] ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     return ( ( sumsq - sum*sum/n ) / ( n - 1 ) ) ;  /* replace n-1 by n for n degrees of freedom */
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float stddev()    /* computed with n-1 degrees of freedom */
tomwalters@0: {
tomwalters@0:     return ( sqrt( variance( buf, n ) ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float min()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float min = 99999999999. ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	if ( buf[i] < min ) min = buf[i] ;
tomwalters@0:     return ( min ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float max()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float max = ( -99999999999. ) ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	if ( buf[i] > max ) max = buf[i] ;
tomwalters@0: 
tomwalters@0:     return ( max ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float absmax()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float min =    99999999999.   ;
tomwalters@0:     float max = ( -99999999999. ) ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 	if ( buf[i] < min ) min = buf[i] ;
tomwalters@0: 	if ( buf[i] > max ) max = buf[i] ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( ( max = fabs( (double)max ) ) < ( min = fabs( (double)min ) ) )
tomwalters@0: 	max = min ;
tomwalters@0: 
tomwalters@0:     return ( max ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float absrange()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float min =    99999999999.   ;
tomwalters@0:     float max = ( -99999999999. ) ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 	if ( buf[i] < min ) min = buf[i] ;
tomwalters@0: 	if ( buf[i] > max ) max = buf[i] ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     return ( max - min ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: The buf contains n/2 pairs of numbers. Calculate how many of each pair have
tomwalters@0: the first number greater than the last, and return this as a "percentage
tomwalters@0: correct".
tomwalters@0: */
tomwalters@0: 
tomwalters@0: float pc()
tomwalters@0: {
tomwalters@0:     int   i ;
tomwalters@0:     float p = 0 ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i+=2 )
tomwalters@0: 	if ( buf[i] > buf[i+1] )
tomwalters@0: 	    p++ ;
tomwalters@0: 
tomwalters@0:     return ( ( p / (n/2) ) * 100 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: