annotate tools/audim.c @ 0:5242703e91d3 tip

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

/*
  audim.c           auditory images
 ---------

   This program provides methods of constructing time-varying auditory
   images from the output of the cochlear model which are alternatives to
   the gensai (stabilized auditory image) program. Correlograms (ie. row-
   wise autocorrelation) and row-wise Fourier transform are provided.  The
   program reads an AIM header and the output from the genbmm (basilar mem-
   brane motion) or gennap (neural activity pattern) programs. This is
   divided into contiguous time frames and written on the stdout with an
   appropriate header as if output from the gensai program. This enables
   genbmm or gennap output to be divided into time frames and replayed as a
   cartoon by gensai using the "useprevious" option.  Additional processing
   to each frame is optionally available to compute alternative forms of
   auditory image according to the "image" option.


   The options "start" and "length" specify the size of the input.  The
   options "width" and "frstep"  specify the frames which are output.  The
   input is divided into frames according to the width option and the
   frstep option.

   Special option values are:

     length=max      Read all the input from the given start to its end.
     width=max       Output framewidth is set equal to the given input
	       length, and if this is also "max", then the
	       framewidth is the remainder of the input.
     frames=1-max    Select the 1st to the last frame inclusively.
     frame=4         Select the 4th frame only.

     image=off       Divide the input into frames for replay as a cartoon.
     image=acgram    Compute correlogram of each frame.
     image=ftgram    Compute power spectrum of each channel of each frame.
     image=phgram    Compute phase spectrum of each channel of each frame.

   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 input nap is divided into frames according to the width option (called
  framewidth internally) and the frstep option (called frameshift internally).
  The framewidth is used as the basis for output image frame width:
  the output width is actually the next power-of-2 above twice the framewidth.
  This is because:
  a. fft and acf output is always half the given frame size, due to symmetry,
     so the input must be padded to twice its width.
  b. To use fft, the input must also be a power of two.
  So, each input row is padded out to the next power of two above twice its
  width.

  Certain display parameters have different default values for different
  applications. The gensai display parameters should be set to the
  appropriate values, in order to plot the cartoon on the same scale. For
  example: when the source application is gennap, set gensai top=1000,
  when the source application is genbmm, set gensai bottom=-100.


Examples:

1. An animated normalized autocorrelogram from a nap:

gennap len=64ms output=stdout display=off file | audim norm=on anim=on > file.sai
gensai useprev=on top=1000 cegc             -(for landscape plot)
genspl useprev=on top=1000 pensize=2 cegc   -(for spiral plot)


2. To convert a nap to multiple animated frames:

gennap len=16ms display=off output=stdout file | audim image=off width=8ms frstep=0.2ms anim=on > file.sai
gensai useprev=on top=1000 cegc             -(for landscape plot)
genspl useprev=on top=1000 pensize=2 cegc   -(for spiral plot)

    (Note: spirals look better in a square box, so you might use options:
	   dencf=1 width=500 height=500 )

3. To view the basilar membrane from a cross section, animating the waves on it.

genbmm mincf=220 maxcf=660 len=8ms output=stdout display=off file | audim image=off width=1p frstep=1p display=on anim=on | edframe Tran=on > foo.sai
gensai bott=-100 useprev=on mag=.2 foo

or:

genbmm mincf=220 maxcf=660 len=32ms output=stdout display=off file | \
audim image=off width=1p frstep=1p display=on anim=on | \
edframe Tran=on Header=off > foo

x11play -n75 -a500 foo

*/


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

char applic[]     = "Auditory images of neural activity patterns (NAP). " ;

static char *helpstr, *debugstr, *dispstr, *anistr, *startstr, *imstr ;
static char *lengthstr, *widthstr, *shiftstr, *headerstr, *cfstr, *framestr ;
static char *normstr, *wstr, *sacfstr, *smagstr, *sphstr ;
static char *blstr, *erbstr ;

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     },
    {   "image"     ,   "acgram"    ,  &imstr       ,   "Form of auditory image."                    , VAL     },
    {   "header"    ,   "on"        ,  &headerstr   ,   "Header (for gensai useprevious=on)."        , VAL     },
    {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p NAP"                     , VAL     },
    {   "length"    ,   "max"       ,  &lengthstr   ,   "Length of i/p NAP to process."              , VAL     },
    {   "frstep"    ,   "16ms"      ,  &shiftstr    ,   "Step between frames of the auditory image." , VAL     },
    {   "width"     ,   "32ms"      ,  &widthstr    ,   "Width of auditory image."                   , VAL     },
    {   "normalize" ,   "off"       ,  &normstr     ,   "Normalization"                              , SETFLAG },
    {   "window"    ,   "off"       ,  &wstr        ,   "Hamming window"                             , SETFLAG },
    {   "bandlimit" ,   "off"       ,  &blstr       ,   "Zero outside bandwidth."                    , SETFLAG },
    {   "erbs"      ,   "2"         ,  &erbstr      ,   "Width (in ERBS) of band either side of cf." , SVAL    },
    {   "align_cf"  ,   "off"       ,  &cfstr       ,   "Align on centre frequency."                 , SETFLAG },
    {   "scale_acf" ,   "0.025"     ,  &sacfstr     ,   "Scale factor for acf output"                , SVAL    },
    {   "scale_fft" ,   "0.08"      ,  &smagstr     ,   "Scale factor for fft magnitude output"      , SVAL    },
    {   "scale_phs" ,   "100"       ,  &sphstr      ,   "Scale factor for fft phase output"          , SVAL    },
   ( char * ) 0 } ;


int     frameheight ;   /* Nap parameters read from header */
int     frames, samplerate ;
int     framewidth  ;   /* width option */

int     start, length ;

int     Newframeheight, Newframewidth ; /* Sai parameters to write  */
int     Newframes ;
int     Newframeshift ;
int     Newframebytes ;


char   *limitstr ;      /* Sai parameters used to get centre-frequencies */
char   *qualstr  ;
char   *minstr   ;
char   *maxstr   ;
char   *denstr   ;
char   *chansstr ;
double *frequencies ;
int     nerbs ;         /* bandwidth defined as number of erbs either side of cf */
int     bandlimit ;     /* flag for zeroing outside bandwidth around centre-freqs */
int     align ;         /* flag for aligning centre-freqs */

#define NONE    0       /* forms of auditory image */
#define ACGRAM  1
#define FTGRAM  2
#define PHGRAM  3
int     image ;

float  *buf ;           /* fft working space */
int     zeroes ;        /* padding */
float  *W ;             /* Hamming window for fft */
short   window ;        /* flag for window */

short   debug ;         /* flag for debug printf's */

main(argc, argv)
int   argc ;
char *argv[] ;
{
    FILE   *fp ;
    int     i, framebytes, startbytes, frameshift, framebyteshift;
    char   *header, *SaiHeader();
    short  *nap, *frame, *endframe;
    int     a, b ;
    char   *val1, *val2 ;

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

    window = ison( wstr ) ;
    nerbs  = atoi( erbstr ) ;
    debug  = ison( debugstr ) ;

    /**** Form of image ****/

    if ( iststr( imstr, "acgram" ) )  {
	image = ACGRAM ;
    }
    else if ( iststr( imstr, "ftgram" ) )  {
	image = FTGRAM ;
    }
    else if ( iststr( imstr, "phgram" ) )  {
	image = PHGRAM ;
    }
    else if ( isoff( imstr ) ) {
	image = NONE ;
    }
    else {
	fprintf(stderr, "audim: unknown form of auditory image [%s]\n", imstr ) ;
	exit( 1 ) ;
    }

    align     = ison( cfstr ) ;
    bandlimit = ison( blstr ) ;


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

    if ( getvals( framestr, &val1, &val2, "-" ) == BADVAL ) {
	fprintf(stderr,"audim: bad frame selector [%s]\n", framestr ) ;
	exit( 1 ) ;
    }
    a = atoi( val1 ) ;
    if ( isempty( val2 ) )    b = a ;
    else if ( ismax( val2 ) ) b = 0 ;
    else                      b = atoi( val2 ) ;

    if (b<a && b>0) fprintf(stderr,"audim warning: bad frame specifiers\n");


    /**** Read header and options to find dimensions of auditory image ****/

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

    frameheight = HeaderInt( header, "frameheight" );
    frames      = HeaderInt( header, "frames"      );
    samplerate  = HeaderSamplerate( header );
    start       = to_p( startstr,  samplerate );
    frameshift  = to_p( shiftstr, samplerate ) ;

    if ( ismax( lengthstr ) )   /* Special case for remainder of input */
	length = frames - start ;
    else
	length = to_p( lengthstr, samplerate );

    if ( start + length > frames ) {
	fprintf(stderr,"audim: nap too small (%d ms) for requested start and length \n", to_ms( frames, samplerate ) );
	exit(1);
    }

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

    if ( ismax( widthstr ) )    /* Special case for single frame of max width */
	framewidth = length ;
    else
	framewidth = to_p( widthstr, samplerate ) ;

    frames = ( length - framewidth ) / frameshift ;
    zeroes = getpower( framewidth << 1 ) - framewidth ;


    /**** Get centre-frequencies information using info given in header ****/

    limitstr =    HeaderStringOnly( header, "bwmin_afb"    );
    qualstr =     HeaderStringOnly( header, "quality_afb"  );
    minstr =      HeaderStringOnly( header, "mincf_afb"    );
    maxstr =      HeaderStringOnly( header, "maxcf_afb"    );
    denstr =      HeaderStringOnly( header, "dencf_afb"    );
    chansstr =    HeaderStringOnly( header, "channels_afb" );

    SetErbParameters( to_Hz( limitstr ), atof( qualstr ) ) ;
    if( OptionInt( chansstr ) == 0 )
	frequencies = GenerateCenterFrequencies( to_Hz( minstr ), to_Hz( maxstr ), atof( denstr ) )   ;
    else
	frequencies = NumberedCenterFrequencies( to_Hz( minstr ), to_Hz( maxstr ), OptionInt( chansstr ) ) ;


    /**** Output frame dimensions (to write into new header) ****/

    Newframeheight = frameheight ;

    if ( !align ) {
	if ( image == NONE )
	    Newframewidth = framewidth ;
	else
	    Newframewidth = ( framewidth + zeroes ) >> 1 ;  /* spectral bins */
    }
    else  {      /* number of bins corresponding to the filter bandwidth */
	switch ( image ) {
	    case NONE   :   Newframewidth = framewidth ;
			    break ;
	    case ACGRAM :   Newframewidth = maxlagwidth( to_Hz( minstr ), nerbs ) ;
			    break ;
	    case FTGRAM :   Newframewidth = maxbandwidth( to_Hz( maxstr ), nerbs, (framewidth+zeroes)>>1 ) ;
			    break ;
	}
    }
    Newframeshift =  frameshift ;
    Newframebytes =  Newframeheight * Newframewidth * sizeof(short) ;

    if ( frameshift > 0 ) {
	if ( b == 0 )
	    Newframes = frames - (a-1) ;
	else
	    Newframes = b - (a-1) ;
    }
    else {
	fprintf(stderr,"audim: non-positive frstep [%d]\n", frameshift);
	exit(1);
    }


    /**** Check limits etc.. ****/

    if ( b > frames ) {
	fprintf(stderr,"audim: can't select frame %d out of %d frames\n", b, frames ) ;
	exit(1);
    }

    if ( length < framewidth ) {
	fprintf(stderr,"audim: nap too small (%d ms) for requested width\n", to_ms( length, samplerate ) );
	exit(1);
    }

    if (frames<=0) {
	if (frames==0)  fprintf(stderr,"audim: zero frames input\n");
	if (frames<0)   fprintf(stderr,"audim: garbled number of frames, (start set too big?)\n");
	exit(1);
    }

    if (a>frames || b>frames)
	fprintf(stderr,"audim warning: bad frame selectors\n");



    if ( debug ) {
	fprintf(stderr, "Header:  frames=%d frameheight=%d\n", frames, frameheight ) ;
	fprintf(stderr, "Options: a=%d b=%d start=%d length=%d frameshift=%d framewidth=%d\n", a, b, start, length, frameshift, framewidth ) ;
	fprintf(stderr, "Output:  zeroes=%d Newframewidth=%d Newframes=%d \n", zeroes, Newframewidth, Newframes ) ;
    }


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

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


    /**** Allocate space for framebytes of nap data ****/

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


    /**** Allocate and assign window coeffs ****/

    if ( window )
	W = hamming( framewidth ) ;


    /**** Allocate fft working space ****/

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


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

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

    /**** Read framebytes of i/p nap data ****/

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

    framebyteshift = frameshift * frameheight ;
    nap = nap + (a-1) * framebyteshift ;
    endframe = nap + Newframes * framebyteshift ;

    switch ( image ) {
	case NONE   : fprintf(stderr,"Output framewidth = %dp \n", Newframewidth ) ;  break ;
	case ACGRAM : fprintf(stderr,"Output framewidth = %dp Lag range = 0-%dms\n", Newframewidth, (Newframewidth*1000)/samplerate ) ; break ;
	case FTGRAM : fprintf(stderr,"Output framewidth = %dp Frequency range = 0-%dkHz\n", Newframewidth, samplerate/2000 ) ; break ;
	case PHGRAM : break ;
    }

    for ( frame = nap ; frame < endframe ; frame += framebyteshift ) {

	fprintf(stderr,"audim: %d / %d\n", a++, frames );
	switch ( image ) {

	    case NONE   : writeframe( frame ) ; break ;
	    case ACGRAM : acgram( frame ) ;     break ;
	    case FTGRAM : ftgram( frame ) ;     break ;
	    case PHGRAM : phgram( frame ) ;     break ;

	}
    }
    fprintf(stderr,"audim: done\n" );
}





int OptionInt( str )
char *str ;
{
    if( strcmp( str, "on" ) == 0 )
	return( 1 ) ;
    else if( strcmp( str, "Not_used" ) == 0 )
	return( 0 ) ;
    else
	return( atoi( str ) ) ;
}


/*********************** Read header and build new header *****************/

/*
   Copy the original nap header to a new sai header, changing in order:
     frames
     frameshift
     framewidth
     frameheight
     framebytes
     animate
     display
   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) + 256 ) ;

    p0 = saiheader ;
    p1 = napheader ;

    /* copy to first item after the header_bytes */

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

    /* insert frstep_aid at top of new header */

    sprintf( str,"frstep_aid=%s\n", shiftstr ) ;
    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++ ;

    if ( ( p2 = ApplicString( napheader ) ) == (char *)0 )  {
	fprintf(stderr,"audim: application name not found in header\n");
	exit( 1 ) ;
    }
    while ( p1 < p2 )
	*p0++ = *p1++ ;

    *p0++ = 's' ; *p0++ = 'a' ; *p0++ = 'i' ;   /* copy sai into applic name */
    p1 += 3 ;

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


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

/* Call acf for each row in auditory image */

acgram( frame )
short *frame ;
{
    static   float  scale ;
    static   int    normalize ;
    static   int    first = 1 ;
    register int    i, j, row, col ;
    short    p ;

    if (first) {
	normalize = ison( normstr ) ;
	scale = 1./(double)( framewidth + zeroes ) * atof( sacfstr ) ;
	first=0;
    }

    for ( row=0 ; row < frameheight ; row++ ) {
	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 ( align || bandlimit )
	    writelags( buf, row, (framewidth+zeroes)>>1, Newframewidth, scale ) ;
	else {
	    if ( normalize == 0 )
		writebuf( buf, Newframewidth, scale ) ;
	    else
		writebuf( buf, Newframewidth, 5000./buf[0] ) ;
	}
    }
}


/********************************** FFT ************************************/

/* Call fft and magnitude spectrum for each row in auditory image */

ftgram( frame )
short *frame ;
{
    static   float  scale ;
    static   int    first = 1 ;
    register int    i, j, row, col ;

    if (first) {
	scale = atof( smagstr ) ;
	first=0;
    }

    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 ;

	fft( buf, framewidth+zeroes, 1 ) ;
	mag( buf, framewidth+zeroes ) ;

	if ( align || bandlimit )
	    writebins( buf, row, (framewidth+zeroes)>>1, Newframewidth, scale ) ;
	else
	    writebuf( buf, Newframewidth, scale ) ;
    }
}


/********************************** Phase **********************************/

/* Call fft and phase spectrum for each row in auditory image */

phgram( frame )
short *frame ;
{
    static   float  scale ;
    static   int    first = 1 ;
    register int    i, j, row, col ;

    if (first) {
	scale = atof( sphstr ) ;
	first=0;
    }

    for ( row=0 ; row < frameheight ; row++ ) {
	for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
	    buf[col] = frame[i] ;
	for ( j=0 ; j<zeroes ; j++ )    /* padding */
	    buf[col++] = 0 ;

	fft( buf, framewidth+zeroes, 1 ) ;
	phase( buf, framewidth+zeroes ) ;

	writebuf( buf, Newframewidth, scale ) ;
    }
}


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

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

    for (i=0 ; i < n ; i++)  {
	p = (short)( buf[i] * scale ) ;
	fwrite( &p, sizeof(short), 1, stdout ) ;
    }
}




writeframe( frame )     /* Write a nap frame in sai format */
short  *frame ;
{
    int    i, row, col;

    for (row=0 ; row < frameheight ; row++)
	for (col=0 , i=row ; col < framewidth ; col++, i+=frameheight)
	    fwrite( &frame[i], sizeof(short), 1, stdout );
}


/*
  Write a band of spectral bins around the centre-freq of the given channel.

  The o/p spectrum is 0 to samplerate/2 Hz, resolved into spectral bins:
	bins = (framewidth + zeroes) / 2

  For the i'th channel, the centre-frequency cf = frequencies[i].
  Filters are symmetrical bandpass, so considering a bandwidth of 4 ERB
  the side-frequencies at the band limits are:

    side1 = cf - 2*ERB     [Hz]
    side2 = cf + 2*ERB     [Hz]

  Note: side1 < side2

  Convert frequency to bin number as follows:
  Let centre-frequency cf have bin number j, so that ratios:
	cf/(samplerate/2) = j/bins
  and so:
	j  = ( cf*bins ) / ( samplerate/2 )

  (Similarly for each of the side-frequencies at the 2-ERB band limits).

*/

writebins( buf, chan, bins, maxwidth, scale )
float *buf ;
int    chan ;
int    bins ;       /* number of spectral bins */
int    maxwidth ;
float  scale ;
{
    double ERB, cf, side1, side2, max ;
    register int  i, j, s1, s2 ;
    short         p ;

    ERB   = frequencies[chan]/9.26449 + 24.7 ;
    cf    = frequencies[chan]          ;
    side1 = frequencies[chan] - nerbs*ERB  ;
    side2 = frequencies[chan] + nerbs*ERB  ;

    j  = (int)( ( cf    * bins ) / ( samplerate >> 1 ) + 0.5 ) ;   /* rounding */
    s1 = (int)( ( side1 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;
    s2 = (int)( ( side2 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;

    if ( debug )
	fprintf(stderr,"chan=%d bins=%d maxwidth=%d j=%d side1=%d side2=%d\n",chan,bins,maxwidth,j,s1,s2);

    /*  Plot the given bandwidth, zeroing elsewhere:     */

    if ( bandlimit && !align ) {

	for ( i=0, p=0 ; i < s1 ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;

	for (  ; i < s2 && i < bins ; i++ )  {
	    p = (short)( buf[i] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}

	for ( p=0 ; i < bins ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;
    }

    /* Plot the given bandwidth, shifting it to the centre */

    if ( align ) {

/* align on max bin in band
	max = 0 ;
	for ( i=s1 ; i<s2 ; i++ )
	    if ( buf[i] >= max ) {
		max = buf[i] ;
		j = i ;
	    }
	for ( i=0, p=0 ; i < s1+(maxwidth/2-j) ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	for ( j = s1 ; j < s2 && i < maxwidth ; i++, j++ )  {
	    p = (short)( buf[j] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}
	for ( p=0 ; i < maxwidth ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;
*/

	j = (maxwidth+s1-s2)/2 ;
	for ( i=0, p=0 ; i < j ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;

	for ( j = s1 ; j < s2 ; i++, j++ )  {
	    p = (short)( buf[j] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}
	for ( p=0 ; i < maxwidth ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;

    }
}



/*
  Write a band of acf lags around the centre-period of the given channel.

  For the i'th channel, the centre-frequency
	Centre-frequency   cf  = frequencies[i]     [Hz]
	Bandwidth          ERB = cf/9.26449 + 24.7  [Hz]
	Centre-period       cp = samplerate / cf    [samples]

  Filters are symmetrical bandpass, so considering a bandwidth of 4 ERB

    side1 = samplerate / ( cf - 2*ERB )    [samples]
    side2 = samplerate / ( cf + 2*ERB )    [samples]

  Note: side2 < side1 since periods

  Also, successive bins ARE successive lags (in samples).
*/



writelags( buf, chan, lags, maxwidth, scale )
float *buf ;
int    chan ;
int    lags ;       /* number of acf lags */
int    maxwidth ;
float  scale ;
{
    double ERB, cp, side1, side2 ;
    register int  i, j, s1, s2 ;
    short         p ;

    ERB   = frequencies[chan]/9.26449 + 24.7 ;
    cp    = (double)samplerate /   frequencies[chan]           ;
    side1 = (double)samplerate / ( frequencies[chan] - nerbs*ERB ) ;
    side2 = (double)samplerate / ( frequencies[chan] + nerbs*ERB ) ;

    j  = (int)( cp    + 0.5 ) ;         /* rounding */
    s1 = (int)( side1 + 0.5 ) ;
    s2 = (int)( side2 + 0.5 ) ;

    if ( debug )
	fprintf(stderr,"chan=%d lags=%d side1=%d side2=%d\n",chan,lags,s1,s2);

    /*  Plot the given bandwidth, zeroing elsewhere:     */

    if ( bandlimit && !align ) {

	for ( i=0, p=0 ; i < s2 ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;

	for (  ; i < s1 && i < lags ; i++ )  {
	    p = (short)( buf[i] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}

	for ( p=0 ; i < lags ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;
    }

    /* Plot the given bandwidth, shifting it to the centre */

    if ( align ) {

	j = (maxwidth+s2-s1)/2 ;
	for ( i=0, p=0 ; i < j ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;

	for ( j = s2 ; j < s1 ; i++, j++ )  {
	    p = (short)( buf[j] * scale ) ;
	    fwrite( &p, sizeof(short), 1, stdout ) ;
	}
	for ( p=0 ; i < maxwidth ; i++ )
	    fwrite( &p, sizeof(short), 1, stdout ) ;
    }
}


/*
  Return the number of lags corresponding to a bandwidth of nerbs ERBS either
  side of a given mincf. This will be the max required framewidth.
*/

maxlagwidth( mincf, nerbs )
double mincf ;
int    nerbs ;
{
    double ERB, side1, side2 ;

    ERB   = mincf/9.26449 + 24.7 ;
    side1 = (double)samplerate / ( mincf - nerbs*ERB ) ;
    side2 = (double)samplerate / ( mincf + nerbs*ERB ) ;

    return (int)( side1 - side2 ) ;
}


/*
  Return the number of bins corresponding to a bandwidth of nerbs ERBS either
  side of a given maxcf. This will be the max required framewidth.
*/

maxbandwidth( maxcf, nerbs, bins )
double maxcf ;
int    nerbs ;
int    bins  ;
{
    double ERB, side1, side2 ;
    int    s1, s2 ;

    ERB   = maxcf/9.26449 + 24.7 ;
    side1 = maxcf - nerbs*ERB  ;
    side2 = maxcf + nerbs*ERB  ;

    /* convert frequencies to bin numbers */

    s1 = (int)( ( side1 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;
    s2 = (int)( ( side2 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;

    return (int)( s2 - s1 ) ;
}