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