tomwalters@0: /*
tomwalters@0:   audim.c           auditory images
tomwalters@0:  ---------
tomwalters@0: 
tomwalters@0:    This program provides methods of constructing time-varying auditory
tomwalters@0:    images from the output of the cochlear model which are alternatives to
tomwalters@0:    the gensai (stabilized auditory image) program. Correlograms (ie. row-
tomwalters@0:    wise autocorrelation) and row-wise Fourier transform are provided.  The
tomwalters@0:    program reads an AIM header and the output from the genbmm (basilar mem-
tomwalters@0:    brane motion) or gennap (neural activity pattern) programs. This is
tomwalters@0:    divided into contiguous time frames and written on the stdout with an
tomwalters@0:    appropriate header as if output from the gensai program. This enables
tomwalters@0:    genbmm or gennap output to be divided into time frames and replayed as a
tomwalters@0:    cartoon by gensai using the "useprevious" option.  Additional processing
tomwalters@0:    to each frame is optionally available to compute alternative forms of
tomwalters@0:    auditory image according to the "image" option.
tomwalters@0: 
tomwalters@0: 
tomwalters@0:    The options "start" and "length" specify the size of the input.  The
tomwalters@0:    options "width" and "frstep"  specify the frames which are output.  The
tomwalters@0:    input is divided into frames according to the width option and the
tomwalters@0:    frstep option.
tomwalters@0: 
tomwalters@0:    Special option values are:
tomwalters@0: 
tomwalters@0:      length=max      Read all the input from the given start to its end.
tomwalters@0:      width=max       Output framewidth is set equal to the given input
tomwalters@0: 	       length, and if this is also "max", then the
tomwalters@0: 	       framewidth is the remainder of the input.
tomwalters@0:      frames=1-max    Select the 1st to the last frame inclusively.
tomwalters@0:      frame=4         Select the 4th frame only.
tomwalters@0: 
tomwalters@0:      image=off       Divide the input into frames for replay as a cartoon.
tomwalters@0:      image=acgram    Compute correlogram of each frame.
tomwalters@0:      image=ftgram    Compute power spectrum of each channel of each frame.
tomwalters@0:      image=phgram    Compute phase spectrum of each channel of each frame.
tomwalters@0: 
tomwalters@0:    Most options in the input header are copied to the output header.  This
tomwalters@0:    enables options which are needed for the eventual display to pass
tomwalters@0:    straight through. Some options are set so that they can override the
tomwalters@0:    input header. For example, the display option is set on to enable
tomwalters@0:    display even when input has display=off. The animate option can be set
tomwalters@0:    on even when the input has animate=off.  Parts of the header are changed
tomwalters@0:    for the new sai format: (frames, frameshift, framewidth, frameheight,
tomwalters@0:    framebytes).
tomwalters@0: 
tomwalters@0:   Padding:
tomwalters@0:   The input nap is divided into frames according to the width option (called
tomwalters@0:   framewidth internally) and the frstep option (called frameshift internally).
tomwalters@0:   The framewidth is used as the basis for output image frame width:
tomwalters@0:   the output width is actually the next power-of-2 above twice the framewidth.
tomwalters@0:   This is because:
tomwalters@0:   a. fft and acf output is always half the given frame size, due to symmetry,
tomwalters@0:      so the input must be padded to twice its width.
tomwalters@0:   b. To use fft, the input must also be a power of two.
tomwalters@0:   So, each input row is padded out to the next power of two above twice its
tomwalters@0:   width.
tomwalters@0: 
tomwalters@0:   Certain display parameters have different default values for different
tomwalters@0:   applications. The gensai display parameters should be set to the
tomwalters@0:   appropriate values, in order to plot the cartoon on the same scale. For
tomwalters@0:   example: when the source application is gennap, set gensai top=1000,
tomwalters@0:   when the source application is genbmm, set gensai bottom=-100.
tomwalters@0: 
tomwalters@0: 
tomwalters@0: Examples:
tomwalters@0: 
tomwalters@0: 1. An animated normalized autocorrelogram from a nap:
tomwalters@0: 
tomwalters@0: gennap len=64ms output=stdout display=off file | audim norm=on anim=on > file.sai
tomwalters@0: gensai useprev=on top=1000 cegc             -(for landscape plot)
tomwalters@0: genspl useprev=on top=1000 pensize=2 cegc   -(for spiral plot)
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 2. To convert a nap to multiple animated frames:
tomwalters@0: 
tomwalters@0: gennap len=16ms display=off output=stdout file | audim image=off width=8ms frstep=0.2ms anim=on > file.sai
tomwalters@0: gensai useprev=on top=1000 cegc             -(for landscape plot)
tomwalters@0: genspl useprev=on top=1000 pensize=2 cegc   -(for spiral plot)
tomwalters@0: 
tomwalters@0:     (Note: spirals look better in a square box, so you might use options:
tomwalters@0: 	   dencf=1 width=500 height=500 )
tomwalters@0: 
tomwalters@0: 3. To view the basilar membrane from a cross section, animating the waves on it.
tomwalters@0: 
tomwalters@0: 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
tomwalters@0: gensai bott=-100 useprev=on mag=.2 foo
tomwalters@0: 
tomwalters@0: or:
tomwalters@0: 
tomwalters@0: genbmm mincf=220 maxcf=660 len=32ms output=stdout display=off file | \
tomwalters@0: audim image=off width=1p frstep=1p display=on anim=on | \
tomwalters@0: edframe Tran=on Header=off > foo
tomwalters@0: 
tomwalters@0: x11play -n75 -a500 foo
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "header.h"
tomwalters@0: #include "options.h"
tomwalters@0: #include "units.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: #include "sigproc.h"
tomwalters@0: #include "freqs.c"
tomwalters@0: 
tomwalters@0: char applic[]     = "Auditory images of neural activity patterns (NAP). " ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *dispstr, *anistr, *startstr, *imstr ;
tomwalters@0: static char *lengthstr, *widthstr, *shiftstr, *headerstr, *cfstr, *framestr ;
tomwalters@0: static char *normstr, *wstr, *sacfstr, *smagstr, *sphstr ;
tomwalters@0: static char *blstr, *erbstr ;
tomwalters@0: 
tomwalters@0: static Options option[] = {
tomwalters@0:     {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                       , DEBUG   },
tomwalters@0:     {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                           , DEBUG   },
tomwalters@0:     {   "display"   ,   "on"        ,  &dispstr     ,   "Display output image."                      , SETFLAG },
tomwalters@0:     {   "animate"   ,   "off"       ,  &anistr      ,   "Animate cartoon."                           , SETFLAG },
tomwalters@0:     {   "frames"    ,   "1-max"     ,  &framestr    ,   "Select frames inclusively"                  , VAL     },
tomwalters@0:     {   "image"     ,   "acgram"    ,  &imstr       ,   "Form of auditory image."                    , VAL     },
tomwalters@0:     {   "header"    ,   "on"        ,  &headerstr   ,   "Header (for gensai useprevious=on)."        , VAL     },
tomwalters@0:     {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p NAP"                     , VAL     },
tomwalters@0:     {   "length"    ,   "max"       ,  &lengthstr   ,   "Length of i/p NAP to process."              , VAL     },
tomwalters@0:     {   "frstep"    ,   "16ms"      ,  &shiftstr    ,   "Step between frames of the auditory image." , VAL     },
tomwalters@0:     {   "width"     ,   "32ms"      ,  &widthstr    ,   "Width of auditory image."                   , VAL     },
tomwalters@0:     {   "normalize" ,   "off"       ,  &normstr     ,   "Normalization"                              , SETFLAG },
tomwalters@0:     {   "window"    ,   "off"       ,  &wstr        ,   "Hamming window"                             , SETFLAG },
tomwalters@0:     {   "bandlimit" ,   "off"       ,  &blstr       ,   "Zero outside bandwidth."                    , SETFLAG },
tomwalters@0:     {   "erbs"      ,   "2"         ,  &erbstr      ,   "Width (in ERBS) of band either side of cf." , SVAL    },
tomwalters@0:     {   "align_cf"  ,   "off"       ,  &cfstr       ,   "Align on centre frequency."                 , SETFLAG },
tomwalters@0:     {   "scale_acf" ,   "0.025"     ,  &sacfstr     ,   "Scale factor for acf output"                , SVAL    },
tomwalters@0:     {   "scale_fft" ,   "0.08"      ,  &smagstr     ,   "Scale factor for fft magnitude output"      , SVAL    },
tomwalters@0:     {   "scale_phs" ,   "100"       ,  &sphstr      ,   "Scale factor for fft phase output"          , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int     frameheight ;   /* Nap parameters read from header */
tomwalters@0: int     frames, samplerate ;
tomwalters@0: int     framewidth  ;   /* width option */
tomwalters@0: 
tomwalters@0: int     start, length ;
tomwalters@0: 
tomwalters@0: int     Newframeheight, Newframewidth ; /* Sai parameters to write  */
tomwalters@0: int     Newframes ;
tomwalters@0: int     Newframeshift ;
tomwalters@0: int     Newframebytes ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: char   *limitstr ;      /* Sai parameters used to get centre-frequencies */
tomwalters@0: char   *qualstr  ;
tomwalters@0: char   *minstr   ;
tomwalters@0: char   *maxstr   ;
tomwalters@0: char   *denstr   ;
tomwalters@0: char   *chansstr ;
tomwalters@0: double *frequencies ;
tomwalters@0: int     nerbs ;         /* bandwidth defined as number of erbs either side of cf */
tomwalters@0: int     bandlimit ;     /* flag for zeroing outside bandwidth around centre-freqs */
tomwalters@0: int     align ;         /* flag for aligning centre-freqs */
tomwalters@0: 
tomwalters@0: #define NONE    0       /* forms of auditory image */
tomwalters@0: #define ACGRAM  1
tomwalters@0: #define FTGRAM  2
tomwalters@0: #define PHGRAM  3
tomwalters@0: int     image ;
tomwalters@0: 
tomwalters@0: float  *buf ;           /* fft working space */
tomwalters@0: int     zeroes ;        /* padding */
tomwalters@0: float  *W ;             /* Hamming window for fft */
tomwalters@0: short   window ;        /* flag for window */
tomwalters@0: 
tomwalters@0: short   debug ;         /* flag for debug printf's */
tomwalters@0: 
tomwalters@0: main(argc, argv)
tomwalters@0: int   argc ;
tomwalters@0: char *argv[] ;
tomwalters@0: {
tomwalters@0:     FILE   *fp ;
tomwalters@0:     int     i, framebytes, startbytes, frameshift, framebyteshift;
tomwalters@0:     char   *header, *SaiHeader();
tomwalters@0:     short  *nap, *frame, *endframe;
tomwalters@0:     int     a, b ;
tomwalters@0:     char   *val1, *val2 ;
tomwalters@0: 
tomwalters@0:     fp = openopts( option,argc,argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts( helpstr, argv[0], applic, option ) ;
tomwalters@0: 
tomwalters@0:     window = ison( wstr ) ;
tomwalters@0:     nerbs  = atoi( erbstr ) ;
tomwalters@0:     debug  = ison( debugstr ) ;
tomwalters@0: 
tomwalters@0:     /**** Form of image ****/
tomwalters@0: 
tomwalters@0:     if ( iststr( imstr, "acgram" ) )  {
tomwalters@0: 	image = ACGRAM ;
tomwalters@0:     }
tomwalters@0:     else if ( iststr( imstr, "ftgram" ) )  {
tomwalters@0: 	image = FTGRAM ;
tomwalters@0:     }
tomwalters@0:     else if ( iststr( imstr, "phgram" ) )  {
tomwalters@0: 	image = PHGRAM ;
tomwalters@0:     }
tomwalters@0:     else if ( isoff( imstr ) ) {
tomwalters@0: 	image = NONE ;
tomwalters@0:     }
tomwalters@0:     else {
tomwalters@0: 	fprintf(stderr, "audim: unknown form of auditory image [%s]\n", imstr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     align     = ison( cfstr ) ;
tomwalters@0:     bandlimit = ison( blstr ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** parse bounds on number of frames ****/
tomwalters@0: 
tomwalters@0:     if ( getvals( framestr, &val1, &val2, "-" ) == BADVAL ) {
tomwalters@0: 	fprintf(stderr,"audim: bad frame selector [%s]\n", framestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     a = atoi( val1 ) ;
tomwalters@0:     if ( isempty( val2 ) )    b = a ;
tomwalters@0:     else if ( ismax( val2 ) ) b = 0 ;
tomwalters@0:     else                      b = atoi( val2 ) ;
tomwalters@0: 
tomwalters@0:     if (b<a && b>0) fprintf(stderr,"audim warning: bad frame specifiers\n");
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Read header and options to find dimensions of auditory image ****/
tomwalters@0: 
tomwalters@0:     if ( (header = ReadHeader(fp)) == (char *) 0 ) {
tomwalters@0: 	fprintf(stderr,"audim: header not found\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     frameheight = HeaderInt( header, "frameheight" );
tomwalters@0:     frames      = HeaderInt( header, "frames"      );
tomwalters@0:     samplerate  = HeaderSamplerate( header );
tomwalters@0:     start       = to_p( startstr,  samplerate );
tomwalters@0:     frameshift  = to_p( shiftstr, samplerate ) ;
tomwalters@0: 
tomwalters@0:     if ( ismax( lengthstr ) )   /* Special case for remainder of input */
tomwalters@0: 	length = frames - start ;
tomwalters@0:     else
tomwalters@0: 	length = to_p( lengthstr, samplerate );
tomwalters@0: 
tomwalters@0:     if ( start + length > frames ) {
tomwalters@0: 	fprintf(stderr,"audim: nap too small (%d ms) for requested start and length \n", to_ms( frames, samplerate ) );
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     framebytes = frameheight * length * sizeof(short) ;
tomwalters@0:     startbytes = frameheight * start  * sizeof(short) ;
tomwalters@0: 
tomwalters@0:     if ( ismax( widthstr ) )    /* Special case for single frame of max width */
tomwalters@0: 	framewidth = length ;
tomwalters@0:     else
tomwalters@0: 	framewidth = to_p( widthstr, samplerate ) ;
tomwalters@0: 
tomwalters@0:     frames = ( length - framewidth ) / frameshift ;
tomwalters@0:     zeroes = getpower( framewidth << 1 ) - framewidth ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Get centre-frequencies information using info given in header ****/
tomwalters@0: 
tomwalters@0:     limitstr =    HeaderStringOnly( header, "bwmin_afb"    );
tomwalters@0:     qualstr =     HeaderStringOnly( header, "quality_afb"  );
tomwalters@0:     minstr =      HeaderStringOnly( header, "mincf_afb"    );
tomwalters@0:     maxstr =      HeaderStringOnly( header, "maxcf_afb"    );
tomwalters@0:     denstr =      HeaderStringOnly( header, "dencf_afb"    );
tomwalters@0:     chansstr =    HeaderStringOnly( header, "channels_afb" );
tomwalters@0: 
tomwalters@0:     SetErbParameters( to_Hz( limitstr ), atof( qualstr ) ) ;
tomwalters@0:     if( OptionInt( chansstr ) == 0 )
tomwalters@0: 	frequencies = GenerateCenterFrequencies( to_Hz( minstr ), to_Hz( maxstr ), atof( denstr ) )   ;
tomwalters@0:     else
tomwalters@0: 	frequencies = NumberedCenterFrequencies( to_Hz( minstr ), to_Hz( maxstr ), OptionInt( chansstr ) ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Output frame dimensions (to write into new header) ****/
tomwalters@0: 
tomwalters@0:     Newframeheight = frameheight ;
tomwalters@0: 
tomwalters@0:     if ( !align ) {
tomwalters@0: 	if ( image == NONE )
tomwalters@0: 	    Newframewidth = framewidth ;
tomwalters@0: 	else
tomwalters@0: 	    Newframewidth = ( framewidth + zeroes ) >> 1 ;  /* spectral bins */
tomwalters@0:     }
tomwalters@0:     else  {      /* number of bins corresponding to the filter bandwidth */
tomwalters@0: 	switch ( image ) {
tomwalters@0: 	    case NONE   :   Newframewidth = framewidth ;
tomwalters@0: 			    break ;
tomwalters@0: 	    case ACGRAM :   Newframewidth = maxlagwidth( to_Hz( minstr ), nerbs ) ;
tomwalters@0: 			    break ;
tomwalters@0: 	    case FTGRAM :   Newframewidth = maxbandwidth( to_Hz( maxstr ), nerbs, (framewidth+zeroes)>>1 ) ;
tomwalters@0: 			    break ;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0:     Newframeshift =  frameshift ;
tomwalters@0:     Newframebytes =  Newframeheight * Newframewidth * sizeof(short) ;
tomwalters@0: 
tomwalters@0:     if ( frameshift > 0 ) {
tomwalters@0: 	if ( b == 0 )
tomwalters@0: 	    Newframes = frames - (a-1) ;
tomwalters@0: 	else
tomwalters@0: 	    Newframes = b - (a-1) ;
tomwalters@0:     }
tomwalters@0:     else {
tomwalters@0: 	fprintf(stderr,"audim: non-positive frstep [%d]\n", frameshift);
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Check limits etc.. ****/
tomwalters@0: 
tomwalters@0:     if ( b > frames ) {
tomwalters@0: 	fprintf(stderr,"audim: can't select frame %d out of %d frames\n", b, frames ) ;
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( length < framewidth ) {
tomwalters@0: 	fprintf(stderr,"audim: nap too small (%d ms) for requested width\n", to_ms( length, samplerate ) );
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if (frames<=0) {
tomwalters@0: 	if (frames==0)  fprintf(stderr,"audim: zero frames input\n");
tomwalters@0: 	if (frames<0)   fprintf(stderr,"audim: garbled number of frames, (start set too big?)\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if (a>frames || b>frames)
tomwalters@0: 	fprintf(stderr,"audim warning: bad frame selectors\n");
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( debug ) {
tomwalters@0: 	fprintf(stderr, "Header:  frames=%d frameheight=%d\n", frames, frameheight ) ;
tomwalters@0: 	fprintf(stderr, "Options: a=%d b=%d start=%d length=%d frameshift=%d framewidth=%d\n", a, b, start, length, frameshift, framewidth ) ;
tomwalters@0: 	fprintf(stderr, "Output:  zeroes=%d Newframewidth=%d Newframes=%d \n", zeroes, Newframewidth, Newframes ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Write new sai header ****/
tomwalters@0: 
tomwalters@0:     if ( ison( headerstr ) )
tomwalters@0: 	WriteHeader( SaiHeader(header), stdout );
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Allocate space for framebytes of nap data ****/
tomwalters@0: 
tomwalters@0:     if ( (nap = (short *)malloc( framebytes )) == NULL ) {
tomwalters@0: 	fprintf(stderr,"audim: malloc out of space\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Allocate and assign window coeffs ****/
tomwalters@0: 
tomwalters@0:     if ( window )
tomwalters@0: 	W = hamming( framewidth ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Allocate fft working space ****/
tomwalters@0: 
tomwalters@0:     buf = (float *)malloc( ( framewidth + zeroes ) * sizeof(float) );
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /**** Seek to start in blocks of framebytes ****/
tomwalters@0: 
tomwalters@0:     for (i=framebytes ; i < startbytes ; i += framebytes)
tomwalters@0: 	if ( fread( nap, framebytes, 1, fp ) == NULL ) {
tomwalters@0: 	    fprintf(stderr,"audim: missing data after header\n");
tomwalters@0: 	    exit(1);
tomwalters@0: 	}
tomwalters@0:     if ( (startbytes -= (i - framebytes)) > 0 )
tomwalters@0: 	if ( fread( nap, startbytes, 1, fp ) == NULL ) {
tomwalters@0: 	    fprintf(stderr,"audim: missing data after header\n");
tomwalters@0: 	    exit(1);
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     /**** Read framebytes of i/p nap data ****/
tomwalters@0: 
tomwalters@0:     if ( fread( nap, framebytes, 1, fp ) == NULL ) {
tomwalters@0: 	fprintf(stderr,"audim: missing data after header\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     framebyteshift = frameshift * frameheight ;
tomwalters@0:     nap = nap + (a-1) * framebyteshift ;
tomwalters@0:     endframe = nap + Newframes * framebyteshift ;
tomwalters@0: 
tomwalters@0:     switch ( image ) {
tomwalters@0: 	case NONE   : fprintf(stderr,"Output framewidth = %dp \n", Newframewidth ) ;  break ;
tomwalters@0: 	case ACGRAM : fprintf(stderr,"Output framewidth = %dp Lag range = 0-%dms\n", Newframewidth, (Newframewidth*1000)/samplerate ) ; break ;
tomwalters@0: 	case FTGRAM : fprintf(stderr,"Output framewidth = %dp Frequency range = 0-%dkHz\n", Newframewidth, samplerate/2000 ) ; break ;
tomwalters@0: 	case PHGRAM : break ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for ( frame = nap ; frame < endframe ; frame += framebyteshift ) {
tomwalters@0: 
tomwalters@0: 	fprintf(stderr,"audim: %d / %d\n", a++, frames );
tomwalters@0: 	switch ( image ) {
tomwalters@0: 
tomwalters@0: 	    case NONE   : writeframe( frame ) ; break ;
tomwalters@0: 	    case ACGRAM : acgram( frame ) ;     break ;
tomwalters@0: 	    case FTGRAM : ftgram( frame ) ;     break ;
tomwalters@0: 	    case PHGRAM : phgram( frame ) ;     break ;
tomwalters@0: 
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0:     fprintf(stderr,"audim: done\n" );
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int OptionInt( str )
tomwalters@0: char *str ;
tomwalters@0: {
tomwalters@0:     if( strcmp( str, "on" ) == 0 )
tomwalters@0: 	return( 1 ) ;
tomwalters@0:     else if( strcmp( str, "Not_used" ) == 0 )
tomwalters@0: 	return( 0 ) ;
tomwalters@0:     else
tomwalters@0: 	return( atoi( str ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*********************** Read header and build new header *****************/
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:    Copy the original nap header to a new sai header, changing in order:
tomwalters@0:      frames
tomwalters@0:      frameshift
tomwalters@0:      framewidth
tomwalters@0:      frameheight
tomwalters@0:      framebytes
tomwalters@0:      animate
tomwalters@0:      display
tomwalters@0:    Finally, update the new header_bytes, and return the new header.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: char *SaiHeader( napheader )
tomwalters@0: char *napheader ;
tomwalters@0: {
tomwalters@0:     char *saiheader;
tomwalters@0:     char *p0, *p1, *p2, *s, str[64];
tomwalters@0: 
tomwalters@0:     saiheader = (char *)malloc( strlen(napheader) + 256 ) ;
tomwalters@0: 
tomwalters@0:     p0 = saiheader ;
tomwalters@0:     p1 = napheader ;
tomwalters@0: 
tomwalters@0:     /* copy to first item after the header_bytes */
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "header_bytes" ) ;
tomwalters@0:     while( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0:     *p0++ = *p1++;
tomwalters@0: 
tomwalters@0:     /* insert frstep_aid at top of new header */
tomwalters@0: 
tomwalters@0:     sprintf( str,"frstep_aid=%s\n", shiftstr ) ;
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0: 
tomwalters@0:     /** copy up to frames **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "frames" ) ;
tomwalters@0:     while( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%d\n", Newframes);
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy up to frameshift **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "frameshift" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%d\n", Newframeshift);
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy up to framewidth **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "framewidth" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%d\n", Newframewidth);
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy up to frameheight **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "frameheight" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%d\n", Newframeheight);
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy up to framebytes **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "framebytes" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%d\n", Newframebytes);
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy up to animate_ctn **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "animate_ctn" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%s\n", anistr );
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy up to display **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "display" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     sprintf(str,"%s\n", dispstr );
tomwalters@0:     for (s = str ; *s != '\n' ; )
tomwalters@0: 	*p0++ = *s++;
tomwalters@0:     *p0++ = *s;
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p1++;
tomwalters@0:     *p1++;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** copy rest of header **/
tomwalters@0: 
tomwalters@0:     p2 = HeaderString( napheader , "Version" ) ;
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     if ( ( p2 = ApplicString( napheader ) ) == (char *)0 )  {
tomwalters@0: 	fprintf(stderr,"audim: application name not found in header\n");
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     while ( p1 < p2 )
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0:     *p0++ = 's' ; *p0++ = 'a' ; *p0++ = 'i' ;   /* copy sai into applic name */
tomwalters@0:     p1 += 3 ;
tomwalters@0: 
tomwalters@0:     while (*p1 != '\n')
tomwalters@0: 	*p0++ = *p1++ ;
tomwalters@0:     *p0++ = *p1++ ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /** update header_bytes **/
tomwalters@0: 
tomwalters@0:     sprintf(str, "%0*d", 7, p0-saiheader);
tomwalters@0:     p0 = HeaderString( saiheader , "header_bytes" ) ;
tomwalters@0:     s = str;
tomwalters@0:     while(*p0 != '\n')
tomwalters@0: 	*p0++ = *s++ ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     return saiheader;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /************************** Autocorrelation function ***********************/
tomwalters@0: 
tomwalters@0: /* Call acf for each row in auditory image */
tomwalters@0: 
tomwalters@0: acgram( frame )
tomwalters@0: short *frame ;
tomwalters@0: {
tomwalters@0:     static   float  scale ;
tomwalters@0:     static   int    normalize ;
tomwalters@0:     static   int    first = 1 ;
tomwalters@0:     register int    i, j, row, col ;
tomwalters@0:     short    p ;
tomwalters@0: 
tomwalters@0:     if (first) {
tomwalters@0: 	normalize = ison( normstr ) ;
tomwalters@0: 	scale = 1./(double)( framewidth + zeroes ) * atof( sacfstr ) ;
tomwalters@0: 	first=0;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for ( row=0 ; row < frameheight ; row++ ) {
tomwalters@0: 	for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
tomwalters@0: 	    buf[col] = frame[i] ;
tomwalters@0: 	for ( j=0 ; j<zeroes ; j++ )    /* padding */
tomwalters@0: 	    buf[col++] = 0 ;
tomwalters@0: 
tomwalters@0: 	acf( buf, framewidth+zeroes ) ;
tomwalters@0: 
tomwalters@0: 	if ( align || bandlimit )
tomwalters@0: 	    writelags( buf, row, (framewidth+zeroes)>>1, Newframewidth, scale ) ;
tomwalters@0: 	else {
tomwalters@0: 	    if ( normalize == 0 )
tomwalters@0: 		writebuf( buf, Newframewidth, scale ) ;
tomwalters@0: 	    else
tomwalters@0: 		writebuf( buf, Newframewidth, 5000./buf[0] ) ;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /********************************** FFT ************************************/
tomwalters@0: 
tomwalters@0: /* Call fft and magnitude spectrum for each row in auditory image */
tomwalters@0: 
tomwalters@0: ftgram( frame )
tomwalters@0: short *frame ;
tomwalters@0: {
tomwalters@0:     static   float  scale ;
tomwalters@0:     static   int    first = 1 ;
tomwalters@0:     register int    i, j, row, col ;
tomwalters@0: 
tomwalters@0:     if (first) {
tomwalters@0: 	scale = atof( smagstr ) ;
tomwalters@0: 	first=0;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for ( row=0 ; row < frameheight ; row++ ) {
tomwalters@0: 	if ( window )
tomwalters@0: 	    for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
tomwalters@0: 		buf[col] = frame[i] * W[col] ;
tomwalters@0: 	else
tomwalters@0: 	    for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
tomwalters@0: 		buf[col] = frame[i] ;
tomwalters@0: 
tomwalters@0: 	for ( j=0 ; j<zeroes ; j++ )    /* padding */
tomwalters@0: 	    buf[col++] = 0 ;
tomwalters@0: 
tomwalters@0: 	fft( buf, framewidth+zeroes, 1 ) ;
tomwalters@0: 	mag( buf, framewidth+zeroes ) ;
tomwalters@0: 
tomwalters@0: 	if ( align || bandlimit )
tomwalters@0: 	    writebins( buf, row, (framewidth+zeroes)>>1, Newframewidth, scale ) ;
tomwalters@0: 	else
tomwalters@0: 	    writebuf( buf, Newframewidth, scale ) ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /********************************** Phase **********************************/
tomwalters@0: 
tomwalters@0: /* Call fft and phase spectrum for each row in auditory image */
tomwalters@0: 
tomwalters@0: phgram( frame )
tomwalters@0: short *frame ;
tomwalters@0: {
tomwalters@0:     static   float  scale ;
tomwalters@0:     static   int    first = 1 ;
tomwalters@0:     register int    i, j, row, col ;
tomwalters@0: 
tomwalters@0:     if (first) {
tomwalters@0: 	scale = atof( sphstr ) ;
tomwalters@0: 	first=0;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for ( row=0 ; row < frameheight ; row++ ) {
tomwalters@0: 	for ( col=0 , i=row ; col < framewidth ; col++, i+=frameheight )
tomwalters@0: 	    buf[col] = frame[i] ;
tomwalters@0: 	for ( j=0 ; j<zeroes ; j++ )    /* padding */
tomwalters@0: 	    buf[col++] = 0 ;
tomwalters@0: 
tomwalters@0: 	fft( buf, framewidth+zeroes, 1 ) ;
tomwalters@0: 	phase( buf, framewidth+zeroes ) ;
tomwalters@0: 
tomwalters@0: 	writebuf( buf, Newframewidth, scale ) ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /**** Write o/p routines ****/
tomwalters@0: 
tomwalters@0: writebuf( buf, n, scale )       /* write buffer as scaled shorts */
tomwalters@0: float *buf ;
tomwalters@0: int    n ;
tomwalters@0: float  scale ;
tomwalters@0: {
tomwalters@0:     register int  i ;
tomwalters@0:     short         p ;
tomwalters@0: 
tomwalters@0:     for (i=0 ; i < n ; i++)  {
tomwalters@0: 	p = (short)( buf[i] * scale ) ;
tomwalters@0: 	fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: writeframe( frame )     /* Write a nap frame in sai format */
tomwalters@0: short  *frame ;
tomwalters@0: {
tomwalters@0:     int    i, row, col;
tomwalters@0: 
tomwalters@0:     for (row=0 ; row < frameheight ; row++)
tomwalters@0: 	for (col=0 , i=row ; col < framewidth ; col++, i+=frameheight)
tomwalters@0: 	    fwrite( &frame[i], sizeof(short), 1, stdout );
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:   Write a band of spectral bins around the centre-freq of the given channel.
tomwalters@0: 
tomwalters@0:   The o/p spectrum is 0 to samplerate/2 Hz, resolved into spectral bins:
tomwalters@0: 	bins = (framewidth + zeroes) / 2
tomwalters@0: 
tomwalters@0:   For the i'th channel, the centre-frequency cf = frequencies[i].
tomwalters@0:   Filters are symmetrical bandpass, so considering a bandwidth of 4 ERB
tomwalters@0:   the side-frequencies at the band limits are:
tomwalters@0: 
tomwalters@0:     side1 = cf - 2*ERB     [Hz]
tomwalters@0:     side2 = cf + 2*ERB     [Hz]
tomwalters@0: 
tomwalters@0:   Note: side1 < side2
tomwalters@0: 
tomwalters@0:   Convert frequency to bin number as follows:
tomwalters@0:   Let centre-frequency cf have bin number j, so that ratios:
tomwalters@0: 	cf/(samplerate/2) = j/bins
tomwalters@0:   and so:
tomwalters@0: 	j  = ( cf*bins ) / ( samplerate/2 )
tomwalters@0: 
tomwalters@0:   (Similarly for each of the side-frequencies at the 2-ERB band limits).
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: writebins( buf, chan, bins, maxwidth, scale )
tomwalters@0: float *buf ;
tomwalters@0: int    chan ;
tomwalters@0: int    bins ;       /* number of spectral bins */
tomwalters@0: int    maxwidth ;
tomwalters@0: float  scale ;
tomwalters@0: {
tomwalters@0:     double ERB, cf, side1, side2, max ;
tomwalters@0:     register int  i, j, s1, s2 ;
tomwalters@0:     short         p ;
tomwalters@0: 
tomwalters@0:     ERB   = frequencies[chan]/9.26449 + 24.7 ;
tomwalters@0:     cf    = frequencies[chan]          ;
tomwalters@0:     side1 = frequencies[chan] - nerbs*ERB  ;
tomwalters@0:     side2 = frequencies[chan] + nerbs*ERB  ;
tomwalters@0: 
tomwalters@0:     j  = (int)( ( cf    * bins ) / ( samplerate >> 1 ) + 0.5 ) ;   /* rounding */
tomwalters@0:     s1 = (int)( ( side1 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;
tomwalters@0:     s2 = (int)( ( side2 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;
tomwalters@0: 
tomwalters@0:     if ( debug )
tomwalters@0: 	fprintf(stderr,"chan=%d bins=%d maxwidth=%d j=%d side1=%d side2=%d\n",chan,bins,maxwidth,j,s1,s2);
tomwalters@0: 
tomwalters@0:     /*  Plot the given bandwidth, zeroing elsewhere:     */
tomwalters@0: 
tomwalters@0:     if ( bandlimit && !align ) {
tomwalters@0: 
tomwalters@0: 	for ( i=0, p=0 ; i < s1 ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 
tomwalters@0: 	for (  ; i < s2 && i < bins ; i++ )  {
tomwalters@0: 	    p = (short)( buf[i] * scale ) ;
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	for ( p=0 ; i < bins ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* Plot the given bandwidth, shifting it to the centre */
tomwalters@0: 
tomwalters@0:     if ( align ) {
tomwalters@0: 
tomwalters@0: /* align on max bin in band
tomwalters@0: 	max = 0 ;
tomwalters@0: 	for ( i=s1 ; i<s2 ; i++ )
tomwalters@0: 	    if ( buf[i] >= max ) {
tomwalters@0: 		max = buf[i] ;
tomwalters@0: 		j = i ;
tomwalters@0: 	    }
tomwalters@0: 	for ( i=0, p=0 ; i < s1+(maxwidth/2-j) ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 	for ( j = s1 ; j < s2 && i < maxwidth ; i++, j++ )  {
tomwalters@0: 	    p = (short)( buf[j] * scale ) ;
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 	for ( p=0 ; i < maxwidth ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 	j = (maxwidth+s1-s2)/2 ;
tomwalters@0: 	for ( i=0, p=0 ; i < j ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 
tomwalters@0: 	for ( j = s1 ; j < s2 ; i++, j++ )  {
tomwalters@0: 	    p = (short)( buf[j] * scale ) ;
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 	for ( p=0 ; i < maxwidth ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:   Write a band of acf lags around the centre-period of the given channel.
tomwalters@0: 
tomwalters@0:   For the i'th channel, the centre-frequency
tomwalters@0: 	Centre-frequency   cf  = frequencies[i]     [Hz]
tomwalters@0: 	Bandwidth          ERB = cf/9.26449 + 24.7  [Hz]
tomwalters@0: 	Centre-period       cp = samplerate / cf    [samples]
tomwalters@0: 
tomwalters@0:   Filters are symmetrical bandpass, so considering a bandwidth of 4 ERB
tomwalters@0: 
tomwalters@0:     side1 = samplerate / ( cf - 2*ERB )    [samples]
tomwalters@0:     side2 = samplerate / ( cf + 2*ERB )    [samples]
tomwalters@0: 
tomwalters@0:   Note: side2 < side1 since periods
tomwalters@0: 
tomwalters@0:   Also, successive bins ARE successive lags (in samples).
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: writelags( buf, chan, lags, maxwidth, scale )
tomwalters@0: float *buf ;
tomwalters@0: int    chan ;
tomwalters@0: int    lags ;       /* number of acf lags */
tomwalters@0: int    maxwidth ;
tomwalters@0: float  scale ;
tomwalters@0: {
tomwalters@0:     double ERB, cp, side1, side2 ;
tomwalters@0:     register int  i, j, s1, s2 ;
tomwalters@0:     short         p ;
tomwalters@0: 
tomwalters@0:     ERB   = frequencies[chan]/9.26449 + 24.7 ;
tomwalters@0:     cp    = (double)samplerate /   frequencies[chan]           ;
tomwalters@0:     side1 = (double)samplerate / ( frequencies[chan] - nerbs*ERB ) ;
tomwalters@0:     side2 = (double)samplerate / ( frequencies[chan] + nerbs*ERB ) ;
tomwalters@0: 
tomwalters@0:     j  = (int)( cp    + 0.5 ) ;         /* rounding */
tomwalters@0:     s1 = (int)( side1 + 0.5 ) ;
tomwalters@0:     s2 = (int)( side2 + 0.5 ) ;
tomwalters@0: 
tomwalters@0:     if ( debug )
tomwalters@0: 	fprintf(stderr,"chan=%d lags=%d side1=%d side2=%d\n",chan,lags,s1,s2);
tomwalters@0: 
tomwalters@0:     /*  Plot the given bandwidth, zeroing elsewhere:     */
tomwalters@0: 
tomwalters@0:     if ( bandlimit && !align ) {
tomwalters@0: 
tomwalters@0: 	for ( i=0, p=0 ; i < s2 ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 
tomwalters@0: 	for (  ; i < s1 && i < lags ; i++ )  {
tomwalters@0: 	    p = (short)( buf[i] * scale ) ;
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	for ( p=0 ; i < lags ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* Plot the given bandwidth, shifting it to the centre */
tomwalters@0: 
tomwalters@0:     if ( align ) {
tomwalters@0: 
tomwalters@0: 	j = (maxwidth+s2-s1)/2 ;
tomwalters@0: 	for ( i=0, p=0 ; i < j ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 
tomwalters@0: 	for ( j = s2 ; j < s1 ; i++, j++ )  {
tomwalters@0: 	    p = (short)( buf[j] * scale ) ;
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 	for ( p=0 ; i < maxwidth ; i++ )
tomwalters@0: 	    fwrite( &p, sizeof(short), 1, stdout ) ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:   Return the number of lags corresponding to a bandwidth of nerbs ERBS either
tomwalters@0:   side of a given mincf. This will be the max required framewidth.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: maxlagwidth( mincf, nerbs )
tomwalters@0: double mincf ;
tomwalters@0: int    nerbs ;
tomwalters@0: {
tomwalters@0:     double ERB, side1, side2 ;
tomwalters@0: 
tomwalters@0:     ERB   = mincf/9.26449 + 24.7 ;
tomwalters@0:     side1 = (double)samplerate / ( mincf - nerbs*ERB ) ;
tomwalters@0:     side2 = (double)samplerate / ( mincf + nerbs*ERB ) ;
tomwalters@0: 
tomwalters@0:     return (int)( side1 - side2 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:   Return the number of bins corresponding to a bandwidth of nerbs ERBS either
tomwalters@0:   side of a given maxcf. This will be the max required framewidth.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: maxbandwidth( maxcf, nerbs, bins )
tomwalters@0: double maxcf ;
tomwalters@0: int    nerbs ;
tomwalters@0: int    bins  ;
tomwalters@0: {
tomwalters@0:     double ERB, side1, side2 ;
tomwalters@0:     int    s1, s2 ;
tomwalters@0: 
tomwalters@0:     ERB   = maxcf/9.26449 + 24.7 ;
tomwalters@0:     side1 = maxcf - nerbs*ERB  ;
tomwalters@0:     side2 = maxcf + nerbs*ERB  ;
tomwalters@0: 
tomwalters@0:     /* convert frequencies to bin numbers */
tomwalters@0: 
tomwalters@0:     s1 = (int)( ( side1 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;
tomwalters@0:     s2 = (int)( ( side2 * bins ) / ( samplerate >> 1 ) + 0.5 ) ;
tomwalters@0: 
tomwalters@0:     return (int)( s2 - s1 ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: