view tools/racf.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

/*
  Autocorrelation function using recursive "leaky integrator" filter
  ie. for the m'th lag, the n'th recursive update is:

	y[n] = a.y[n-1] + (1-a).f[n].f[n-m]

  where decay constant a = exp(-Ts/T)
  and   Ts is sample interval in seconds
	T  is decay time constant in seconds

  The recursion computes the mean value of the product f[n].f[n-m]
  within an exponential window which weights past values.
  The window width parameter (eg. corresponding with the width of the
  fft analysis window) is the time constant parameter.
  The half life of the exponential window is given by -T.ln(0.5) = 0.693T
  ie. for a given time constant, T secs, the window decays to half of its
  current value in a period of 0.693T secs.
  So T can be set to accomodate an expected period.
  As a rough guide, in a period of 0.5T the window decays to about 0.6
  (ie by 60%), so set T to twice the expected period of the waveform to
  get 60% of the window over this period.
  Eg. if the waveform has an expected period of t ms, then to set 60% of
  the window to this period, set 0.5T = t ms, ie. T = 2t ms.

  The routine generates the autocovariance function.
  An optional normalization (dividing each coefficient by the zeroth
  coefficient) produces the autocorrelation function,
  (and this is then scaled up to the range [0-1000]).

  The output framewidth is the given maximum acf lag.

  Frames are selected from the input stream using the "frames" option,
  which has syntax: [[-]frame=a[-b]]. Input frames  are numbered 1,2,...
  For example:
    frame=a      Select just the a'th frame.
    frame=a-b    Select frames from the a'th to b'th inclusive.
  The strings "min" and "max" can be used as specifiers, meaning eg:
    frame=min    Select the first frame.
    frame=max    Select the last frame.
    frame=a-max  Select frames from the a'th to the last inclusive.
    frame=min-b  Select frames from the first to the b'th inclusive.
  The default selects all frames, (ie frame=min-max).


  Examples:

1. To print the input and output frame sizes in sample points, eg for a
   subsequent plotting program, use the size option:

acf ... size=on

2. An acf of a waveform sampled at 10kHz, computed to a max lag of 12.8ms,
   plotting the 2nd frame in a sequence of frames with frameshift 12.8ms.

acf samp=10kHz frstep=12.8ms frame=2 lag=12.8ms file | x11plot

3. An animated plot of successive acf's of a waveform sampled at 10kHz,
   shifted by 2 sample points.

acf samp=10kHz frstep=2p lag=12.8ms file | x11play -n128

*/

#include <stdio.h>
#include <math.h>
#include "options.h"
#include "units.h"
#include "strmatch.h"
#include "sigproc.h"

char applic[]     = "autocorrelation function of contiguous frames.\n      i/p and o/p data in binary shorts." ;

static char *helpstr, *debugstr, *sampstr,  *startstr, *shiftstr, *sizestr ;
static char *lagstr,  *tcstr,    *scalestr, *normstr,  *framestr, *echostr ;
static char *negstr,  *upstr,    *slopestr ;

static Options option[] = {
    {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                  , DEBUG   },
    {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                      , DEBUG   },
    {   "samplerate",   "20kHz"     ,  &sampstr     ,   "Samplerate "                           , VAL     },
    {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p file."              , VAL     },
    {   "frames"    ,   "1-max"     ,  &framestr    ,   "Select frames inclusively"             , VAL     },
    {   "frstep"    ,   "1ms"       ,  &shiftstr    ,   "Step between output frames."           , VAL     },
    {   "lag"       ,   "32ms"      ,  &lagstr      ,   "Maximum acf lag (output frame width)." , VAL     },
    {   "update"    ,   "1p"        ,  &upstr       ,   "Step between recursive updates."       , VAL     },
    {   "timeconst" ,   "3ms"       ,  &tcstr       ,   "Decay time constant"                   , VAL     },
    {   "normalize" ,   "off"       ,  &normstr     ,   "Normalized acf"                        , SETFLAG },
    {   "scale"     ,   "0.003"     ,  &scalestr    ,   "scale output"                          , SVAL    },
    {   "sigmoid"   ,   "off"       ,  &slopestr    ,   "Sigmoidal weighting slope parameter"   , SVAL    },
    {   "size"      ,   "off"       ,  &sizestr     ,   "Print input/output frame size in samples"      , SETFLAG },
    {   "echo"      ,   "off"       ,  &echostr     ,   "echo buffered input without processing", SVAL    },
    {   "negative"  ,   "off"       ,  &negstr      ,   "write negative half of acf (zeroth lag on right)", SVAL    },
   ( char * ) 0 } ;


int     samplerate  ;
int     framestep   ;
int     updatestep  ;
int     maxlag      ;
double  scale       ;
int     normalize   ;
int     echo        ;
int     negative    ;
double  alpha       ;   /* decay constant */

float  *vec   ;
short  *obuf  ;

main (argc, argv)
int     argc;
char  **argv;
{
    FILE     *fp  ;
    int       a, b ;

    fp = openopts( option,argc,argv ) ;
    if ( !isoff( helpstr ) )
	helpopts( helpstr, argv[0], applic, option ) ;

    samplerate = to_Hz( sampstr ) ;
    normalize  = ison(  normstr ) ;
    echo       = ison(  echostr ) ;
    negative   = ison(  negstr  ) ;

    framestep  = to_p( shiftstr, samplerate )  ;
    updatestep = to_p( upstr   , samplerate )  ;
    maxlag     = to_p( lagstr  , samplerate )  ;

    scale      = atof( scalestr ) ;

    alpha      = exp( -(double)( 1. / ( samplerate * to_s( tcstr, samplerate ) ) ) ) ;

    if ( framestep % updatestep != 0 ) {
	fprintf(stderr,"frstep [%dp] must be an integer multiple of update [%dp] \n", framestep, updatestep ) ;
	exit( 1 ) ;
    }

    /* frame size printout */

    if ( ison( sizestr ) ) {
	fprintf(stderr,"acf sizes in sample points:\n" ) ;
	fprintf(stderr,"    output frame size = %d \n", maxlag ) ;
	exit( 0 ) ;
    }

    /* parse bounds on number of frames */

    if ( selector( framestr, &a, &b ) == 0 ) {
	fprintf(stderr,"acf: bad frame selector [%s]\n", framestr ) ;
	exit( 1 ) ;
    }

    /* Allocate working space */

    obuf  = (short *)malloc( maxlag * sizeof(short) ) ;
    vec   = (float *)malloc( maxlag * sizeof(float) ) ;

    /* Compute acf for each frame of width shorts in the input stream */

    if ( seekbytes( fp, (int)( to_p( startstr, samplerate ) * sizeof(short) ) ) == 0 ) {
	fprintf(stderr,"improper seek\n") ;
	exit( 1 ) ;
    }

    if ( isoff( slopestr ) )
	process( fp, a, b ) ;                   /* conventional acf recursion */
    else
	process2( fp, a, b, atof( slopestr) ) ; /* sigmoidally weighted recursion */
						/* (recommend set sigmoid=.02)    */

    fclose( fp ) ;
}


process( fp, a, b )
FILE *fp ;
int   a, b ;
{
    register  short  *buf, *bptr, *endptr ;
    register  float  *vptr ;
    register  double  beta ;
    register  int     update_count = 0 ;
    register  int     frame_count  = 0 ;

    /* Initialize process with the first maxlag samples */

    if ( ( buf = getframe( fp, 1, maxlag, updatestep ) ) != (short *)0 ) {
	for ( endptr = buf ; endptr < buf + maxlag ; endptr++ ) {
	    bptr = endptr ;
	    vptr = vec ;
	    beta = ( 1-alpha ) * *bptr ;
	    while ( bptr >= buf )
		*vptr++ =  alpha * *vptr + beta * *bptr-- ;
	}

	update_count = 1 ;
	frame_count  = 1 ;

	if ( frame_count >= a ) {
	    if ( echo )
		fwrite( buf, sizeof(short), maxlag, stdout ) ;
	    else
		writeframe() ;
	}
    }

    /* Continue process with successive buffers of maxlag samples */

    while ( ( buf = getframe( fp, ++update_count, maxlag, updatestep ) ) != (short *)0 && ( frame_count<=b || b==0 ) ) {
	bptr = buf + maxlag - 1 ;
	vptr = vec ;
	beta = ( 1-alpha ) * *bptr ;

	while ( bptr >= buf )     /* recursive update */
	    *vptr++ =  alpha * *vptr + beta * *bptr-- ;

	if ( update_count % framestep == 1 ) {
	    if ( ++frame_count >= a && ( frame_count <= b || b == 0 ) ) {
		if ( echo )
		    fwrite( buf, sizeof(short), maxlag, stdout ) ;
		else
		    writeframe() ;
	    }
	}
    }
    if ( frame_count <= b && b > 0 )
	fprintf(stderr,"warning: not enough frames for request\n");
}


writeframe()
{
    register  float  *vptr   = vec ;
    register  short  *optr   = obuf ;
    register  short  *endlag = obuf + maxlag ;

    if ( normalize )
	scale = 1000. / *vec ;

    while ( optr < endlag )
	*optr++ = (short)( scale * *vptr++ ) ;

    if ( !negative )
	fwrite( obuf, sizeof(short), maxlag, stdout ) ;
    else {
	for ( optr = endlag-1 ; optr >= obuf ; optr-- )
	    fwrite( optr, sizeof(short), 1, stdout ) ;
    }
}


/*************************************************************************/
/*
  same as `process', but with added sigmoidal weighting
*/

process2( fp, a, b, slope )
FILE *fp ;
int   a, b  ;
float slope ;
{
    register  short  *buf, *bptr, *endptr ;
    register  float  *vptr ;
    register  double  beta ;
    register  int     update_count = 0 ;
    register  int     frame_count  = 0 ;
    double    sigmoid() ;

    /* Initialize process with the first maxlag samples */

    if ( ( buf = getframe( fp, 1, maxlag, updatestep ) ) != (short *)0 ) {
	for ( endptr = buf ; endptr < buf + maxlag ; endptr++ ) {
	    bptr = endptr ;
	    vptr = vec ;
	    beta = ( 1-alpha ) * *bptr ;
	    while ( bptr >= buf )
		*vptr++ =  alpha * *vptr + beta * sigmoid( fabs((double)*endptr) - fabs((double)*bptr ), slope ) * *bptr-- ;
	}

	update_count = 1 ;
	frame_count  = 1 ;

	if ( frame_count >= a ) {
	    if ( echo )
		fwrite( buf, sizeof(short), maxlag, stdout ) ;
	    else
		writeframe() ;
	}
    }

    /* Continue process with successive buffers of maxlag samples */

    while ( ( buf = getframe( fp, ++update_count, maxlag, updatestep ) ) != (short *)0 && ( frame_count<=b || b==0 ) ) {
	endptr = buf + maxlag - 1 ;
	bptr   = endptr ;
	vptr   = vec ;
	beta   = ( 1-alpha ) * *bptr ;

	while ( bptr >= buf )     /* recursive update */
	    *vptr++ =  alpha * *vptr + beta * sigmoid( fabs((double)*endptr) - fabs((double)*bptr ), slope ) * *bptr-- ;

	if ( update_count % framestep == 1 ) {
	    if ( ++frame_count >= a && ( frame_count <= b || b == 0 ) ) {
		if ( echo )
		    fwrite( buf, sizeof(short), maxlag, stdout ) ;
		else
		    writeframe() ;
	    }
	}
    }
    if ( frame_count <= b && b > 0 )
	fprintf(stderr,"warning: not enough frames for request\n");
}


double sigmoid( x, slope )
double x ;
float  slope ;
{
    return (double)( 1. / ( exp( -(double)( 4.595 + slope * x ) ) + 1. ) ) ;
}