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

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

/*
  Autocorrelation function using FFT routines adapted from routines given in
  Numerical Recipes.
  The routine acf() uses an fft to compute the array of lagged products
  which is 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, which must not be
  greater than the input framewidth.

  Padding with zeroes is recommended to counter end-effects (See NR, p416).
  It is also necessary since:
1) framesize + padding must be a power of two (sample points at given rate).
2) max possible acf lag = ( framesize + padding ) / 2

  Each input frame is padded with zeroes so that it is at least twice the
  required max acf lag, and is also a power of 2.
  This is because the fft method is "radix-2", and because the output acf is
  symmetrical so that just half the points are unique.
  Therefore each input frame is padded with zeroes to the next power of 2
  larger than either the input framewidth, or twice the required max acf lag,
  (whichever is  the larger).

  If necessary, extra padding can be enforced using the padding option to
  add extra zeroes, padding to a larger power of 2.
  The amount of extra padding is "exponential", expanding the basic size to:
	( framesize + padding ) * 2**n
  where the padding option is n.
  (n=0 by default, so that no extra padding is added. When n=1 then padding is
  added to double the size, and when n=2 the size is quadrupled, etc.).

  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).

  Other options (set using <option>=on):
    size         Print the input and output frame sizes in sample points.
    echo         Echo buffered input without processing.
    negative     Write negative half of acf (ie. with zeroth lag on right).


  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
   within a frame of 12.8ms, plotting the 2nd frame in a sequence of frames
   with frameshift 12.8ms.

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

3. An animated plot of successive acf's of a waveform sampled at 10kHz,
   each computed within a frame of 12.8ms, and shifted by 2 sample points.

acf samp=10kHz width=12.8ms 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, *widthstr, *sizestr ;
static char *padstr, *wstr, *lagstr, *scalestr, *normstr, *framestr, *echostr, *negstr ;

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"    ,   "16ms"      ,  &shiftstr    ,   "Step between input frames."            , VAL     },
    {   "width"     ,   "32ms"      ,  &widthstr    ,   "Width of input frames."                , VAL     },
    {   "lag"       ,   "32ms"      ,  &lagstr      ,   "Maximum acf lag (output frame width)." , VAL     },
    {   "padding"   ,   "0"         ,  &padstr      ,   "Exponential frame padding"             , SVAL    },
    {   "window"    ,   "on"        ,  &wstr        ,   "Hamming window"                        , SETFLAG },
    {   "normalize" ,   "off"       ,  &normstr     ,   "normalized acf"                        , SETFLAG },
    {   "scale"     ,   "0.025"     ,  &scalestr    ,   "scale output"                          , 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     width       ;
int     step        ;
int     zeroes      ;
int     maxlag      ;
int     dowindow    ;
double  scale       ;
int     normalize   ;
int     echo        ;
int     negative    ;

float  *vec, *window ;
short  *obuf ;

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

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

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

    width      = to_p( widthstr, samplerate )  ;
    step       = to_p( shiftstr, samplerate )  ;
    maxlag     = to_p( lagstr  , samplerate )  ;

    if ( ( maxlag << 1 ) < width )
	zeroes = ( getpower( width ) << atoi( padstr ) ) - width ;
    else if ( maxlag <= width )
	zeroes = ( getpower( maxlag << 1 ) << atoi( padstr ) ) - width ;
    else {
	fprintf(stderr,"acf: required max lag [%dms] greater than framewidth [%dms]\n", maxlag*1000/samplerate, width*1000/samplerate);
	exit( 1 ) ;
    }

    scale = ( 1./(double)(width+zeroes) ) * atof( scalestr ) ;

    /* frame size printout */

    if ( ison( sizestr ) ) {
	fprintf(stderr,"acf sizes in sample points:\n" ) ;
	fprintf(stderr,"    input  frame size = %d  (framewidth=%d + padding=%d)\n", width+zeroes, width, zeroes ) ;
	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( width+zeroes * sizeof(float) ) ;
    if ( dowindow )
	window = hamming( width ) ;

    /* 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 ) ;
    }

    while ( ( buf = getframe( fp, a, width, step ) ) != (short *)0  && ( a<=b || b==0 ) ) {
	if ( echo ) fwrite( buf, sizeof(short), width, stdout ) ;
	else process( buf ) ;
	a++ ;
    }
    if ( a<=b && b>0 )
	fprintf(stderr,"warning: not enough frames for request\n");

    fclose( fp ) ;
}


process( buf )
short *buf ;
{
    short *bptr   = buf ;
    short *endptr = buf + width ;
    short *optr   = obuf ;
    short *endlag = obuf + maxlag ;
    float *wptr   = window ;
    float *vptr   = vec    ;
    float *endvec = vec + width + zeroes ;

    if ( dowindow )
	while ( bptr < endptr )
	    *vptr++ =  *bptr++  *  *wptr++ ;
    else
	while ( bptr < endptr )
	    *vptr++ =  *bptr++ ;

    while ( vptr < endvec )         /* padding */
	*vptr++ =  0 ;

    acf( vec, width+zeroes ) ;

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

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

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