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