fft.c

Fast Fourier Transform in C

/* file: fft.c		G. Moody	24 February 1988
		   Last revised:	  10 June 2005

-------------------------------------------------------------------------------
fft: Fast Fourier transform of real data
Copyright (C) 1988-2005 George B. Moody

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.

You may contact the author by e-mail ([e-mail]) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA).  For updates to this software,
please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________
   
The default input to this program is a text file containing a uniformly sampled
time series, presented as a column of numbers.  The default standard output is
the magnitude of the discrete Fourier transform of the input series, normalized
by the length of the FFT.

The functions `realft' and `four1' are based on those in Press, W.H., et al.,
Numerical Recipes in C: the Art of Scientific Computing (Cambridge Univ. Press,
1989;  2nd ed., 1992).
*/

#include <stdio.h>
#include <math.h>
#ifdef __STDC__
#include <stdlib.h>
#else
extern void exit();
#endif

#ifndef BSD
# include <string.h>
#else           /* for Berkeley UNIX only */
# include <strings.h>
# define strchr index
#endif

#define	LEN	16384		/* maximum points in FFT */
#ifdef i386
#define strcasecmp strcmp
#endif
#define PI	M_PI	/* pi to machine precision, defined in math.h */
#define TWOPI	(2.0*PI)

void four1();
void realft();

double wsum;

/* See Oppenheim & Schafer, Digital Signal Processing, p. 241 (1st ed.) */
double win_bartlett(j, n)
int j, n;
{
    double a = 2.0/(n-1), w;

    if ((w = j*a) > 1.0) w = 2.0 - w;
    wsum += w;
    return (w);
}

/* See Oppenheim & Schafer, Digital Signal Processing, p. 242 (1st ed.) */
double win_blackman(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.42 - 0.5*cos(a*j) + 0.08*cos(2*a*j);
    wsum += w;
    return (w);
}

/* See Harris, F.J., "On the use of windows for harmonic analysis with the
   discrete Fourier transform", Proc. IEEE, Jan. 1978 */
double win_blackman_harris(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.35875 - 0.48829*cos(a*j) + 0.14128*cos(2*a*j) - 0.01168*cos(3*a*j);
    wsum += w;
    return (w);
}

/* See Oppenheim & Schafer, Digital Signal Processing, p. 242 (1st ed.) */
double win_hamming(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.54 - 0.46*cos(a*j);
    wsum += w;
    return (w);
}

/* See Oppenheim & Schafer, Digital Signal Processing, p. 242 (1st ed.)
   The second edition of Numerical Recipes calls this the "Hann" window. */
double win_hanning(j, n)
int j, n;
{
    double a = 2.0*PI/(n-1), w;

    w = 0.5 - 0.5*cos(a*j);
    wsum += w;
    return (w);
}

/* See Press, Flannery, Teukolsky, & Vetterling, Numerical Recipes in C,
   p. 442 (1st ed.) */
double win_parzen(j, n)
int j, n;
{
    double a = (n-1)/2.0, w;

    if ((w = (j-a)/(a+1)) > 0.0) w = 1 - w;
    else w = 1 + w;
    wsum += w;
    return (w);
}

/* See any of the above references. */
double win_square(j, n)
int j, n;
{
    wsum += 1.0;
    return (1.0);
}

/* See Press, Flannery, Teukolsky, & Vetterling, Numerical Recipes in C,
   p. 442 (1st ed.) or p. 554 (2nd ed.) */
double win_welch(j, n)
int j, n;
{
    double a = (n-1)/2.0, w;

    w = (j-a)/(a+1);
    w = 1 - w*w;
    wsum += w;
    return (w);
}

char *pname;
double rsum;
FILE *ifile;
int m, n;
int cflag;
int decimation = 1;
int iflag;
int fflag;
int len = LEN;
int nflag;
int Nflag;
int nchunks;
int pflag;
int Pflag;
int smooth = 1;
int wflag;
int zflag;
static float *c;
double freq, fstep, norm, rmean, (*window)();

main(argc, argv)
int argc;
char *argv[];
{
    int i;
    char *prog_name();
    double atof();
    void help();

    pname = prog_name(argv[0]);
    if (--argc < 1) {
	help();
	exit(1);
    }
    else if (strcmp(argv[argc], "-") == 0)
	ifile = stdin;		/* read data from standard input */
    else if ((ifile = fopen(argv[argc], "rt")) == NULL) {
	fprintf(stderr, "%s: can't open %s\n", pname, argv[argc]);
	exit(2);
    }
    for (i = 1; i < argc; i++) {
	if (*argv[i] == '-') switch (*(argv[i]+1)) {
	  case 'c':	/* print complex FFT (rectangular form) */
	    cflag = 1;
	    break;
	  case 'f':	/* print frequencies */
	    if (++i >= argc) {
		fprintf(stderr, "%s: sampling frequency must follow -f\n",
			pname);
		exit(1);
	    }
	    freq = atof(argv[i]);
	    fflag = 1;
	    break;
	  case 'h':	/* print help and exit */
	    help();
	    exit(0);
	    break;
	  case 'i':	/* calculate inverse FFT from `fft -c' format input */
	    if (argc > 2) {
		fprintf(stderr, "%s: no other option may be used with -i\n",
			pname);
		exit(1);
	    }
	    iflag = 1;
	    break;
	  case 'I':	/* calculate inverse FFT from `fft -p' format input */
	    if (argc > 2) {
		fprintf(stderr, "%s: no other option may be used with -I\n",
			pname);
		exit(1);
	    }
	    iflag = -1;
	    break;
	  case 'l':	/* perform up to n-point transforms */
	    if (++i >= argc) {
		fprintf(stderr, "%s: transform size must follow -l\n", pname);
		exit(1);
	    }
	    len = atoi(argv[i]);
	    break;
	  case 'n':	/* process in overlapping n-point chunks, output avg */
	    if (++i >= argc) {
		fprintf(stderr, "%s: chunk size must follow -n\n", pname);
		exit(1);
	    }
	    nflag = atoi(argv[i]);
	    break;
	  case 'N':	/* process in overlapping n-point chunks, output raw */
	    if (++i >= argc) {
		fprintf(stderr, "%s: chunk size must follow -N\n", pname);
		exit(1);
	    }
	    Nflag = atoi(argv[i]);
	    break;
	  case 'p':	/* print phases (polar form) */
	    pflag = 1;
	    break;
	  case 'P':	/* print power spectrum (squared magnitudes) */
	    Pflag = 1;
	    break;
	  case 's':	/* smooth spectrum */
	    if (++i >= argc || ((smooth=atoi(argv[i])) < 2) || smooth > 1024) {
		fprintf(stderr, "%s: smoothing parameter must follow -s\n",
			pname);
		exit(1);
	    }
	    break;
	  case 'S':	/* smooth and decimate spectrum */
	    if (++i >= argc || ((smooth=atoi(argv[i])) < 2) || smooth > 1024) {
		fprintf(stderr, "%s: smoothing parameter must follow -s\n",
			pname);
		exit(1);
	    }
	    decimation = smooth;
	    break;
	  case 'w':	/* apply windowing function to input */
	    if (++i >= argc) {
		fprintf(stderr, "%s: window type must follow -w\n",
			pname);
		exit(1);
	    }
	    if (strcasecmp(argv[i], "Bartlett") == 0)
		window = win_bartlett;
	    else if (strcasecmp(argv[i], "Blackman") == 0)
		window = win_blackman;
	    else if (strcasecmp(argv[i], "Blackman-Harris") == 0)
		window = win_blackman_harris;
	    else if (strcasecmp(argv[i], "Hamming") == 0)
		window = win_hamming;
	    /* Numerical Recipes 2nd ed. calls Hanning window "Hann window" */
	    else if (strcasecmp(argv[i], "Hann") == 0)
		window = win_hanning;
	    else if (strcasecmp(argv[i], "Hanning") == 0)
		window = win_hanning;
	    else if (strcasecmp(argv[i], "Parzen") == 0)
		window = win_parzen;
	    else if (strcasecmp(argv[i], "Square") == 0 ||
		     strcasecmp(argv[i], "Rectangular") == 0 ||
		     strcasecmp(argv[i], "Dirichlet") == 0)
		window = win_square;
	    else if (strcasecmp(argv[i], "Welch") == 0)
		window = win_welch;
	    else {
		fprintf(stderr, "%s: unrecognized window type %s\n",
			pname, argv[i]);
		exit(1);
	    }
	    wflag = 1;
	    break;
	  case 'z':	/* zero-mean the input */
	    zflag = 1;
	    break;
	  case 'Z':	/* zero-mean and detrend the input */
	    zflag = 2;
	    break;
	  default:
	    fprintf(stderr, "%s: unrecognized option %s ignored\n",
		    pname, argv[i]);
	    break;
	}
    }
    if (cflag) {
	if (fflag) {
	    fprintf(stderr, "%s: -c and -f are incompatible\n", pname);
	    exit(1);
	}
	if (pflag) {
	    fprintf(stderr, "%s: -c and -p are incompatible\n", pname);
	    exit(1);
	}
	if (Pflag) {
	    fprintf(stderr, "%s: -c and -P are incompatible\n", pname);
	    exit(1);
	}
    }
    if (nflag & pflag & Pflag) {
	fprintf(stderr, "%s: -n, -p, and -P are incompatible\n", pname);
	exit(1);
    }
    if (Nflag) {
	if (nflag) {
	    fprintf(stderr, "%s: -n and -N are incompatible\n", pname);
	    exit(1);
	}
	else
	    nflag = Nflag;
    }
    if (smooth > 1) {
        if (cflag) {
	    fprintf(stderr, "%s: -c and -s or -S are incompatible\n", pname);
	    exit(1);
	}
        if (pflag) {
	    fprintf(stderr, "%s: -p and -s or -S are incompatible\n", pname);
	    exit(1);
	}
    }

    /* Make sure that len is a power of two. */
    if (len < 1) len = 1;
    if (len < LEN) {
	for (m = LEN; m >= len; m >>= 1)
	    ;
	m <<= 1;
    }
    else {
	for (m = LEN; m < len; m <<= 1)
	    ;
    }
    len = m;

    if ((c = (float *)calloc(len, sizeof(float))) == NULL) {
	fprintf(stderr, "%s: insufficient memory\n", pname);
	exit(2);
    }