fft.c.txt

This program produces a Frequency Domain display from the Time Domain * data input using the Fast

/* 
 *  fft.c 
 * 
 *  Unix Version 2.4 by Steve Sampson, Public Domain, September 1988 
 * 
 *  This program produces a Frequency Domain display from the Time Domain 
 *  data input; using the Fast Fourier Transform. 
 * 
 *  The Real data is generated by the in-phase (I) channel, and the 
 *  Imaginary data is produced by the quadrature-phase (Q) channel of 
 *  a Doppler Radar receiver.  The middle filter is zero Hz.  Closing 
 *  targets are displayed to the right, and Opening targets to the left. 
 * 
 *  Note: With Imaginary data set to zero the output is a mirror image. 
 * 
 *  Usage:  fft samples input output 
 *      1. samples is a variable power of two. 
 *      2. Input is (samples * sizeof(double)) characters. 
 *      3. Output is (samples * sizeof(double)) characters. 
 *      4. Standard error is help or debugging messages. 
 * 
 *  See also: readme.doc, pulse.c, and sine.c. 
 */  
  
/* Includes */  
  
#include <STDIO.H>   
#include <MALLOC.H>   
#include <MATH.H>   
  
/* Defines */  
  
#define TWO_PI  (2.0 * 3.14159265358979324) /* alias 360 deg */   
#define Chunk   (Samples * sizeof(double))   
  
#define sine(x)     Sine[(x)]   
#define cosine(x)   Sine[((x) + (Samples >> 2)) % Samples]   
  
/* Globals, Forward declarations */  
  
static int  Samples, Power;  
static double   *Real, *Imag, Maxn, magnitude();  
static void scale(), fft(), max_amp(), display(), err_report();  
  
static int  permute();  
static double   *Sine;  
static void build_trig();  
  
static FILE *Fpi, *Fpo;  
  
  
/* The program */  
  
main(argc, argv)  
int argc;  
char    **argv;  
{  
    if (argc != 4)  
        err_report(0);  
  
    Samples = abs(atoi(*++argv));  
    Power = log10((double)Samples) / log10(2.0);  
  
    /* Allocate memory for the variable arrays */  
  
    if ((Real = (double *)malloc(Chunk)) == NULL)  
        err_report(1);  
  
    if ((Imag = (double *)malloc(Chunk)) == NULL)  
        err_report(2);  
  
    if ((Sine = (double *)malloc(Chunk)) == NULL)  
        err_report(3);  
  
    /* open the data files */  
  
    if ((Fpi = fopen(*++argv, "r")) == NULL)  
        err_report(4);  
  
    if ((Fpo = fopen(*++argv, "w")) == NULL)  
        err_report(5);  
  
    /* read in the data from the input file */  
  
    fread(Real, sizeof(double), Samples, Fpi);  
    fread(Imag, sizeof(double), Samples, Fpi);  
    fclose(Fpi);  
  
    build_trig();  
    scale();  
    fft();  
    display();  
  
    /* write the frequency domain data to the standard output */  
  
    fwrite(Real, sizeof(double), Samples, Fpo);  
    fwrite(Imag, sizeof(double), Samples, Fpo);  
    fclose(Fpo);  
  
    /* free up memory and let's get back to our favorite shell */  
  
    free((char *)Real);  
    free((char *)Imag);  
    free((char *)Sine);  
  
    exit(0);  
}  
  
  
static void err_report(n)  
int n;  
{  
    switch (n)  {  
    case 0:  
        fprintf(stderr, "Usage: fft samples in_file out_file\n");  
        fprintf(stderr, "Where samples is a power of two\n");  
        break;  
    case 1:  
        fprintf(stderr, "fft: Out of memory getting real space\n");  
        break;  
    case 2:  
        fprintf(stderr, "fft: Out of memory getting imag space\n");  
        free((char *)Real);  
        break;  
    case 3:  
        fprintf(stderr, "fft: Out of memory getting sine space\n");  
        free((char *)Real);  
        free((char *)Imag);  
        break;  
    case 4:  
        fprintf(stderr,"fft: Unable to open data input file\n");  
        free((char *)Real);  
        free((char *)Imag);  
        free((char *)Sine);  
        break;  
    case 5:  
        fprintf(stderr,"fft: Unable to open data output file\n");  
        fclose(Fpi);  
        free((char *)Real);  
        free((char *)Imag);  
        free((char *)Sine);  
        break;  
    }  
  
    exit(1);  
}  
  
  
static void scale()  
{  
    register int    loop;  
  
    for (loop = 0; loop < Samples; loop++)  {  
        Real[loop] /= Samples;  
        Imag[loop] /= Samples;  
    }  
}  
  
  
static void fft()  
{  
    register int    loop, loop1, loop2;  
    unsigned    i1;         /* going to right shift this */  
    int     i2, i3, i4, y;  
    double      a1, a2, b1, b2, z1, z2;  
  
    i1 = Samples >> 1;  
    i2 = 1;  
  
    /* perform the butterfly's */  
  
    for (loop = 0; loop < Power; loop++)  {  
        i3 = 0;  
        i4 = i1;  
  
        for (loop1 = 0; loop1 < i2; loop1++)  {  
            y = permute(i3 / (int)i1);  
            z1 =  cosine(y);  
            z2 = -sine(y);  
  
            for (loop2 = i3; loop2 < i4; loop2++)  {  
  
                a1 = Real[loop2];  
                a2 = Imag[loop2];  
  
                b1  = z1*Real[loop2+i1] - z2*Imag[loop2+i1];  
                b2  = z2*Real[loop2+i1] + z1*Imag[loop2+i1];  
  
                Real[loop2] = a1 + b1;  
                Imag[loop2] = a2 + b2;  
  
                Real[loop2+i1] = a1 - b1;  
                Imag[loop2+i1] = a2 - b2;  
            }  
  
            i3 += (i1 << 1);  
            i4 += (i1 << 1);  
        }  
  
        i1 >>= 1;  
        i2 <<= 1;  
    }  
}  
  
/* 
 *  Display the frequency domain. 
 * 
 *  The filters are aranged so that DC is in the middle filter. 
 *  Thus -Doppler is on the left, +Doppler on the right. 
 */  
  
static void display()  
{  
    register int    loop, offset;  
    int     n, x;  
  
    max_amp();  
    n = Samples >> 1;  
  
    for (loop = n; loop < Samples; loop++)  {  
        x = (int)(magnitude(loop) * 59.0 / Maxn);  
        printf("%d\t|", loop - n);  
        offset = 0;  
        while (++offset <= x)  
            putchar('=');  
  
        putchar('\n');  
    }  
  
    for (loop = 0; loop < n; loop++)  {  
        x = (int)(magnitude(loop) * 59.0 / Maxn);  
        printf("%d\t|", loop + n);  
        offset = 0;  
        while (++offset <= x)  
            putchar('=');  
  
        putchar('\n');  
    }  
}  
  
/* 
 *  Find maximum amplitude 
 */  
  
static void max_amp()  
{  
    register int    loop;  
    double      mag;  
  
    Maxn = 0.0;  
    for (loop = 0; loop < Samples; loop++)  {  
        if ((mag = magnitude(loop)) > Maxn)  
            Maxn = mag;  
    }  
}  
  
/* 
 *  Calculate Power Magnitude 
 */  
  
static double magnitude(n)  
int n;  
{  
    n = permute(n);  
    return hypot(Real[n], Imag[n]);  
}  
  
/* 
 *  Bit reverse the number 
 * 
 *  Change 11100000b to 00000111b or vice-versa 
 */  
  
static int permute(index)  
int index;  
{  
    register int    loop;  
    unsigned    n1;  
    int     result;  
  
    n1 = Samples;  
    result = 0;  
  
    for (loop = 0; loop < Power; loop++)  {  
        n1 >>= 1;  
        if (index < n1)  
            continue;  
  
        /* Unix has a truncation problem with pow() */  
  
        result += (int)(pow(2.0, (double)loop) + .05);  
        index -= n1;  
    }  
  
    return result;  
}  
  
/* 
 *  Pre-compute the sine and cosine tables 
 */  
  
static void build_trig()  
{  
    register int    loop;  
    double      angle, increment;  
  
    angle = 0.0;  
    increment = TWO_PI / (double)Samples;  
  
    for (loop = 0; loop < Samples; loop++)  {  
        Sine[loop] = sin(angle);  
        angle += increment;  
    }  
}  
  
/* EOF */