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

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

/*
 ftgram.c       Row-wise fft program by Michael Allerhand.
----------

  Short-time Fourier transform applied to each row of input frames.
  Output frames consist of row-wise magnitude spectra, computed using fft.

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

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

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

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

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

Padding:
  The output framewidth is the next power-of-2 above twice the input
  framewidth.
  The rows of each input frame are padded with zeroes so that the input
  framewidth is at least twice its original size, and also a power of 2.
  This is because the fft method is "radix-2", and because each output
  spectrum is symmetrical so that just half the points are unique.
  It is advisable to pad input frames for the additional reason that
  padding counters frame end-effects.

  Therefore each row of each input frame is padded with zeroes to the next
  power of 2 larger than the original input framewidth.
  If necessary, extra padding can be enforced using the padding option to
  add extra zeroes, padding to a larger power of 2.
  The amount of extra padding is "exponential" since if the basic padded
  framesize is N, and the padding option is n, then extra padding is to:
  N * 2**n.
  (n=0 by default, so that N is not changed. When n=1 then N is doubled, and
  when n=2 then N is quadrupled, etc.).

  The output framewidth (and hence the frequency resolution) depends upon the
  input framewidth and the padding.

Examples:

1. ftgram of a NAP, animated, with framewidth 16ms

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

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

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


*/

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

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

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

static Options option[] = {
    {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                       , DEBUG   },
    {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                           , DEBUG   },
    {   "display"   ,   "on"        ,  &dispstr     ,   "Display output image."                      , SETFLAG },
    {   "animate"   ,   "off"       ,  &anistr      ,   "Animate cartoon."                           , SETFLAG },
    {   "frames"    ,   "1-max"     ,  &framestr    ,   "Select frames inclusively"                  , VAL     },
    {   "header"    ,   "on"        ,  &headerstr   ,   "Header (for gensai useprevious=on)."        , VAL     },
    {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p file."                   , VAL     },
    {   "length"    ,   "max"       ,  &lengthstr   ,   "Length of i/p file to process."             , VAL     },
    {   "frstep"    ,   "16ms"      ,  &shiftstr    ,   "Step between input frames."                 , VAL     },
    {   "width"     ,   "32ms"      ,  &widthstr    ,   "Width of input frames."                     , VAL     },
    {   "padding"   ,   "0"         ,  &padstr      ,   "Exponential frame padding"                  , SVAL    },
    {   "window"    ,   "on"        ,  &wstr        ,   "Hamming window"                             , SETFLAG },
    {   "scale"     ,   "0.08"      ,  &scalestr    ,   "Scale factor for output"                    , SVAL    },
   ( char * ) 0 } ;


int     samplerate ;

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

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

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

int     zeroes ;
int     window ;        /* Flag for Hamming window */
float   scale  ;        /* Scale factor for output */

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


main(argc, argv)
int   argc ;
char *argv[] ;
{
    FILE   *fp ;
    char   *header, *SaiHeader();
    char   *versionstr, c ;
    int     startbytes, framebyteshift ;
    short  *frame, *endframe ;
    int     i, a, b ;

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

    window  = ison( wstr ) ;


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

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


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

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

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

    samplerate  = HeaderSamplerate( header );

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

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

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

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

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

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

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

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

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

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

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

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

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


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

    Newframeheight =  frameheight ;
    Newframewidth  =  ( framewidth + zeroes ) >> 1 ;  /* spectral bins */
    Newframeshift  =  frameshift ;
    Newframebytes  =  Newframeheight * Newframewidth * sizeof(short) ;


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

    zeroes = ( getpower( framewidth << 1 ) << atoi( padstr ) ) - framewidth ;
    scale  = atof( scalestr ) ;


    /**** Debug ****/

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


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

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

    if ( window )
	W = hamming( framewidth ) ;

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


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

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


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

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


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

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


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

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

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

    fclose(fp);
    fprintf(stderr,"ftgram: done\n" ) ;
}



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

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

ftgram( frame )
short *frame ;
{
    register int  i, j, row, col ;

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

	writebuf( buf, Newframewidth, scale ) ;
    }
}


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

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


/************************ Output 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) + 64 ) ;

    p0 = saiheader ;
    p1 = napheader ;


    /** copy up to frames **/

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

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


    /** copy up to frameshift **/

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

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


    /** copy up to framewidth **/

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

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


    /** copy up to frameheight **/

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

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


    /** copy up to framebytes **/

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

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


    /** copy up to animate_ctn **/

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

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


    /** copy up to display **/

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

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


    /** copy rest of header **/

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

    while (*p1 != '\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;
}