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

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

/*
  FFT program adapted from routines given in Numerical Recipes.

  The output framewidth = ( input framewidth + padding ) / 2  spectral points.

  Padding with zeroes is recommended to counter end-effects (See NR, p416).
  It is also necessary since framewidth + padding must be a power of 2
  (sample points at given rate).
  This is because the fft method is "radix-2", and because the output is
  symmetrical so that just half the points are unique.
  Therefore each input frame is padded with zeroes to the next power of 2
  larger than the 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", expanding the basic size to:
	( framewidth + padding ) * 2**n
  where the padding option is n.
  (n=0 by default, so that no extra padding is added. When n=1 then padding is
  added to double the size, and when n=2 the size is quadrupled, etc.).

  Frames are selected from the input stream using the "frames" option,
  which has syntax: [[-]frame=a[-b]]. Input frames  are numbered 1,2,...
  For example:
    frame=a      Select just the a'th frame.
    frame=a-b    Select frames from the a'th to b'th inclusive.
  The strings "min" and "max" can be used as specifiers, meaning eg:
    frame=min    Select the first frame.
    frame=max    Select the last frame.
    frame=a-max  Select frames from the a'th to the last inclusive.
    frame=min-b  Select frames from the first to the b'th inclusive.
  The default selects all frames, (ie frame=min-max).

  Several forms of fft processing are provided, according to the "spectrum"
  option:  log/magnitude/phase/complex/inverse/verbose.

  Log and magnitude are both magnitude spectra (log is log10 of magnitude).
  Phase is the phase spectrum.
  Complex is the full complex spectrum, in <real,imag> pairs.
  Inverse transform reads framewidth numbers which are interpreted as
  <real,imag> pairs, (ie framewidth/2 complex numbers),
  and outputs the inverse transform scaled by 1/framewidth.
  Verbose prints the spectrum in ASCII on the stdout.


  Examples:

1. To print the input and output frame sizes in sample points, eg for a
   subsequent plotting program, use the size option:

fft ... size=on

2. An fft of a waveform sampled at 10kHz, computed within a frame of 12.8ms,
   plotting the 2nd frame in a sequence of frames with half-frame overlap.

fft samp=10kHz width=12.8ms frstep=6.4ms frame=2 file | x11plot

3. An animated plot of successive fft spectra of a waveform sampled at 10kHz,
   each computed within a frame of 12.8ms, and shifted by 2 sample points.

fft samp=10kHz width=12.8ms frstep=2p file | x11play -n64

4. Using the complex output from fft, and inverse transform without windowing
   to recover original input.

fft samp=10kHz frame=2 spec=complex window=off file > foo
fft samp=10kHz frame=1 spec=inverse window=off foo | x11plot

*/

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

char applic[]     = "fast Fourier transform of contiguous frames.\n      i/p and o/p data in binary shorts." ;

static char *helpstr, *debugstr, *sampstr, *widthstr, *padstr, *wstr, *sstr, *sizestr ;
static char *startstr, *shiftstr, *framestr, *smagstr, *slogstr, *sphasestr, *echostr ;

static Options option[] = {
    {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                          , DEBUG   },
    {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                              , DEBUG   },
    {   "samplerate",   "20kHz"     ,  &sampstr     ,   "samplerate "                                   , VAL     },
    {   "start"     ,   "0"         ,  &startstr    ,   "Start point in i/p file."                      , VAL     },
    {   "frames"    ,   "1-max"     ,  &framestr    ,   "Select frames inclusively"                     , 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 },
    {   "spectrum"  ,   "magnitude" ,  &sstr        ,   "log/magnitude/phase/complex/inverse/verbose"   , SETFLAG },
    {   "scalemag"  ,   "0.08"      ,  &smagstr     ,   "scale magnitude  spectrum"                     , SVAL    },
    {   "scalelog"  ,   "10"        ,  &slogstr     ,   "scale log  spectrum"                           , SVAL    },
    {   "scalephase",   "100"       ,  &sphasestr   ,   "scale phase spectrum"                          , SVAL    },
    {   "size"      ,   "off"       ,  &sizestr     ,   "print input/output frame size in samples"      , SETFLAG },
    {   "echo"      ,   "off"       ,  &echostr     ,   "echo buffered input without processing"        , SVAL    },
   ( char * ) 0 } ;


int     samplerate  ;
int     width       ;
int     step        ;
int     zeroes      ;
int     dowindow    ;
double  scale       ;
int     echo        ;

float  *vec, *window ;
short  *obuf ;

int     allbins  ;
int     halfbins ;

#define MAG     1       /* magnitude spectrum */
#define PHASE   2       /* phase spectrum     */
#define INV     3       /* inverse transform  */
#define VERB    4       /* verbose (ASCII)    */
#define LOG     5       /* log magnitude      */
#define COMP    6       /* complex spectrum   */

int     output ;

main (argc, argv)
int     argc;
char  **argv;
{
    FILE     *fp  ;
    short    *buf ;
    int       a, b ;
    int       isign ;

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

    samplerate = to_Hz( sampstr ) ;
    dowindow   = ison( wstr ) ;
    echo       = ison( echostr ) ;
    width      = to_p( widthstr, samplerate )  ;
    step       = to_p( shiftstr, samplerate )  ;
    zeroes     = ( getpower( width ) << atoi( padstr ) ) - width ;

    allbins    = width + zeroes ;
    halfbins   = allbins / 2    ;

    if ( iststr( sstr, "magnitude" ) )  {
	output = MAG ;
	scale = atof( smagstr ) ;
	isign = 1 ;
    }
    else if ( iststr( sstr, "log" ) )  {
	output = LOG ;
	scale = atof( slogstr ) ;
	isign = 1 ;
    }
    else if ( iststr( sstr, "phase" ) ) {
	output = PHASE ;
	scale = atof( sphasestr ) ;
	isign = 1 ;
    }
    else if ( iststr( sstr, "complex" ) ) {
	output = COMP ;
	isign = 1 ;
    }
    else if ( iststr( sstr, "inverse" ) ) {
	output = INV ;
	scale = 1./(double)width ;
	isign = (-1) ;
    }
    else if ( iststr( sstr, "verbose" ) ) {
	output = VERB ;
	isign = 1 ;
    }
    else {
	fprintf(stderr,"unknown spectrum \n");
	exit( 1 ) ;
    }

    /* frame size printout */

    if ( ison( sizestr ) ) {
	fprintf(stderr,"fft sizes in sample points:\n" ) ;
	fprintf(stderr,"    input  frame size = %d  (framewidth=%d + padding=%d)\n", width+zeroes, width, zeroes ) ;
	fprintf(stderr,"    output frame size = %d \n", (width + zeroes)/2 ) ;
	exit( 0 ) ;
    }

    /* parse bounds on number of frames */

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

    /* Allocate working space */

    obuf  = (short *)malloc( halfbins * sizeof(short) ) ;
    vec   = (float *)malloc( allbins * sizeof(float) ) ;
    if ( dowindow )
	window = hamming( width ) ;


    /* Compute fft for each frame of width shorts in the input stream */

    if ( seekbytes( fp, (int)( to_p( startstr, samplerate ) * sizeof(short) ) ) == 0 ) {
	fprintf(stderr,"improper seek\n") ;
	exit( 1 ) ;
    }

    while ( ( buf = getframe( fp, a, width, step ) ) != (short *)0  && ( a<=b || b==0 ) ) {
	if ( echo ) fwrite( buf, sizeof(short), width, stdout ) ;
	else process( buf, isign ) ;
	a++ ;
    }
    if ( a<=b && b>0 )
	fprintf(stderr,"warning: not enough frames for request\n");

    fclose( fp ) ;
}


process( buf, isign )
short *buf ;
int    isign ;
{
    short *bptr     = buf ;
    short *endptr   = buf + width ;
    short *optr     = obuf ;
    short *endbin   = obuf + halfbins ;
    float *wptr     = window ;
    float *vptr     = vec    ;
    float *endvec   = vec + allbins ;

    if ( dowindow )
	while ( bptr < endptr )
	    *vptr++ =  *bptr++  *  *wptr++ ;
    else
	while ( bptr < endptr )
	    *vptr++ =  *bptr++ ;

    while ( vptr < endvec )         /* padding */
	*vptr++ =  0 ;

    fft( vec, width+zeroes, isign ) ;

    vptr = vec ;
    switch ( output ) {

	case COMP  :  while ( vptr < endvec )
			  *optr++ = (short)( *vptr++ ) ;
		      fwrite( obuf, sizeof(short), allbins, stdout ) ;
		      break;
	case PHASE :  phase( vec, allbins ) ;
		      while ( optr < endbin )
			  *optr++ = (short)( scale * *vptr++ ) ;
		      fwrite( obuf, sizeof(short), halfbins, stdout ) ;
		      break;
	case MAG   :  mag( vec, allbins ) ;
		      while ( optr < endbin )
			  *optr++ = (short)( scale * *vptr++ ) ;
		      fwrite( obuf, sizeof(short), halfbins, stdout ) ;
		      break;
	case LOG   :  mag( vec, allbins ) ;
		      while ( optr < endbin )
			  *optr++ = (short)( scale * 10 * log10(*vptr++) ) ;
		      fwrite( obuf, sizeof(short), halfbins, stdout ) ;
		      break;
	case INV   :  while ( vptr < endvec )
			  *optr++ = (short)( scale * *vptr++ ) ;
		      fwrite( obuf, sizeof(short), allbins, stdout ) ;
		      break;
	case VERB  :  print_complex_spectrum( (complex *)vec, halfbins, samplerate ) ;
		      break;

    }
}


/*
 Pretty-print complex k-spectrum on the stdout.
*/

print_complex_spectrum(C,k,samplerate)
complex *C;
int      k, samplerate;
{
    float res;      /* frequency resolution. */
    int   i;

    res = (float)(samplerate>>1) / (float)k;
    for (i=0; i<k ; i++) {
	printf ("%3d.  %6.1fHz.  Re=%9.2f  Im=%9.2f  mod=%8.2f  arg=%4.2f\n",
		i, i*res, Re(&C[i]), Im(&C[i]), mod(&C[i]), arg(&C[i]) );
    }
}