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