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

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

/*
 acgram.c       Autocorrelogram program by Michael Allerhand.
----------

  Short-time autocorrelation function applied to each row of input frames.
  Output frames consist of row-wise autocorrelation coefficients.
  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-5000]).

Input/output:
  Read a header and frame in NAP format.
  Write a header, (constructed from the original input header), and a
  succession of frames in SAI format.

  The input file is interpreted as one large frame in NAP format to be
  divided into input frames, each to be processed and output in SAI format.

  The options "start" and "length" specify the input file.
  The special option value length=max specifies all the input file from
  the given start to its end.

  The options "width" and "frstep"  specify the input frames.
  The width option is the framewidth of each input frame) and the frstep
  option is the frameshift between input frames in the input file.
  The special option value width=max specifies the input framewidth as equal
  to the given input file length, (and if this is also "max", then the
  input framewidth is the remainder of the input file).

  Most options in the input header are copied to the output header.
  This enables options which are needed for the eventual display
  to pass straight through. Some options are set so that they can override
  the input header. For example, the display option is set on to enable
  display even when input has display=off. The animate option can be set on
  even when the input has animate=off.
  Parts of the header are changed for the new sai format:
  (frames, frameshift, framewidth, frameheight, framebytes).

Padding:
  The output framewidth is the given maximum acf lag, which is also the
  number of acf coefficients computed from each row of each input frame.
  The rows of each input frame are padded with zeroes so that the input
  framewidth is at least twice the required output framewidth (max acf lag),
  and 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. It is advisable to pad input frames for the additional reason that
  padding counters end-effects (which can be significant with the fft method
  for computing acf).

  Therefore each row of each input frame is padded with zeroes to the next
  power of 2 larger than either the original input framewidth or twice the
  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" since if the basic padded
  framesize is N, and the padding option is n, then extra padding is to:
  N * 2**n.
  (n=0 by default, so that N is not changed. When n=1 then N is doubled, and
  when n=2 then N is quadrupled, etc.).


Examples:

1. Autocorrelogram of a NAP, animated and normalized, with max lag 16ms:

gennap len=128ms output=stdout display=off file | acgram lag=16ms norm=on anim=on > file.sai
gensai useprev=on headr=5 top=1000 file             -(for landscape plot)
genspl useprev=on headr=5 top=1000 pensize=2 file   -(for spiral plot)

2. Autocorrelogram of an SAI:
  (Note that gensai removes file.sai, so you must use some other name, eg foo.sai).

gensai len=64 pwidth=64 nwidth=0 output=stdout display=off file |  \
saitonap frame=3 | acgram lag=32ms frame=1 > foo.sai
gensai useprev=on top=1000 headr=5 mag=2 foo

3. Autocorrelogram of a NAP, followed by a cross-channel summary:

gennap len=128ms output=stdout display=off file | acgram lag=16ms scale=0.001 > file.sai
edframe frame=1 file.sai | integframe var=freq > file.sum

*/


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

char applic[]     = "Autocorrelogram auditory image." ;

static char *helpstr  , *debugstr, *dispstr , *anistr   , *startstr  ;
static char *lengthstr, *widthstr, *shiftstr, *headerstr, *framestr  ;
static char *normstr  , *wstr    , *scalestr, *lagstr   , *padstr    ;
static char *negstr   , *vstr ;

static Options option[] = {
    {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                       , DEBUG   },
    {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                           , DEBUG   },
    {   "display"   ,   "on"        ,  &dispstr     ,   "Display output image."                      , SETFLAG },
    {   "animate"   ,   "off"       ,  &anistr      ,   "Animate cartoon."                           , SETFLAG },
    {   "frames"    ,   "1-max"     ,  &framestr    ,   "Select frames inclusively"                  , VAL     },
    {   "Header"    ,   "on"        ,  &headerstr   ,   "Header (for gensai useprevious=on)."        , VAL     },
    {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p file."                   , VAL     },
    {   "length"    ,   "max"       ,  &lengthstr   ,   "Length of i/p file to process."             , 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     ,   "Normalization"                              , SETFLAG },
    {   "scale"     ,   "1."        ,  &scalestr    ,   "Scale factor for output"                    , VAL     },
    {   "negative"  ,   "off"       ,  &negstr      ,   "write negative half of acf (zeroth lag on right)", SVAL    },
    {   "verbose"   ,   "off"       ,  &vstr        ,   "print at each frame"                        , SVAL    },
   ( char * ) 0 } ;


int     samplerate ;

int     frameheight, framewidth ;
int     frames, framebytes ;
int     frameshift ;

int     Newframeheight, Newframewidth ;
int     Newframes, Newframebytes ;
int     Newframeshift ;

int     start, length ; /* Size of input file (num cols) */
int     frameslength ;  /* Length (num cols) of input file to process */

int     zeroes ;
int     window ;        /* Flag for Hamming window */
float   scale  ;        /* Scale factor for output */
int     normalize ;     /* Flag to normalize acf   */
int     negative  ;
int     Verbose   ;

short  *file ;          /* Input file (NAP format) */
float  *W    ;          /* Hamming window for fft  */
float  *buf  ;          /* fft working space       */


main(argc, argv)
int   argc ;
char *argv[] ;
{
    FILE   *fp ;
    char   *header, *SaiHeader();
    char   *versionstr, c ;
    int     startbytes, framebyteshift ;
    short  *frame, *endframe ;
    int     i, a, b ;
    float   acgram(), f, newscale = 1. ;

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

    window    = ison( wstr    ) ;
    normalize = ison( normstr ) ;
    negative  = ison( negstr  ) ;
    Verbose   = ison( vstr    ) ;


    /**** parse bounds on number of frames ****/

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


    /**** Input frame dimensions (from header and options) ****/

    if ((header = ReadHeader(fp)) == (char *) 0 ) {
	fprintf(stderr,"acgram: header not found\n");
	exit(1);
    }

    if ( (versionstr = HeaderString( header, "Version" )) == (char *)0 ) {
	fprintf(stderr,"acgram: model version-number not found in header\n");
	exit(1);
    }

    samplerate  = HeaderSamplerate( header );

    /* In NAP format, the size of the 2-D image is given by:  */
    /*   frameheight = number of rows (filterbank channels)   */
    /*   frames      = number of columns (sample points)      */

    frameheight = HeaderInt( header, "frameheight" );
    frames      = HeaderInt( header, "frames"      );

    if ( ( frameshift  = to_p( shiftstr, samplerate ) ) <= 0 ) {
	fprintf(stderr,"acgram: non-positive frstep [%d]\n", frameshift);
	exit(1);
    }

    /* calculate start of first frame (in cols) */

    if ( ( start = to_p( startstr, samplerate ) + ( a-1 ) * frameshift ) >= frames ) {
	fprintf(stderr,"acgram: input file too small (%dms) for given framing parameters\n", frames * 1000 / samplerate );
	exit(1);
    }

    /* get length of input file (in cols) and framewidth of input frames */

    if ( ismax( lengthstr ) ) length = frames - start ;
    else                      length = to_p( lengthstr, samplerate ) - start ;

    if ( ismax( widthstr ) )  framewidth = length ;
    else                      framewidth = to_p( widthstr, samplerate ) ;

    if ( length < framewidth ) {
	fprintf(stderr,"acgram: input file too small (%dms) for given framing parameters\n", frames * 1000 / samplerate );
	exit(1);
    }

    /* calculate length (num cols) to process, and the number of frames */

    if ( b==0 ) {
	frameslength = length ;
	Newframes = 1 + ( length - framewidth ) / frameshift ;
    }
    else {
	frameslength = framewidth + ( b-a ) * frameshift ;
	Newframes = b - ( a-1 ) ;
    }

    if ( start + frameslength > frames ) {
	fprintf(stderr,"acgram: input file too small (%dms) for requested start and length \n", frames * 1000 / samplerate );
	exit(1);
    }

    framebytes = frameheight * frameslength * sizeof(short) ;
    startbytes = frameheight * start  * sizeof(short) ;


    /**** Output frame dimensions ****/

    Newframeheight =  frameheight ;
    Newframewidth  =  to_p( lagstr, samplerate );   /* Max acf lag   */
    Newframeshift  =  frameshift ;
    Newframebytes  =  Newframeheight * Newframewidth * sizeof(short) ;


    /**** Padding and output scale factor ****/

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

    if ( normalize == 0 )   /* no normalize, so divide by length to get average of products */
	scale = 1./(double)(framewidth + zeroes) * atof( scalestr ) ;
    else
	scale = atof( scalestr ) ;


    /**** Debug ****/

    if ( ison( debugstr ) ) {
	fprintf(stderr, "Input:   frames=%d frameheight=%d frameshift=%d\n", frames, frameheight, frameshift ) ;
	fprintf(stderr, "Sizes:   start=%d length=%d framewidth=%d frameslength=%d\n", start, length, framewidth, frameslength ) ;
	fprintf(stderr, "Output:  zeroes=%d Newframewidth=%d Newframes=%d \n", zeroes, Newframewidth, Newframes ) ;
	exit( 1 ) ;
    }


    /**** Allocate space (input file, window coeffs, and fft working space) ****/

    if ( ( file = (short *)malloc( framebytes ) ) == NULL ) {
	fprintf(stderr,"acgram: malloc out of space\n");
	exit(1);
    }

    if ( window )
	W = hamming( framewidth ) ;

    buf = (float *)malloc( ( framewidth + zeroes ) * sizeof(float) );


    /**** Write new sai header ****/

    if ( ison( headerstr ) )
	WriteHeader( SaiHeader(header), stdout );


    /**** Seek to start of input file in blocks of framebytes ****/

    for (i=framebytes ; i < startbytes ; i += framebytes)
	if ( fread( file, framebytes, 1, fp ) == NULL ) {
	    fprintf(stderr,"acgram: missing data after header\n");
	    exit(1);
	}
    if ( (startbytes -= (i - framebytes)) > 0 )
	if ( fread( file, startbytes, 1, fp ) == NULL ) {
	    fprintf(stderr,"acgram: missing data after header\n");
	    exit(1);
	}


    /**** Read framebytes of i/p file ****/

    if ( fread( file, framebytes, 1, fp ) == NULL ) {
	fprintf(stderr,"acgram: missing data after header\n");
	exit(1);
    }


    /**** Process frames ****/

    framebyteshift = frameshift * frameheight ;
    endframe = file + Newframes * framebyteshift ;

    for ( frame = file, i=1 ; frame < endframe ; frame += framebyteshift, i++ ) {
	if ( Verbose )
	    fprintf(stderr,"acgram: %d / %d  [%dms]\n", i, Newframes, start*1000/samplerate );
	if ( ( f = acgram( frame ) ) < newscale )
	    newscale = f ;
	start += frameshift ;
    }

    fclose(fp);

    if ( newscale < 1. ) {
	fprintf( stderr, "Warning: 16-bit overflow during acgram.  " ) ;
	if ( normalize == 0 )
	    fprintf( stderr, "Try scale<%f\n", newscale * (float)(framewidth + zeroes) ) ;
	else
	    fprintf( stderr, "Try scale<%f\n", newscale ) ;
    }

    if ( Verbose )
	fprintf(stderr,"acgram: done\n" ) ;
}


/************************** Autocorrelation function ***********************/

/* Call acf for each row in the SAI frame */

float acgram( frame )
short *frame ;
{
    register int  i, j, row, col ;
    float    writebuf(), f, newscale = 1. ;

    for ( row=0 ; row < frameheight ; row++ ) {

	if ( window )
	    for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
		buf[col] = frame[i] * W[col] ;
	else
	    for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
		buf[col] = frame[i] ;

	for ( j=0 ; j<zeroes ; j++ )    /* padding */
	    buf[col++] = 0 ;

	acf( buf, framewidth+zeroes ) ;

	if ( normalize == 0 )
	    f = writebuf( buf, Newframewidth, scale ) ;
	else
	    f = writebuf( buf, Newframewidth, scale/buf[0] ) ;

	if ( f < newscale )
	    newscale = f ;
    }

    return ( newscale ) ;
}


/********************** Write o/p ****************************************/

float writebuf( buf, n, scale )       /* write buffer as scaled shorts */
float *buf ;
int    n ;
float  scale ;
{
    register int  i ;
    short         p ;
    float newscale  ;

    if ( !negative ) {  /* usual acf (ie. positive half, with zeroth lag on left) */
	for (i=0 ; i < n ; i++)  {

	    newscale = check_overflow( buf[i], scale, 1 ) ;

	    p = (short)( buf[i] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}
    }

    else {      /* write negative half of acf (zeroth lag on right) */
	for ( i = n-1 ; i >= 0 ; --i )  {

	    newscale = check_overflow( buf[i], scale, 1 ) ;

	    p = (short)( buf[i] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}
    }
    return ( newscale ) ;
}



/************************ Output header ************************************/
/*
   First create a new header with
     nwidth_aid
     pwidth_aid
     frstep_aid
   Then copy the original nap header into the sai header, changing in order:
     frames
     frameshift
     framewidth
     frameheight
     framebytes
     animate
     display
   Then change applic name [gennap] to [gensai] in the Version string.
   Finally, update the new header_bytes, and return the new header.
*/

char *SaiHeader( napheader )
char *napheader ;
{
    char *saiheader;
    char *p0, *p1, *p2, *s, str[64];

    saiheader = (char *)malloc( strlen(napheader) + 128 ) ;

    p0 = saiheader ;
    p1 = napheader ;


    /** copy initial header bytes **/

    p2 = HeaderString( napheader , "header_bytes" ) ;
    while( p1 < p2 )
	*p0++ = *p1++ ;
    while (*p1 != '\n')
	*p0++ = *p1++ ;
    *p0++ = *p1++ ;


    /** insert nwidth_aid at top of new header **/

    sprintf(str,"nwidth_aid=0\n");
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;


    /** insert pwidth_aid into new header **/

    sprintf(str,"pwidth_aid=%dms\n", (int)(1000*((float)Newframewidth/samplerate)) );
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;


    /** insert frstep_aid into new header **/

    sprintf(str,"frstep_aid=%dms\n", (int)(1000*((float)frameshift/samplerate)) );
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;


    /** copy up to frames **/

    p2 = HeaderString( napheader , "frames" ) ;
    while( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%d\n", Newframes);
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy up to frameshift **/

    p2 = HeaderString( napheader , "frameshift" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%d\n", Newframeshift);
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy up to framewidth **/

    p2 = HeaderString( napheader , "framewidth" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%d\n", Newframewidth);
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy up to frameheight **/

    p2 = HeaderString( napheader , "frameheight" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%d\n", Newframeheight);
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy up to framebytes **/

    p2 = HeaderString( napheader , "framebytes" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%d\n", Newframebytes);
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy up to animate_ctn **/

    p2 = HeaderString( napheader , "animate_ctn" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%s\n", anistr );
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy up to display **/

    p2 = HeaderString( napheader , "display" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    sprintf(str,"%s\n", dispstr );
    for (s = str ; *s != '\n' ; )
	*p0++ = *s++;
    *p0++ = *s;
    while (*p1 != '\n')
	*p1++;
    *p1++;


    /** copy rest of header **/

    p2 = HeaderString( napheader , "Version" ) ;
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    while (*p1 != ']')          /* change applic name [gennap] to [gensai] */
	*p0++ = *p1++ ;
    *(p0-3) = 's' ;
    *(p0-2) = 'a' ;
    *(p0-1) = 'i' ;
    while (*p1 != '\n')
	*p0++ = *p1++ ;
    *p0++ = *p1++ ;


    /** update header_bytes **/

    sprintf(str, "%0*d", 7, p0-saiheader);
    p0 = HeaderString( saiheader , "header_bytes" ) ;
    s = str;
    while(*p0 != '\n')
	*p0++ = *s++ ;


    return saiheader;
}