comparison 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;
+}
+
+