tomwalters@0: /*
tomwalters@0:   scale.c       shift and scale a data stream by given constants or to
tomwalters@0:   -------       fit a given range (top - bottom).
tomwalters@0: 
tomwalters@0:   Shifting and scaling is done for all input x using the formula:
tomwalters@0: 	y = ( x - shift ) * scale
tomwalters@0: 
tomwalters@0:   If only a shift term   is given, then the default scale is 1.0
tomwalters@0:   If only a scale factor is given, then the default shift is 0.
tomwalters@0: 
tomwalters@0:   If neither shift nor scale are given, then appropriate values are
tomwalters@0:   chosen so that the output fits the range specified by `top' and `bottom'.
tomwalters@0:   For an input which ranges between extreme values of `max' and `min', the
tomwalters@0:   shift and scale factors to give a required range between `top' and `bottom'
tomwalters@0:   are found by simultaneously solving:
tomwalters@0: 	top    = ( max - shift ) * scale
tomwalters@0: 	bottom = ( min - shift ) * scale
tomwalters@0: 
tomwalters@0:   If top < bottom then the input is inverted in the resulting range.
tomwalters@0: 
tomwalters@0:   The input stream depends upon the `type' option. The output datatype is
tomwalters@0:   the same as that selected for input.
tomwalters@0:   With no filename arguments, data is expected on the stdin, and the
tomwalters@0:   scaled result is written on the stdout. Otherwise each given filename is
tomwalters@0:   processed in turn. When the `output' option is "off" output overwrites the
tomwalters@0:   respective input file. (The respective output for the stdin is the stdout).
tomwalters@0:   Otherwise all scaled input files are written to the given output file
tomwalters@0:   (which is the stdout by default).
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:   Examples:
tomwalters@0: 
tomwalters@0: 1. Invert the range of a given pulse train to obtain negative-going pulses.
tomwalters@0:    (Two ways):
tomwalters@0: 
tomwalters@0: ptrain amp=512 | scale top=-512 bot=0 | x11plot
tomwalters@0: ptrain amp=512 | scale scale=-1       | x11plot
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[]     = "waveform shifting and scaling:  y = (x - shift) * scale" ;
tomwalters@0: 
tomwalters@0: static char *helpstr,   *debugstr,  *sampstr,   *rangestr,  *outstr     ;
tomwalters@0: static char *topstr,    *botstr,    *scalestr,  *typestr,   *sizestr    ;
tomwalters@0: static char *shiftstr,  *normstr                                        ;
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:     {   "range"     ,   "0-max"     ,  &rangestr    ,   "start-finish limits in data"   , VAL     },
tomwalters@0:     {   "top"       ,   "1000"      ,  &topstr      ,   "max of scaled data"            , VAL     },
tomwalters@0:     {   "bottom"    ,   "-1000"     ,  &botstr      ,   "min of scaled data"            , VAL     },
tomwalters@0:     {   "shift"     ,   "off"       ,  &shiftstr    ,   "shift term"                    , VAL     },
tomwalters@0:     {   "scale"     ,   "off"       ,  &scalestr    ,   "scale factor"                  , VAL     },
tomwalters@0:     {   "normalize" ,   "off"       ,  &normstr     ,   "set zero mean and unit std dev", VAL     },
tomwalters@0:     {   "type"      ,   "short"     ,  &typestr     ,   "datatype"                      , VAL     },
tomwalters@0:     {   "output"    ,   "stdout"    ,  &outstr      ,   "output filename (off = overwrite input)"   , VAL     },
tomwalters@0:     {   "SIZE"      ,   "262144p"   ,  &sizestr     ,   "buffer size (s, ms, or p)"     , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int    samplerate  ;
tomwalters@0: int    type        ;    /* datatype index */
tomwalters@0: int    bytes       ;
tomwalters@0: int    SIZE        ;    /* buffer size    */
tomwalters@0: 
tomwalters@0: int    FIT_TO_RANGE = 0 ;
tomwalters@0: int    NORMALIZE    = 0 ;
tomwalters@0: 
tomwalters@0: float  shift, scale ;
tomwalters@0: 
tomwalters@0: float *data  ;
tomwalters@0: 
tomwalters@0: FILE  *ofp   ;
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int     argc;
tomwalters@0: char  **argv;
tomwalters@0: {
tomwalters@0:     FILE   *fp ;
tomwalters@0:     int     i, j, a, b, n ;
tomwalters@0: 
tomwalters@0:     i = getopts( option, argc, argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts( helpstr, argv[0], applic, option ) ;
tomwalters@0: 
tomwalters@0:     samplerate = to_Hz( sampstr ) ;
tomwalters@0:     SIZE       = to_p( sizestr, samplerate ) ;
tomwalters@0: 
tomwalters@0:     if ( ( type = typeindex( typestr ) ) < 0 ) {
tomwalters@0: 	fprintf( stderr, "scale: 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,"scale: bad range [%s]\n", rangestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( ison( normstr ) )  NORMALIZE = 1 ;
tomwalters@0:     else if ( isoff( shiftstr ) && isoff( scalestr ) ) FIT_TO_RANGE = 1 ;
tomwalters@0:     else if ( isoff( shiftstr ) ) {
tomwalters@0: 	shift = 0  ;
tomwalters@0: 	scale = atof( scalestr ) ;
tomwalters@0:     }
tomwalters@0:     else if ( isoff( scalestr ) ) {
tomwalters@0: 	shift = atof( shiftstr ) ;
tomwalters@0: 	scale = 1. ;
tomwalters@0:     }
tomwalters@0:     else {
tomwalters@0: 	shift = atof( shiftstr ) ;
tomwalters@0: 	scale = atof( scalestr ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( ( data = (float *)malloc( SIZE * sizeof(float) ) ) == NULL ) {
tomwalters@0: 	fprintf( stderr, "malloc out of space\n" ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
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,"scale: can't open %s\n", argv[argc-i] ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	if ( seekstart( a, bytes, fp ) < a ) {
tomwalters@0: 	    fprintf( stderr,"scale: insufficient data in file\n" ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	for ( j = a, n = 0 ; ( j <= b || b == (-1) ) && n < SIZE && readitem( &data[n], type, 1, fp ) ; j++, n++ )
tomwalters@0: 	    ;
tomwalters@0: 	if ( n == SIZE )
tomwalters@0: 	    fprintf( stderr, "scale warning: file %s exceeded buffer size\n",  argv[argc-i] ) ;
tomwalters@0: 
tomwalters@0: 	fclose( fp ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 	if ( NORMALIZE )
tomwalters@0: 	    getnorm( data, n, &shift, &scale ) ;
tomwalters@0: 	else if ( FIT_TO_RANGE )
tomwalters@0: 	    getrange( data, n, atof( topstr ), atof( botstr ), &shift, &scale ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 	if ( isstr( outstr, "print" ) )
tomwalters@0: 	    printf( "%8.3f  %8.3f \n", shift, scale ) ;
tomwalters@0: 	else {
tomwalters@0: 	    open_output( argc, argv , i ) ;
tomwalters@0: 	    for ( j = 0 ; j < n ; j++ ) {
tomwalters@0: 		data[j] = ( data[j] - shift ) * scale ;
tomwalters@0: 		writeitem( &data[j], type, 1, ofp ) ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     } while ( --i > 0 ) ;
tomwalters@0: 
tomwalters@0:     fclose( ofp ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Open ouput file pointer, ofp.
tomwalters@0: Use the input filename (argv[argc-i]) if outstr is "off", otherwise use
tomwalters@0: the given name outstr.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: open_output( argc, argv , i )
tomwalters@0: int     argc ;
tomwalters@0: char  **argv ;
tomwalters@0: int     i    ;
tomwalters@0: {
tomwalters@0:     static int first = 1 ;
tomwalters@0: 
tomwalters@0:     if ( first ) {
tomwalters@0: 
tomwalters@0: 	if ( isoff( outstr ) ) {
tomwalters@0: 	    if ( i == 0 ) ofp = stdout ;
tomwalters@0: 	    else if ( ( ofp = fopen( argv[argc-i], "w" ) ) == (FILE *)0 ) {
tomwalters@0: 		fprintf( stderr,"scale: can't create %s\n", argv[argc-i] ) ;
tomwalters@0: 		exit( 1 ) ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0: 	else if ( isstr( outstr, "stdout") ) ofp = stdout ;
tomwalters@0: 	else if ( ( ofp = fopen( outstr, "w" ) ) == (FILE *)0 ) {
tomwalters@0: 	    fprintf( stderr,"scale: can't create %s\n", outstr ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 	first = 0 ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     else {
tomwalters@0: 	if ( isoff( outstr ) ) {
tomwalters@0: 	    fclose ( ofp ) ;
tomwalters@0: 	    if ( i == 0 ) ofp = stdout ;
tomwalters@0: 	    else if ( ( ofp = fopen( argv[argc-i], "w" ) ) == (FILE *)0 ) {
tomwalters@0: 		fprintf( stderr,"scale: can't create %s\n", argv[argc-i] ) ;
tomwalters@0: 		exit( 1 ) ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Return parameters `shift' and `scale' to normalize the data for zero mean
tomwalters@0: and unit std dev:
tomwalters@0: 	shift = mean
tomwalters@0: 	scale = 1 / sqrt(variance)
tomwalters@0: */
tomwalters@0: 
tomwalters@0: getnorm( data, n, shift, scale )
tomwalters@0: float *data ;
tomwalters@0: int    n    ;
tomwalters@0: float *shift, *scale  ;
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   += data[i] ;
tomwalters@0: 	sumsq += data[i] * data[i] ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     *shift = sum / n ;
tomwalters@0:     *scale = 1. / sqrt( sumsq / n - *shift * *shift ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Return parameters `shift' and `scale' to fit the given float data into the
tomwalters@0: amplitude range delimited by `top' and `bottom'.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: getrange( data, n, top, bottom, shift, scale )
tomwalters@0: float *data ;
tomwalters@0: int    n    ;
tomwalters@0: float  top,    bottom ;
tomwalters@0: float *shift, *scale  ;
tomwalters@0: {
tomwalters@0:     float  min = 1e+8, max = (-1e+8) ;
tomwalters@0:     int    i ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 	if ( data[i] > max ) max = data[i] ;
tomwalters@0: 	if ( data[i] < min ) min = data[i] ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     *scale = ( top - bottom ) / ( max - min ) ;
tomwalters@0:     *shift = ( min * top - max * bottom ) / ( top - bottom ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: