tomwalters@0: /*
tomwalters@0:   Autocorrelation function using recursive "leaky integrator" filter
tomwalters@0:   ie. for the m'th lag, the n'th recursive update is:
tomwalters@0: 
tomwalters@0: 	y[n] = a.y[n-1] + (1-a).f[n].f[n-m]
tomwalters@0: 
tomwalters@0:   where decay constant a = exp(-Ts/T)
tomwalters@0:   and   Ts is sample interval in seconds
tomwalters@0: 	T  is decay time constant in seconds
tomwalters@0: 
tomwalters@0:   The recursion computes the mean value of the product f[n].f[n-m]
tomwalters@0:   within an exponential window which weights past values.
tomwalters@0:   The window width parameter (eg. corresponding with the width of the
tomwalters@0:   fft analysis window) is the time constant parameter.
tomwalters@0:   The half life of the exponential window is given by -T.ln(0.5) = 0.693T
tomwalters@0:   ie. for a given time constant, T secs, the window decays to half of its
tomwalters@0:   current value in a period of 0.693T secs.
tomwalters@0:   So T can be set to accomodate an expected period.
tomwalters@0:   As a rough guide, in a period of 0.5T the window decays to about 0.6
tomwalters@0:   (ie by 60%), so set T to twice the expected period of the waveform to
tomwalters@0:   get 60% of the window over this period.
tomwalters@0:   Eg. if the waveform has an expected period of t ms, then to set 60% of
tomwalters@0:   the window to this period, set 0.5T = t ms, ie. T = 2t ms.
tomwalters@0: 
tomwalters@0:   The routine generates the autocovariance function.
tomwalters@0:   An optional normalization (dividing each coefficient by the zeroth
tomwalters@0:   coefficient) produces the autocorrelation function,
tomwalters@0:   (and this is then scaled up to the range [0-1000]).
tomwalters@0: 
tomwalters@0:   The output framewidth is the given maximum acf lag.
tomwalters@0: 
tomwalters@0:   Frames are selected from the input stream using the "frames" option,
tomwalters@0:   which has syntax: [[-]frame=a[-b]]. Input frames  are numbered 1,2,...
tomwalters@0:   For example:
tomwalters@0:     frame=a      Select just the a'th frame.
tomwalters@0:     frame=a-b    Select frames from the a'th to b'th inclusive.
tomwalters@0:   The strings "min" and "max" can be used as specifiers, meaning eg:
tomwalters@0:     frame=min    Select the first frame.
tomwalters@0:     frame=max    Select the last frame.
tomwalters@0:     frame=a-max  Select frames from the a'th to the last inclusive.
tomwalters@0:     frame=min-b  Select frames from the first to the b'th inclusive.
tomwalters@0:   The default selects all frames, (ie frame=min-max).
tomwalters@0: 
tomwalters@0: 
tomwalters@0:   Examples:
tomwalters@0: 
tomwalters@0: 1. To print the input and output frame sizes in sample points, eg for a
tomwalters@0:    subsequent plotting program, use the size option:
tomwalters@0: 
tomwalters@0: acf ... size=on
tomwalters@0: 
tomwalters@0: 2. An acf of a waveform sampled at 10kHz, computed to a max lag of 12.8ms,
tomwalters@0:    plotting the 2nd frame in a sequence of frames with frameshift 12.8ms.
tomwalters@0: 
tomwalters@0: acf samp=10kHz frstep=12.8ms frame=2 lag=12.8ms file | x11plot
tomwalters@0: 
tomwalters@0: 3. An animated plot of successive acf's of a waveform sampled at 10kHz,
tomwalters@0:    shifted by 2 sample points.
tomwalters@0: 
tomwalters@0: acf samp=10kHz frstep=2p lag=12.8ms file | x11play -n128
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: #include "sigproc.h"
tomwalters@0: 
tomwalters@0: char applic[]     = "autocorrelation function of contiguous frames.\n      i/p and o/p data in binary shorts." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *sampstr,  *startstr, *shiftstr, *sizestr ;
tomwalters@0: static char *lagstr,  *tcstr,    *scalestr, *normstr,  *framestr, *echostr ;
tomwalters@0: static char *negstr,  *upstr,    *slopestr ;
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:     {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p file."              , VAL     },
tomwalters@0:     {   "frames"    ,   "1-max"     ,  &framestr    ,   "Select frames inclusively"             , VAL     },
tomwalters@0:     {   "frstep"    ,   "1ms"       ,  &shiftstr    ,   "Step between output frames."           , VAL     },
tomwalters@0:     {   "lag"       ,   "32ms"      ,  &lagstr      ,   "Maximum acf lag (output frame width)." , VAL     },
tomwalters@0:     {   "update"    ,   "1p"        ,  &upstr       ,   "Step between recursive updates."       , VAL     },
tomwalters@0:     {   "timeconst" ,   "3ms"       ,  &tcstr       ,   "Decay time constant"                   , VAL     },
tomwalters@0:     {   "normalize" ,   "off"       ,  &normstr     ,   "Normalized acf"                        , SETFLAG },
tomwalters@0:     {   "scale"     ,   "0.003"     ,  &scalestr    ,   "scale output"                          , SVAL    },
tomwalters@0:     {   "sigmoid"   ,   "off"       ,  &slopestr    ,   "Sigmoidal weighting slope parameter"   , SVAL    },
tomwalters@0:     {   "size"      ,   "off"       ,  &sizestr     ,   "Print input/output frame size in samples"      , SETFLAG },
tomwalters@0:     {   "echo"      ,   "off"       ,  &echostr     ,   "echo buffered input without processing", SVAL    },
tomwalters@0:     {   "negative"  ,   "off"       ,  &negstr      ,   "write negative half of acf (zeroth lag on right)", SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int     samplerate  ;
tomwalters@0: int     framestep   ;
tomwalters@0: int     updatestep  ;
tomwalters@0: int     maxlag      ;
tomwalters@0: double  scale       ;
tomwalters@0: int     normalize   ;
tomwalters@0: int     echo        ;
tomwalters@0: int     negative    ;
tomwalters@0: double  alpha       ;   /* decay constant */
tomwalters@0: 
tomwalters@0: float  *vec   ;
tomwalters@0: short  *obuf  ;
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int     argc;
tomwalters@0: char  **argv;
tomwalters@0: {
tomwalters@0:     FILE     *fp  ;
tomwalters@0:     int       a, b ;
tomwalters@0: 
tomwalters@0:     fp = openopts( 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:     normalize  = ison(  normstr ) ;
tomwalters@0:     echo       = ison(  echostr ) ;
tomwalters@0:     negative   = ison(  negstr  ) ;
tomwalters@0: 
tomwalters@0:     framestep  = to_p( shiftstr, samplerate )  ;
tomwalters@0:     updatestep = to_p( upstr   , samplerate )  ;
tomwalters@0:     maxlag     = to_p( lagstr  , samplerate )  ;
tomwalters@0: 
tomwalters@0:     scale      = atof( scalestr ) ;
tomwalters@0: 
tomwalters@0:     alpha      = exp( -(double)( 1. / ( samplerate * to_s( tcstr, samplerate ) ) ) ) ;
tomwalters@0: 
tomwalters@0:     if ( framestep % updatestep != 0 ) {
tomwalters@0: 	fprintf(stderr,"frstep [%dp] must be an integer multiple of update [%dp] \n", framestep, updatestep ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* frame size printout */
tomwalters@0: 
tomwalters@0:     if ( ison( sizestr ) ) {
tomwalters@0: 	fprintf(stderr,"acf sizes in sample points:\n" ) ;
tomwalters@0: 	fprintf(stderr,"    output frame size = %d \n", maxlag ) ;
tomwalters@0: 	exit( 0 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* parse bounds on number of frames */
tomwalters@0: 
tomwalters@0:     if ( selector( framestr, &a, &b ) == 0 ) {
tomwalters@0: 	fprintf(stderr,"acf: bad frame selector [%s]\n", framestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* Allocate working space */
tomwalters@0: 
tomwalters@0:     obuf  = (short *)malloc( maxlag * sizeof(short) ) ;
tomwalters@0:     vec   = (float *)malloc( maxlag * sizeof(float) ) ;
tomwalters@0: 
tomwalters@0:     /* Compute acf for each frame of width shorts in the input stream */
tomwalters@0: 
tomwalters@0:     if ( seekbytes( fp, (int)( to_p( startstr, samplerate ) * sizeof(short) ) ) == 0 ) {
tomwalters@0: 	fprintf(stderr,"improper seek\n") ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( isoff( slopestr ) )
tomwalters@0: 	process( fp, a, b ) ;                   /* conventional acf recursion */
tomwalters@0:     else
tomwalters@0: 	process2( fp, a, b, atof( slopestr) ) ; /* sigmoidally weighted recursion */
tomwalters@0: 						/* (recommend set sigmoid=.02)    */
tomwalters@0: 
tomwalters@0:     fclose( fp ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: process( fp, a, b )
tomwalters@0: FILE *fp ;
tomwalters@0: int   a, b ;
tomwalters@0: {
tomwalters@0:     register  short  *buf, *bptr, *endptr ;
tomwalters@0:     register  float  *vptr ;
tomwalters@0:     register  double  beta ;
tomwalters@0:     register  int     update_count = 0 ;
tomwalters@0:     register  int     frame_count  = 0 ;
tomwalters@0: 
tomwalters@0:     /* Initialize process with the first maxlag samples */
tomwalters@0: 
tomwalters@0:     if ( ( buf = getframe( fp, 1, maxlag, updatestep ) ) != (short *)0 ) {
tomwalters@0: 	for ( endptr = buf ; endptr < buf + maxlag ; endptr++ ) {
tomwalters@0: 	    bptr = endptr ;
tomwalters@0: 	    vptr = vec ;
tomwalters@0: 	    beta = ( 1-alpha ) * *bptr ;
tomwalters@0: 	    while ( bptr >= buf )
tomwalters@0: 		*vptr++ =  alpha * *vptr + beta * *bptr-- ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	update_count = 1 ;
tomwalters@0: 	frame_count  = 1 ;
tomwalters@0: 
tomwalters@0: 	if ( frame_count >= a ) {
tomwalters@0: 	    if ( echo )
tomwalters@0: 		fwrite( buf, sizeof(short), maxlag, stdout ) ;
tomwalters@0: 	    else
tomwalters@0: 		writeframe() ;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* Continue process with successive buffers of maxlag samples */
tomwalters@0: 
tomwalters@0:     while ( ( buf = getframe( fp, ++update_count, maxlag, updatestep ) ) != (short *)0 && ( frame_count<=b || b==0 ) ) {
tomwalters@0: 	bptr = buf + maxlag - 1 ;
tomwalters@0: 	vptr = vec ;
tomwalters@0: 	beta = ( 1-alpha ) * *bptr ;
tomwalters@0: 
tomwalters@0: 	while ( bptr >= buf )     /* recursive update */
tomwalters@0: 	    *vptr++ =  alpha * *vptr + beta * *bptr-- ;
tomwalters@0: 
tomwalters@0: 	if ( update_count % framestep == 1 ) {
tomwalters@0: 	    if ( ++frame_count >= a && ( frame_count <= b || b == 0 ) ) {
tomwalters@0: 		if ( echo )
tomwalters@0: 		    fwrite( buf, sizeof(short), maxlag, stdout ) ;
tomwalters@0: 		else
tomwalters@0: 		    writeframe() ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0:     if ( frame_count <= b && b > 0 )
tomwalters@0: 	fprintf(stderr,"warning: not enough frames for request\n");
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: writeframe()
tomwalters@0: {
tomwalters@0:     register  float  *vptr   = vec ;
tomwalters@0:     register  short  *optr   = obuf ;
tomwalters@0:     register  short  *endlag = obuf + maxlag ;
tomwalters@0: 
tomwalters@0:     if ( normalize )
tomwalters@0: 	scale = 1000. / *vec ;
tomwalters@0: 
tomwalters@0:     while ( optr < endlag )
tomwalters@0: 	*optr++ = (short)( scale * *vptr++ ) ;
tomwalters@0: 
tomwalters@0:     if ( !negative )
tomwalters@0: 	fwrite( obuf, sizeof(short), maxlag, stdout ) ;
tomwalters@0:     else {
tomwalters@0: 	for ( optr = endlag-1 ; optr >= obuf ; optr-- )
tomwalters@0: 	    fwrite( optr, sizeof(short), 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*************************************************************************/
tomwalters@0: /*
tomwalters@0:   same as `process', but with added sigmoidal weighting
tomwalters@0: */
tomwalters@0: 
tomwalters@0: process2( fp, a, b, slope )
tomwalters@0: FILE *fp ;
tomwalters@0: int   a, b  ;
tomwalters@0: float slope ;
tomwalters@0: {
tomwalters@0:     register  short  *buf, *bptr, *endptr ;
tomwalters@0:     register  float  *vptr ;
tomwalters@0:     register  double  beta ;
tomwalters@0:     register  int     update_count = 0 ;
tomwalters@0:     register  int     frame_count  = 0 ;
tomwalters@0:     double    sigmoid() ;
tomwalters@0: 
tomwalters@0:     /* Initialize process with the first maxlag samples */
tomwalters@0: 
tomwalters@0:     if ( ( buf = getframe( fp, 1, maxlag, updatestep ) ) != (short *)0 ) {
tomwalters@0: 	for ( endptr = buf ; endptr < buf + maxlag ; endptr++ ) {
tomwalters@0: 	    bptr = endptr ;
tomwalters@0: 	    vptr = vec ;
tomwalters@0: 	    beta = ( 1-alpha ) * *bptr ;
tomwalters@0: 	    while ( bptr >= buf )
tomwalters@0: 		*vptr++ =  alpha * *vptr + beta * sigmoid( fabs((double)*endptr) - fabs((double)*bptr ), slope ) * *bptr-- ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	update_count = 1 ;
tomwalters@0: 	frame_count  = 1 ;
tomwalters@0: 
tomwalters@0: 	if ( frame_count >= a ) {
tomwalters@0: 	    if ( echo )
tomwalters@0: 		fwrite( buf, sizeof(short), maxlag, stdout ) ;
tomwalters@0: 	    else
tomwalters@0: 		writeframe() ;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* Continue process with successive buffers of maxlag samples */
tomwalters@0: 
tomwalters@0:     while ( ( buf = getframe( fp, ++update_count, maxlag, updatestep ) ) != (short *)0 && ( frame_count<=b || b==0 ) ) {
tomwalters@0: 	endptr = buf + maxlag - 1 ;
tomwalters@0: 	bptr   = endptr ;
tomwalters@0: 	vptr   = vec ;
tomwalters@0: 	beta   = ( 1-alpha ) * *bptr ;
tomwalters@0: 
tomwalters@0: 	while ( bptr >= buf )     /* recursive update */
tomwalters@0: 	    *vptr++ =  alpha * *vptr + beta * sigmoid( fabs((double)*endptr) - fabs((double)*bptr ), slope ) * *bptr-- ;
tomwalters@0: 
tomwalters@0: 	if ( update_count % framestep == 1 ) {
tomwalters@0: 	    if ( ++frame_count >= a && ( frame_count <= b || b == 0 ) ) {
tomwalters@0: 		if ( echo )
tomwalters@0: 		    fwrite( buf, sizeof(short), maxlag, stdout ) ;
tomwalters@0: 		else
tomwalters@0: 		    writeframe() ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0:     if ( frame_count <= b && b > 0 )
tomwalters@0: 	fprintf(stderr,"warning: not enough frames for request\n");
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: double sigmoid( x, slope )
tomwalters@0: double x ;
tomwalters@0: float  slope ;
tomwalters@0: {
tomwalters@0:     return (double)( 1. / ( exp( -(double)( 4.595 + slope * x ) ) + 1. ) ) ;
tomwalters@0: }