fft.c

Fast Fourier Transform in C

	
    if (iflag) {		/* calculate and print inverse FFT */
	for (n = 0, rsum = 0.; n < len && fscanf(ifile, "%f", &c[n]) == 1; n++)
	    ;
	if (n == 0) {
	    fprintf(stderr, "%s: standard input is empty\n", pname);
	    exit(2);
	}
	ifft();
	exit(0);
    }

    else {			/* calculate and print forward FFT */
	if (nflag) {		/* process input in chunks */
	    float *s, *t;
	    int nf2 = nflag/2;

	    for (m = len; m >= nflag; m >>= 1)
		;
	    m <<= 1;		/* m is now the smallest power of 2 >= nflag */
	    if ((s = (float *)calloc(sizeof(float), m)) == NULL ||
		(t = (float *)calloc(sizeof(float), m/2)) == NULL) {
		fprintf(stderr, "%s: insufficient memory\n", pname);
		exit(2);
	    }
	    for (n = 0; n < nf2 && fscanf(ifile, "%f", &t[n]) == 1; n++)
		;
	    while (1) {
		if (zflag) {
		    for (n = 0, rsum = 0.; n < nf2; n++)
			rsum += c[n] = t[n];
		    for (i=0; n<nflag && fscanf(ifile,"%f",&t[i])==1; i++, n++)
			rsum += c[n] = t[i]; /* read input, accumulate sum */
		}
		else {
		    for (n = 0, rsum = 0.; n < nf2; n++)
			c[n] = t[n];
		    for (i=0; n<nflag && fscanf(ifile,"%f",&t[i])==1; i++, n++)
			c[n] = t[i];
		    ;
		}
		if (n < nflag) break;
		fft();
		if (Nflag)
		    fft_out();
		else if (Pflag)
		    for (i = 0; i < m; i++)
			s[i] += c[i]*c[i];
		else
		    for (i = 0; i < m; i++)
			s[i] += c[i];
		nchunks++;
	    }
	    if (nchunks < 1) {
		fprintf(stderr, "%s: input series is too short\n", pname);
		exit(2);
	    }
	    if (Nflag == 0 && Pflag)
		for (i = 0; i < m; i++)
		    c[i] = sqrt(s[i])/nchunks;
	    else
		for (i = 0; i < m; i++)
		    c[i] = s[i]/nchunks;
	}

	else {
	    read_input();
	    if (n == 0) {
		fprintf(stderr, "%s: standard input is empty\n", pname);
		exit(2);
	    }
	    fft();
	}

	if (!Nflag)
	    fft_out();

	exit(0);
    }
}

/* This function detrends (subtracts a least-squares fitted line from) a
   a sequence of n uniformily spaced ordinates supplied in c. */
detrend(c, n)
float *c;
int n;
{
    int i;
    double a, b = 0.0, tsqsum = 0.0, ysum = 0.0, t;

    for (i = 0; i < n; i++)
	ysum += c[i];
    for (i = 0; i < n; i++) {
	t = i - n/2 + 0.5;
	tsqsum += t*t;
	b += t*c[i];
    }
    b /= tsqsum;
    a = ysum/n - b*(n-1)/2.0;
    for (i = 0; i < n; i++)
	c[i] -= a + b*i;
    if (b < -0.04 || b > 0.04)
	fprintf(stderr,
		"%s: (warning) possibly significant trend in input series\n",
		pname);
}

read_input()
{
    if (zflag)
	for (n = 0, rsum = 0.; n < len && fscanf(ifile, "%f", &c[n]) == 1; n++)
	    rsum += c[n];	/* read input, accumulate sum */
    else
	for (n = 0, rsum = 0.; n < len && fscanf(ifile, "%f", &c[n]) == 1; n++)
	    ;
}

fft()		/* calculate forward FFT */
{
    int i;

    if (zflag) {		/* zero-mean the input array */
	rmean = rsum/n;
	for (i = 0; i < n; i++)
	    c[i] -= rmean;	
	if (zflag == 2)
	    detrend(c, n);
    }
    for (m = len; m >= n; m >>= 1)
	;
    m <<= 1;		/* m is now the smallest power of 2 >= n; this is the
			   length of the input series (including padding) */
    if (wflag)			/* apply the chosen windowing function */
	for (i = 0; i < m; i++)
	    c[i] *= (*window)(i, m);
    else
	wsum = m;
    norm = sqrt(2.0/(wsum*n));
    if (fflag) fstep = freq/(2.*m); /* note that fstep is actually half of
				       the frequency interval;  it is
				       multiplied by the doubled index i
				       to obtain the center frequency for
				       bin (i/2) */
    realft(c-1, m/2, 1);	/* perform the FFT;  see Numerical Recipes */
}

fft_out()	/* print the FFT */
{
    int i;

    c[m] = c[1];		/* unpack the output array */
    c[1] = c[m+1] = 0.;
    for (i = 0; i <= m; i += 2*decimation) {
	int j;
	double pow;

	if (fflag) printf("%g\t", i*fstep);
	if (cflag) printf("%g\t%g\n", c[i], c[i+1]);
	else {
	    for (j = 0, pow = 0.0; j < 2*smooth; j += 2)
	        pow += (c[i+j]*c[i+j] + c[i+j+1]*c[i+j+1])*norm*norm;
	    pow /= smooth/decimation;
	    if (Pflag) printf("%g", pow);
	    else printf("%g", sqrt(pow));
	    if (pflag) printf("\t%g", atan2(c[i+1], c[i]));
	    printf("\n");
	}
    }
}

ifft()		/* calculate and print inverse FFT */
{
    int i;

    n -= 2;
    c[1] = c[n];		/* repack IFFT input array */
    if (iflag < 0) {		/* convert polar form input to rectangular */
	for (i = 2; i < n; i += 2) {
	    float im;

	    im = c[i]*sin(c[i+1]);
	    c[i] *= cos(c[i+1]);
	    c[i+1] = im;
	}
    }
    realft(c-1, n/2, -1);
    if (iflag < 0) {
	norm = sqrt(2.0);
	for (i = 0; i < n; i++)
	    printf("%g\n", c[i]*norm);
    }
    else
	for (i = 0; i < n; i++)
	    printf("%g\n", c[i]/(n/2.0));
}

char *prog_name(s)
char *s;
{
    char *p = s + strlen(s);

#ifdef MSDOS
    while (p >= s && *p != '\\' && *p != ':') {
	if (*p == '.')
	    *p = '\0';		/* strip off extension */
	if ('A' <= *p && *p <= 'Z')
	    *p += 'a' - 'A';	/* convert to lower case */
	p--;
    }
#else
    while (p >= s && *p != '/')
	p--;
#endif
    return (p+1);
}

static char *help_strings[] = {
 "usage: %s [ OPTIONS ...] INPUT-FILE\n",
 " where INPUT-FILE is the name of a text file containing a time series",
 " (use `-' to read the standard input), and OPTIONS may be any of:",
 " -c       Output unnormalized complex FFT (real components in first column,",
 "          imaginary components in second column).",
 " -f FREQ  Show the center frequency for each bin in the first column.  The",
 "          FREQ argument specifies the input sampling frequency;  the center",
 "          frequencies are given in the same units.",
 " -h       Print on-line help.",
 " -i       Perform inverse FFT;  in this case, the standard input should be",
 "          in the form generated by `fft -c', and the standard output is",
 "          a series of samples.",
 " -I       Perform inverse FFT as above, but using input generated by `fft -p'.",
 " -l LEN   Perform up to LEN-point transforms.  `fft' rounds n up to the next",
 "          higher power of two unless LEN is already a power of two.  If the",
 "          input series contains fewer than LEN samples, it is padded with",
 "          zeros up to the next higher power of two.  Any additional input",
 "          samples beyond the first LEN are not read.  Default: LEN = 16384.",
 " -n NN     Process the input in overlapping chunks of N samples and output",
 "          an averaged spectrum.  If used in combination with -P, the output",
 "          is the average of the individual squared magnitudes;  otherwise,",
 "          the output is derived from the averages of the real components and",
 "          of the imaginary components taken separately.  For best results,",
 "          NN should be a power of two.",
 " -N NN    Process the input in overlapping chunks of NN samples and output a",
 "          spectrum for each chunk.  For best results, NN should be a power",
 "          of two.",
 " -p       Show the phase in radians in the last column.",
 " -P       Generate a power spectrum (print squared magnitudes).",
 " -s N     Smooth the output by applying an N-point moving average to each bin.",
 " -S N     Smooth the output by summing sets of N consecutive bins.",
 " -w WINDOW",
 "          Apply the specified WINDOW to the input data.  WINDOW may be one",
 "          of: `Bartlett', `Blackman', `Blackman-Harris', `Hamming',",
 "          `Hanning', `Parzen', `Square', and `Welch'.",
 " -z       Zero-mean the input data.",
 " -Z       Detrend and zero-mean the input data.",
 NULL
};

void help()
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++) {
	(void)fprintf(stderr, "%s\n", help_strings[i]);
	if (i % 23 == 0) {
	    char b[5];
	    (void)fprintf(stderr, "--More--");
	    (void)fgets(b, 5, stdin);
	    (void)fprintf(stderr, "\033[A\033[2K"); /* erase "--More--";
						       assumes ANSI terminal */
	}
    }
}

void realft(data,n,isign)
float data[];
int n,isign;
{
    int i, i1, i2, i3, i4, n2p3;
    float c1 = 0.5, c2, h1r, h1i, h2r, h2i;
    double wr, wi, wpr, wpi, wtemp, theta;
    void four1();

    theta = PI/(double) n;
    if (isign == 1) {
	c2 = -0.5;
	four1(data, n, 1);
    } 
    else {
	c2 = 0.5;
	theta = -theta;
    }
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0+wpr;
    wi = wpi;
    n2p3 = 2*n+3;
    for (i = 2; i <= n/2; i++) {
	i4 = 1 + (i3 = n2p3 - (i2 = 1 + ( i1 = i + i - 1)));
	h1r =  c1*(data[i1] + data[i3]);
	h1i =  c1*(data[i2] - data[i4]);
	h2r = -c2*(data[i2] + data[i4]);
	h2i =  c2*(data[i1] - data[i3]);
	data[i1] =  h1r + wr*h2r - wi*h2i;
	data[i2] =  h1i + wr*h2i + wi*h2r;
	data[i3] =  h1r - wr*h2r + wi*h2i;
	data[i4] = -h1i + wr*h2i + wi*h2r;
	wr = (wtemp = wr)*wpr - wi*wpi+wr;
	wi = wi*wpr + wtemp*wpi + wi;
    }
    if (isign == 1) {
	data[1] = (h1r = data[1]) + data[2];
	data[2] = h1r - data[2];
    } else {
	data[1] = c1*((h1r = data[1]) + data[2]);
	data[2] = c1*(h1r - data[2]);
	four1(data, n, -1);
    }
}

void four1(data, nn, isign)
float data[];
int nn, isign;
{
    int n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    float tempr, tempi;
    
    n = nn << 1;
    j = 1;
    for (i = 1; i < n; i += 2) {
	if (j > i) {
	    tempr = data[j];     data[j] = data[i];     data[i] = tempr;
	    tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;
	}
	m = n >> 1;
	while (m >= 2 && j > m) {
	    j -= m;
	    m >>= 1;
	}
	j += m;
    }
    mmax = 2;
    while (n > mmax) {
	istep = 2*mmax;
	theta = TWOPI/(isign*mmax);
	wtemp = sin(0.5*theta);
	wpr = -2.0*wtemp*wtemp;
	wpi = sin(theta);
	wr = 1.0;
	wi = 0.0;
	for (m = 1; m < mmax; m += 2) {
	    for (i = m; i <= n; i += istep) {
		j =i + mmax;
		tempr = wr*data[j]   - wi*data[j+1];
		tempi = wr*data[j+1] + wi*data[j];
		data[j]   = data[i]   - tempr;
		data[j+1] = data[i+1] - tempi;
		data[i] += tempr;
		data[i+1] += tempi;
	    }
	    wr = (wtemp = wr)*wpr - wi*wpi + wr;
	    wi = wi*wpr + wtemp*wpi + wi;
	}
	mmax = istep;
    }
}