tomwalters@0: /*
tomwalters@0:   Autocorrelation function using FFT routines adapted from routines given in
tomwalters@0:   Numerical Recipes.
tomwalters@0:   The routine acf() uses an fft to compute the array of lagged products
tomwalters@0:   which is the autocovariance function. An optional normalization (dividing
tomwalters@0:   each coefficient by the zeroth coefficient) produces the autocorrelation
tomwalters@0:   function, (and this is then scaled up to the range [0-1000]).
tomwalters@0: 
tomwalters@0:   The output framewidth is the given maximum acf lag, which must not be
tomwalters@0:   greater than the input framewidth.
tomwalters@0: 
tomwalters@0:   Padding with zeroes is recommended to counter end-effects (See NR, p416).
tomwalters@0:   It is also necessary since:
tomwalters@0: 1) framesize + padding must be a power of two (sample points at given rate).
tomwalters@0: 2) max possible acf lag = ( framesize + padding ) / 2
tomwalters@0: 
tomwalters@0:   Each input frame is padded with zeroes so that it is at least twice the
tomwalters@0:   required max acf lag, and is also a power of 2.
tomwalters@0:   This is because the fft method is "radix-2", and because the output acf is
tomwalters@0:   symmetrical so that just half the points are unique.
tomwalters@0:   Therefore each input frame is padded with zeroes to the next power of 2
tomwalters@0:   larger than either the input framewidth, or twice the required max acf lag,
tomwalters@0:   (whichever is  the larger).
tomwalters@0: 
tomwalters@0:   If necessary, extra padding can be enforced using the padding option to
tomwalters@0:   add extra zeroes, padding to a larger power of 2.
tomwalters@0:   The amount of extra padding is "exponential", expanding the basic size to:
tomwalters@0: 	( framesize + padding ) * 2**n
tomwalters@0:   where the padding option is n.
tomwalters@0:   (n=0 by default, so that no extra padding is added. When n=1 then padding is
tomwalters@0:   added to double the size, and when n=2 the size is quadrupled, etc.).
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:   Other options (set using <option>=on):
tomwalters@0:     size         Print the input and output frame sizes in sample points.
tomwalters@0:     echo         Echo buffered input without processing.
tomwalters@0:     negative     Write negative half of acf (ie. with zeroth lag on right).
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:    within a frame of 12.8ms, plotting the 2nd frame in a sequence of frames
tomwalters@0:    with frameshift 12.8ms.
tomwalters@0: 
tomwalters@0: acf samp=10kHz width=12.8ms 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:    each computed within a frame of 12.8ms, and shifted by 2 sample points.
tomwalters@0: 
tomwalters@0: acf samp=10kHz width=12.8ms 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, *widthstr, *sizestr ;
tomwalters@0: static char *padstr, *wstr, *lagstr, *scalestr, *normstr, *framestr, *echostr, *negstr ;
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"    ,   "16ms"      ,  &shiftstr    ,   "Step between input frames."            , VAL     },
tomwalters@0:     {   "width"     ,   "32ms"      ,  &widthstr    ,   "Width of input frames."                , VAL     },
tomwalters@0:     {   "lag"       ,   "32ms"      ,  &lagstr      ,   "Maximum acf lag (output frame width)." , VAL     },
tomwalters@0:     {   "padding"   ,   "0"         ,  &padstr      ,   "Exponential frame padding"             , SVAL    },
tomwalters@0:     {   "window"    ,   "on"        ,  &wstr        ,   "Hamming window"                        , SETFLAG },
tomwalters@0:     {   "normalize" ,   "off"       ,  &normstr     ,   "normalized acf"                        , SETFLAG },
tomwalters@0:     {   "scale"     ,   "0.025"     ,  &scalestr    ,   "scale output"                          , 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     width       ;
tomwalters@0: int     step        ;
tomwalters@0: int     zeroes      ;
tomwalters@0: int     maxlag      ;
tomwalters@0: int     dowindow    ;
tomwalters@0: double  scale       ;
tomwalters@0: int     normalize   ;
tomwalters@0: int     echo        ;
tomwalters@0: int     negative    ;
tomwalters@0: 
tomwalters@0: float  *vec, *window ;
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:     short    *buf ;
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:     dowindow   = ison(  wstr    ) ;
tomwalters@0:     normalize  = ison(  normstr ) ;
tomwalters@0:     echo       = ison(  echostr ) ;
tomwalters@0:     negative   = ison(  negstr  ) ;
tomwalters@0: 
tomwalters@0:     width      = to_p( widthstr, samplerate )  ;
tomwalters@0:     step       = to_p( shiftstr, samplerate )  ;
tomwalters@0:     maxlag     = to_p( lagstr  , samplerate )  ;
tomwalters@0: 
tomwalters@0:     if ( ( maxlag << 1 ) < width )
tomwalters@0: 	zeroes = ( getpower( width ) << atoi( padstr ) ) - width ;
tomwalters@0:     else if ( maxlag <= width )
tomwalters@0: 	zeroes = ( getpower( maxlag << 1 ) << atoi( padstr ) ) - width ;
tomwalters@0:     else {
tomwalters@0: 	fprintf(stderr,"acf: required max lag [%dms] greater than framewidth [%dms]\n", maxlag*1000/samplerate, width*1000/samplerate);
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     scale = ( 1./(double)(width+zeroes) ) * atof( scalestr ) ;
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,"    input  frame size = %d  (framewidth=%d + padding=%d)\n", width+zeroes, width, zeroes ) ;
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( width+zeroes * sizeof(float) ) ;
tomwalters@0:     if ( dowindow )
tomwalters@0: 	window = hamming( width ) ;
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:     while ( ( buf = getframe( fp, a, width, step ) ) != (short *)0  && ( a<=b || b==0 ) ) {
tomwalters@0: 	if ( echo ) fwrite( buf, sizeof(short), width, stdout ) ;
tomwalters@0: 	else process( buf ) ;
tomwalters@0: 	a++ ;
tomwalters@0:     }
tomwalters@0:     if ( a<=b && b>0 )
tomwalters@0: 	fprintf(stderr,"warning: not enough frames for request\n");
tomwalters@0: 
tomwalters@0:     fclose( fp ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: process( buf )
tomwalters@0: short *buf ;
tomwalters@0: {
tomwalters@0:     short *bptr   = buf ;
tomwalters@0:     short *endptr = buf + width ;
tomwalters@0:     short *optr   = obuf ;
tomwalters@0:     short *endlag = obuf + maxlag ;
tomwalters@0:     float *wptr   = window ;
tomwalters@0:     float *vptr   = vec    ;
tomwalters@0:     float *endvec = vec + width + zeroes ;
tomwalters@0: 
tomwalters@0:     if ( dowindow )
tomwalters@0: 	while ( bptr < endptr )
tomwalters@0: 	    *vptr++ =  *bptr++  *  *wptr++ ;
tomwalters@0:     else
tomwalters@0: 	while ( bptr < endptr )
tomwalters@0: 	    *vptr++ =  *bptr++ ;
tomwalters@0: 
tomwalters@0:     while ( vptr < endvec )         /* padding */
tomwalters@0: 	*vptr++ =  0 ;
tomwalters@0: 
tomwalters@0:     acf( vec, width+zeroes ) ;
tomwalters@0: 
tomwalters@0:     if ( normalize )
tomwalters@0: 	scale = 1000. / *vec ;
tomwalters@0: 
tomwalters@0:     vptr = vec ;
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-- ) /* write in reverse */
tomwalters@0: 	    fwrite( optr, sizeof(short), 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: