/*
 *	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>
#ifdef AZTEC_C
#include <strings.h>
#else
#include <malloc.h>
#endif
#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;
	}
}