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