www.pudn.com > 2DFFT-C.zip > FFT.C

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