Source Code 1992 March/Source_Code_CD-ROM_Walnut_Creek_March_1992.iso/msdos/turbo_c/fft.arc/FFT.C

C/C++ Source or Header | 1988-08-07 | 4.3 KB | 229 lines

/*
 *    fft.c
 *
 *    C Version 1.0 by Steve Sampson, Public Domain
 *
 *    This program is based on the work by W. D. Stanley
 *    and S. J. Peterson, Old Dominion University.
 *
 *    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_data  output_data
 *    Where 'samples' is a power of two
 *
 *    Array Version for Turbo C 1.5
 */

/* Includes */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <alloc.h>

/* Defines */

#define    TWO_PI    ((double)2.0 * M_PI)

/* Globals */

int    samples, power;
double    *real, *imag, max;
FILE    *fpi, *fpo;

/* Prototypes and forward declarations */

void    fft(void), max_amp(void), display(void);
int    permute(int);
double    magnitude(int);

/* The program */

main(argc, argv)
int    argc;
char    *argv[];
{
    int    n;

    if (argc != 4)  {
err1:        fprintf(stderr, "Usage: fft samples input output\n");
        fprintf(stderr, "Where samples is a power of 2\n");
        exit(1);
    }

    samples = abs(atoi(argv[1]));
    power =  log10((double)samples) / log10((double)2.0);

    if ((real = (double *)malloc(samples * sizeof(double))) == NULL)
err2:        fprintf(stderr, "Out of memory\n");

    if ((imag = (double *)malloc(samples * sizeof(double))) == NULL)
        goto err2;

    if ((fpo = fopen(argv[3], "wb")) == (FILE *)NULL)
        goto err1;

    if ((fpi = fopen(argv[2], "rb")) != (FILE *)NULL)  {
        fread(real, sizeof(double), samples, fpi);
        fread(imag, sizeof(double), samples, fpi);
        fclose(fpi);
    }
    else
        goto err1;

    fft();
    max_amp();
    display();

    fwrite(real, sizeof(double), samples, fpo);
    fwrite(imag, sizeof(double), samples, fpo);

    fclose(fpo);
}


void fft()
{
    unsigned i1, i2, i3, i4, y;
    int     loop, loop1, loop2;
    double     a1, a2, b1, b2, z1, z2, v;

    /* Scale the data */

    for (loop = 0; loop < samples; loop++)  {
        real[loop] /= (double)samples;
        imag[loop] /= (double)samples;
    }

    i1 = samples >> 1;
    i2 = 1;
    v = TWO_PI * ((double)1.0 / (double)samples);

    for (loop = 0; loop < power; loop++)  {
        i3 = 0;
        i4 = i1;

        for (loop1 = 0; loop1 < i2; loop1++)  {
            y = permute(i3 / i1);
            z1 =  cos(v * y);
            z2 = -sin(v * 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;
    }
}

/*
 *    Find maximum amplitude
 */

void max_amp()
{
    double    mag;
    int    loop;

    max = (double)0.0;
    for (loop = 0; loop < samples; loop++)  {
        if ((mag = magnitude(loop)) > max)
            max = mag;
    }
}

/*
 *    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.
 */

void display()
{
    int    c, n, x, loop;

    n = samples / 2;

    for (loop = n; loop < samples; loop++)  {
        x = (int)(magnitude(loop) * (double)56.0 / max);
        printf("%d\t|", loop - n);
        c = 0;
        while (++c <= x)
            putchar('=');

        putchar('\n');
    }

    for (loop = 0; loop < n; loop++)  {
        x = (int)(magnitude(loop) * (double)56.0 / max);
        printf("%d\t|", loop + n);
        c = 0;
        while (++c <= x)
            putchar('=');

        putchar('\n');
    }
}

/*
 *    Calculate Power Magnitude
 */

double magnitude(n)
int    n;
{
    n = permute(n);
    return (sqrt(real[n] * real[n] + imag[n] * imag[n]));
}

/*
 *    Bit reverse the number
 *
 *    Change 11100000b to 00000111b or vice-versa
 */

int permute(index)
int    index;
{
    int    n1, result, loop;

    n1 = samples;
    result = 0;

    for (loop = 0; loop < power; loop++)  {
        n1 >>= 1;            /* n1 / 2.0 */
        if (index < n1)
            continue;

        result += (int) pow((double)2.0, (double)loop);
        index -= n1;
    }

    return result;
}

/* EOF */