tomwalters@0: /*
tomwalters@0:     merge.c     Merge N streams onto stdout using given operator.
tomwalters@0:     -------
tomwalters@0: 
tomwalters@0:  Usage:  merge  [options]  file1  file2  ...  fileN
tomwalters@0: 
tomwalters@0:  For each vector of items (p1,p2,...,pN) read from file1, file2,...,fileN
tomwalters@0:  respectively write the result of a merging operation on the stdout.
tomwalters@0: 
tomwalters@0:  An input item is a number depending upon the type option.
tomwalters@0:  The output item has the same type as that selected for input.
tomwalters@0: 
tomwalters@0:  Internal processing is in floating point. A scale factor is included to
tomwalters@0:  avoid overflow when casting output data.
tomwalters@0:  If 16-bit over or under-flow occurs a warning is printed on the stderr.
tomwalters@0: 
tomwalters@0:  The `range' option sets the start and duration of the process.
tomwalters@0:  Its arguments are of the form: range=a-b (where start=a and duration=b-a+1)
tomwalters@0: 			    or: range=a   (where start=a and duration=1    )
tomwalters@0:  The arguments can be in time units (ms, s) or samples (no units), and both
tomwalters@0:  "min" (start of file) and "max" (end of file) are recognised.
tomwalters@0:  Samples in a file are numbered 0,1,2,...,max.
tomwalters@0: 
tomwalters@0:  All files are interpreted as functions sampled with origin at the point in
tomwalters@0:  the file given by the first argument of the `range' option.
tomwalters@0:  The functions are assumed to be zero for all samples beyond those specified
tomwalters@0:  in the files.
tomwalters@0: 
tomwalters@0:  The `phase' option shifts the origin to introduce a phase difference between
tomwalters@0:  files. The `phase' option takes as argument a comma-separated list of
tomwalters@0:  arguments which must be the same length as the number of file arguments.
tomwalters@0:  Each argument can be in time units (ms, s) or samples (no units).
tomwalters@0:  A +ve phase advances the origin into the file.
tomwalters@0:  A -ve phase retards the origin. If a negative phase in combination with
tomwalters@0:  the start of the file as set by the `range' option results in a new origin
tomwalters@0:  outside the specified file, then the difference is padded with zeroes.
tomwalters@0: 
tomwalters@0:  The weights option takes as argument a comma-separated list of real numbers
tomwalters@0:  which must be the same length as the number of file arguments.
tomwalters@0:  These number are applied as weights to each item prior to the merging
tomwalters@0:  operation.
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: #include "units.h"
tomwalters@0: 
tomwalters@0: char applic[] = "Merge N streams onto stdout using given operator." ;
tomwalters@0: char usage[]  = "merge [options] file1 file2 ... fileN" ;
tomwalters@0: 
tomwalters@0: static char *helpstr,   *debugstr,  *sampstr,   *opstr      ;
tomwalters@0: static char *phstr,     *rangestr,  *scalestr,  *typestr    ;
tomwalters@0: static char *wtstr ;
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:     {   "operator"  ,   "mean"      ,  &opstr       ,   "operator for merging"                  , VAL     },
tomwalters@0:     {   "phase"     ,   "off"       ,  &phstr       ,   "phase shift (phi1,...,phiN)"           , VAL     },
tomwalters@0:     {   "weights"   ,   "off"       ,  &wtstr       ,   "weight coefficients (w1,...,wN)"       , VAL     },
tomwalters@0:     {   "range"     ,   "0-max"     ,  &rangestr    ,   "inclusive time range for all files"    , VAL     },
tomwalters@0:     {   "scale"     ,   "1.0"       ,  &scalestr    ,   "scale factor for output"               , VAL     },
tomwalters@0:     {   "type"      ,   "short"     ,  &typestr     ,   "datatype"                              , VAL     },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: char *operator[] = {
tomwalters@0:     "cat"      ,
tomwalters@0:     "add"      ,
tomwalters@0:     "subtract" ,
tomwalters@0:     "multiply" ,
tomwalters@0:     "divide"   ,
tomwalters@0:     "abs"      ,
tomwalters@0:     "max"      ,
tomwalters@0:     "min"      ,
tomwalters@0:     "mean"     ,
tomwalters@0:     "norm"     ,
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #define Cat          0  /* Functions (index of operator list) */
tomwalters@0: #define Add          1
tomwalters@0: #define Subtract     2
tomwalters@0: #define Multiply     3
tomwalters@0: #define Divide       4
tomwalters@0: #define Abs          5
tomwalters@0: #define Max          6
tomwalters@0: #define Min          7
tomwalters@0: #define Mean         8
tomwalters@0: #define Norm         9
tomwalters@0: 
tomwalters@0: #define sq(p)           ((p)*(p))
tomwalters@0: 
tomwalters@0: int     samplerate ;
tomwalters@0: int     bytes      ;
tomwalters@0: int     type       ;    /* datatype index */
tomwalters@0: int     op         ;    /* operator index */
tomwalters@0: float   scale      ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: main(argc, argv)
tomwalters@0: int    argc;
tomwalters@0: char **argv;
tomwalters@0: {
tomwalters@0:     FILE  **fp ;
tomwalters@0:     float  *vec, *weight ;
tomwalters@0:     int    *origin ;
tomwalters@0:     char  **list ;
tomwalters@0:     int     i, j, n, morehelp() ;
tomwalters@0:     int     a, b ;
tomwalters@0:     float   writevec(), f, newscale = 1. ;
tomwalters@0: 
tomwalters@0:     if ( ( i = getopts( option, argc, argv ) ) == 0 || !isoff( helpstr ) )
tomwalters@0: 	helpopts2( helpstr, argv[0], applic, usage, option, morehelp ) ;
tomwalters@0: 
tomwalters@0:     samplerate  = to_Hz( sampstr, 0 ) ;
tomwalters@0:     scale       = atof( scalestr ) ;
tomwalters@0:     op          = opindex( opstr ) ;
tomwalters@0: 
tomwalters@0:     if ( ( type = typeindex( typestr ) ) < 0 ) {
tomwalters@0: 	fprintf( stderr, "merge: bad type [%s]\n", typestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     bytes = typebytes( type ) ;
tomwalters@0: 
tomwalters@0:     if ( range( rangestr, &a, &b, samplerate ) == 0 ) {
tomwalters@0: 	fprintf(stderr,"merge: bad range [%s]\n", rangestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     fp     = (FILE **)malloc( i * sizeof( FILE *) ) ;
tomwalters@0:     vec    = (float *)malloc( i * sizeof( float ) ) ;
tomwalters@0:     origin = (int   *)malloc( i * sizeof( int   ) ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( isoff( wtstr ) ) weight = (float *)0 ;
tomwalters@0:     else {
tomwalters@0: 	list   = (char **)malloc( i * sizeof( char *) ) ;
tomwalters@0: 	weight = (float *)malloc( i * sizeof( float ) ) ;
tomwalters@0: 	if ( ( n = tokens( wtstr, list, i, ',' ) ) != i ) {
tomwalters@0: 	    fprintf( stderr,"merge: incorrect number of weights\n" ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 	for ( j = 0 ; j < n ; j++ )
tomwalters@0: 	    weight[j] = atof( list[j] ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( isoff( phstr ) )
tomwalters@0: 	for ( j = 0 ; j < i ; j++ )
tomwalters@0: 	    origin[j] = a ;
tomwalters@0:     else {
tomwalters@0: 	list   = (char **)malloc( i * sizeof( char *) ) ;
tomwalters@0: 	if ( ( n = tokens( phstr, list, i, ',' ) ) != i ) {
tomwalters@0: 	    fprintf( stderr,"merge: incorrect number of phases\n" ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 	for ( j = 0 ; j < n ; j++ )
tomwalters@0: 	    origin[j] = a + to_p( list[j], samplerate ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     for ( n = 0 ; i > 0 ; i--, n++ )  {
tomwalters@0: 	if ( ( fp[n] = fropen( argv[argc-i] ) ) == (FILE *)0 ) {
tomwalters@0: 	    fprintf( stderr,"merge: can't open %s\n", argv[argc-i] ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	if ( origin[i] > 0 )
tomwalters@0: 	    if ( seekstart( origin[i], bytes, fp[i] ) < origin[i] ) {
tomwalters@0: 		fprintf( stderr, "insufficient data in file%d\n", i+1 ) ;
tomwalters@0: 		exit( 1 ) ;
tomwalters@0: 	    }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     for ( i = a ; ( i <= b || b == (-1) ) && readvec( vec, n, fp, type, origin ) ; i++ ) {
tomwalters@0: 
tomwalters@0: 	if ( weight != (float *)0 )
tomwalters@0: 	    do_weights( vec, weight, n ) ;
tomwalters@0: 
tomwalters@0: 	j = do_operation( vec, n, op ) ;
tomwalters@0: 	if ( ( f = writevec( vec, j, type, scale ) ) < newscale )
tomwalters@0: 	    newscale = f ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( newscale < 1. )
tomwalters@0: 	fprintf( stderr, "Warning: 16-bit overflow during merge. Try scale<%.4f\n", newscale ) ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Return the operator index (to the operator list) of the given operator string.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: int opindex( opstr )
tomwalters@0: char *opstr ;
tomwalters@0: {
tomwalters@0:     int  i ;
tomwalters@0: 
tomwalters@0:     if ( ( i = listindex( operator, opstr ) ) < 0 ) {
tomwalters@0: 	if ( i == (-1) )
tomwalters@0: 	    fprintf(stderr,"merge: unknown operator [%s]\n", opstr ) ;
tomwalters@0: 	if ( i == (-2) )
tomwalters@0: 	    fprintf(stderr,"merge: ambiguous operator [%s]\n", opstr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     return i ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Assign the given n-vector of floats using one item of given type (an index
tomwalters@0: to the datatype list in options.h) read from each of the n streams given in
tomwalters@0: the fp array. Return 0 if eof on any stream, otherwise return 1.
tomwalters@0: If the stream origin < 0 read a value of 0, increment the origin, and return
tomwalters@0: 1 (successful read).
tomwalters@0: */
tomwalters@0: 
tomwalters@0: readvec( vec, n, fp, type, origin )
tomwalters@0: float  *vec  ;
tomwalters@0: int     n    ;
tomwalters@0: FILE  **fp   ;
tomwalters@0: int     type ;
tomwalters@0: int    *origin ;
tomwalters@0: {
tomwalters@0:     int  i ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 
tomwalters@0: 	if ( origin[i] < 0 ) {
tomwalters@0: 	    vec[i] = 0  ;
tomwalters@0: 	    origin[i]++ ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	else if ( readitem( &vec[i], type, 1, fp[i] ) == 0 )
tomwalters@0: 	    return 0 ;
tomwalters@0: 
tomwalters@0:     }
tomwalters@0:     return 1 ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Write n elements from the given vector on the stdout in the given type
tomwalters@0: (an index to the datatype list in options.h)
tomwalters@0: Multiply each element by the given scale factor (prior to type conversion).
tomwalters@0: If any element will over or underflow when scaled and cast into type
tomwalters@0: then return a more appropriate scale factor.
tomwalters@0: Return scale factor of 1 if no over or underflow.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: float writevec( vec, n, type, scale )
tomwalters@0: float *vec  ;
tomwalters@0: int    n    ;
tomwalters@0: int    type ;
tomwalters@0: float  scale;
tomwalters@0: {
tomwalters@0:     float newscale ;
tomwalters@0:     int   i ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 
tomwalters@0: 	newscale = check_overflow( vec[i], scale, type ) ;
tomwalters@0: 
tomwalters@0: 	vec[i] *= scale ;
tomwalters@0: 	writeitem( &vec[i], type, 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     return ( newscale ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Perform the given operation (op is an index to the operations list) on the
tomwalters@0: given n-vector, storing the results in the vector.
tomwalters@0: Return the size of the result (the number of elements of vec it occupies).
tomwalters@0: */
tomwalters@0: 
tomwalters@0: do_operation( vec, n, op )
tomwalters@0: float *vec  ;
tomwalters@0: int    n    ;
tomwalters@0: int    op   ;
tomwalters@0: {
tomwalters@0:     int  i ;
tomwalters@0: 
tomwalters@0:     switch ( op ) {
tomwalters@0: 
tomwalters@0: 	case Cat      : return n ;
tomwalters@0: 	case Add      : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] += vec[i] ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Subtract : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] -= vec[i] ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Abs      : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] = fabs(vec[0]-vec[i]) ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Multiply : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] *= vec[i] ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Divide   : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] /= vec[i] ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Max      : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    if ( vec[i] > vec[0] ) vec[0] = vec[i] ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Min      : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    if ( vec[i] < vec[0] ) vec[0] = vec[i] ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Mean     : for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] += vec[i] ;
tomwalters@0: 			vec[0] /= n ;
tomwalters@0: 			return 1 ;
tomwalters@0: 	case Norm:      vec[0] = sq(vec[0]) ;
tomwalters@0: 			for ( i = 1 ; i < n ; i++ )
tomwalters@0: 			    vec[0] += sq(vec[i]) ;
tomwalters@0: 			vec[0] = sqrt(vec[0]) ;
tomwalters@0: 			return 1 ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: do_weights( vec, weight, n )
tomwalters@0: float *vec    ;
tomwalters@0: float *weight ;
tomwalters@0: int    n      ;
tomwalters@0: {
tomwalters@0:     int  i ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	vec[i] *= weight[i] ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: morehelp()
tomwalters@0: {
tomwalters@0:     fprintf(stderr,"\noperators:                       \n");
tomwalters@0:     fprintf(stderr,"            cat         concatenate  p1,p2,p3,... \n");
tomwalters@0:     fprintf(stderr,"            add         p1+p2+p3+...      \n");
tomwalters@0:     fprintf(stderr,"            subtract    p1-p2-p3-...      \n");
tomwalters@0:     fprintf(stderr,"            abs         ||p1-p2|-p3|-...    \n");
tomwalters@0:     fprintf(stderr,"            multiply    p1*p2*p3*...      \n");
tomwalters@0:     fprintf(stderr,"            divide      p1/p2/p3/...      \n");
tomwalters@0:     fprintf(stderr,"            max         max(p1,p2,p3,...) \n");
tomwalters@0:     fprintf(stderr,"            min         min(p1,p2,p3,...) \n");
tomwalters@0:     fprintf(stderr,"            mean        mean(p1,p2,p3,...) \n");
tomwalters@0:     fprintf(stderr,"            norm        sqrt( p1*p1+p2*p2+p3*p3+... )  \n");
tomwalters@0:     exit( 1 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: