/* *	fft.c * *	Version 2.6 by Steve Sampson, Public Domain, November 1988 * *	This program produces a Frequency Domain display from the Time Domain *	data input; using the Fast Fourier Transform. * *	Runs in @ 30 seconds on a 5 MHz PC XT Clone without 8087. * *	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 input *		1. samples is 256. *		2. Input is (samples * sizeof(double)) characters. *		3. Standard error is help or debugging messages. * *	Auto detects CGA, EGA, and VGA in Turbo-C graphics mode. * *	See also: readme.doc, pulse.c, and sine.c. *//* Includes */#include <stdio.h>#include <alloc.h>#include <math.h>#include <conio.h>#include <dos.h>#include <bios.h>#include <graphics.h>#include <string.h>#include <stdlib.h>/* Defines */#define	SAMPLES		256#define	POWER		8	/* 2 to the 8th power is 256 */#define ESC		27	/* exit program */#define LEFTKEY		331	/* cursor left */#define RIGHTKEY	333	/* cursor right */#define CTRLLEFTKEY	371	/* cursor 10 left */#define CTRLRIGHTKEY	372	/* cursor 10 right */#define HOMEKEY		327	/* cursor at filter 0 */#define ENDKEY		335	/* cursor at filter 255 */#define F1		315	/* print screen */#define	TOP		50#define	LEFT		64#define	RIGHT		575#define	MIDDLE		192#define	TEXTX		160#define	TEXTXN		304#define	TEXTY		25/* *	Fixed constants should be used in the macros for speed. * *	A cosine wave leads a sine wave by 90 degrees, so offset into *	the lookup table 1/4 the way into it (256 / 4 = 64).  Using *	modulo 256, lookups will wrap around to zero for numbers greater *	than 255.  (cosine(200) = Sine[264 % 256] = Sine[8]). */#define	permute(x)	Br_table[(x)]#define	sine(x)		Sine[(x)]#define cosine(x)	Sine[((x) + 64) % 256]/* Globals, Forward declarations */double	Real[SAMPLES], Imag[SAMPLES], Maxn, magnitude();int	GraphDriver = DETECT, GraphMode, Primary, Cursor, Length;int	Bottom, Left, PrintChar(), PrintScreen(), getkey();void	*Save, quit(), beep(), build_window(), commands(), bee_bop();void	scale(), fft(), max_amp(), display();/* *	Bit Reverse Table for size 256 *	Lookup saves 20 seconds in Turbo-C over the pow() function. * *	Br_table[x] = x inverted (eg. 00000001 flipped to 10000000) */unsigned char Br_table[] = {	0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0, 0x10, 0x90, 0x50, 0xd0,	0x30, 0xb0, 0x70, 0xf0, 0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8,	0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8, 0x04, 0x84, 0x44, 0xc4,	0x24, 0xa4, 0x64, 0xe4, 0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4,	0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec, 0x1c, 0x9c, 0x5c, 0xdc,	0x3c, 0xbc, 0x7c, 0xfc, 0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2,	0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2, 0x0a, 0x8a, 0x4a, 0xca,	0x2a, 0xaa, 0x6a, 0xea, 0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa,	0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6, 0x16, 0x96, 0x56, 0xd6,	0x36, 0xb6, 0x76, 0xf6, 0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee,	0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe, 0x01, 0x81, 0x41, 0xc1,	0x21, 0xa1, 0x61, 0xe1, 0x11, 0x91, 0x51, 0xd1, 0x31, 0xb1, 0x71, 0xf1,	0x09, 0x89, 0x49, 0xc9, 0x29, 0xa9, 0x69, 0xe9, 0x19, 0x99, 0x59, 0xd9,	0x39, 0xb9, 0x79, 0xf9, 0x05, 0x85, 0x45, 0xc5, 0x25, 0xa5, 0x65, 0xe5,	0x15, 0x95, 0x55, 0xd5, 0x35, 0xb5, 0x75, 0xf5, 0x0d, 0x8d, 0x4d, 0xcd,	0x2d, 0xad, 0x6d, 0xed, 0x1d, 0x9d, 0x5d, 0xdd, 0x3d, 0xbd, 0x7d, 0xfd,	0x03, 0x83, 0x43, 0xc3, 0x23, 0xa3, 0x63, 0xe3, 0x13, 0x93, 0x53, 0xd3,	0x33, 0xb3, 0x73, 0xf3, 0x0b, 0x8b, 0x4b, 0xcb, 0x2b, 0xab, 0x6b, 0xeb,	0x1b, 0x9b, 0x5b, 0xdb, 0x3b, 0xbb, 0x7b, 0xfb, 0x07, 0x87, 0x47, 0xc7,	0x27, 0xa7, 0x67, 0xe7, 0x17, 0x97, 0x57, 0xd7, 0x37, 0xb7, 0x77, 0xf7,	0x0f, 0x8f, 0x4f, 0xcf, 0x2f, 0xaf, 0x6f, 0xef, 0x1f, 0x9f, 0x5f, 0xdf,	0x3f, 0xbf, 0x7f, 0xff};/* *	Sine/Cosine Table for size 256, Lookup saves 9 seconds in Turbo-C * *	Sine[n] = sin(x), x = x + (2Pi / 256) */float	Sine[] = {	0.000000, 0.024541, 0.049068, 0.073565,	0.098017, 0.122411,	0.146730, 0.170962, 0.195090, 0.219101, 0.242980, 0.266713,	0.290285, 0.313682, 0.336890, 0.359895, 0.382683, 0.405241,	0.427555, 0.449611, 0.471397, 0.492898, 0.514103, 0.534998,	0.555570, 0.575808, 0.595699, 0.615232,	0.634393, 0.653173,	0.671559, 0.689541, 0.707107, 0.724247, 0.740951, 0.757209,	0.773010, 0.788346, 0.803208, 0.817585, 0.831470, 0.844854,	0.857729, 0.870087, 0.881921, 0.893224, 0.903989, 0.914210,	0.923880, 0.932993, 0.941544, 0.949528,	0.956940, 0.963776,	0.970031, 0.975702, 0.980785, 0.985278, 0.989177, 0.992480,	0.995185, 0.997290, 0.998795, 0.999699,	1.000000, 0.999699,	0.998795, 0.997290, 0.995185, 0.992480, 0.989177, 0.985278,	0.980785, 0.975702, 0.970031, 0.963776,	0.956940, 0.949528,	0.941544, 0.932993, 0.923880, 0.914210, 0.903989, 0.893224,	0.881921, 0.870087, 0.857729, 0.844854,	0.831470, 0.817585,	0.803208, 0.788346, 0.773010, 0.757209, 0.740951, 0.724247,	0.707107, 0.689541, 0.671559, 0.653173,	0.634393, 0.615232,	0.595699, 0.575808, 0.555570, 0.534998, 0.514103, 0.492898,	0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895,	0.336890, 0.313682, 0.290285, 0.266713, 0.242980, 0.219101,	0.195090, 0.170962, 0.146730, 0.122411,	0.098017, 0.073565,	0.049068, 0.024541, 0.000000, -0.024541, -0.049068, -0.073565,	-0.098017, -0.122411, -0.146730, -0.170962, -0.195090, -0.219101,	-0.242980, -0.266713, -0.290285, -0.313682, -0.336890, -0.359895,	-0.382683, -0.405241, -0.427555, -0.449611, -0.471397, -0.492898,	-0.514103, -0.534998, -0.555570, -0.575808, -0.595699, -0.615232,	-0.634393, -0.653173, -0.671559, -0.689541, -0.707107, -0.724247,	-0.740951, -0.757209, -0.773010, -0.788346, -0.803208, -0.817585,	-0.831470, -0.844854, -0.857729, -0.870087, -0.881921, -0.893224,	-0.903989, -0.914210, -0.923880, -0.932993, -0.941544, -0.949528,	-0.956940, -0.963776, -0.970031, -0.975702, -0.980785, -0.985278,	-0.989177, -0.992480, -0.995185, -0.997290, -0.998795, -0.999699,	-1.000000, -0.999699, -0.998795, -0.997290, -0.995185, -0.992480,	-0.989177, -0.985278, -0.980785, -0.975702, -0.970031, -0.963776,	-0.956940, -0.949528, -0.941544, -0.932993, -0.923880, -0.914210,	-0.903989, -0.893224, -0.881921, -0.870087, -0.857729, -0.844854,	-0.831470, -0.817585, -0.803208, -0.788346, -0.773010, -0.757209,	-0.740951, -0.724247, -0.707107, -0.689541, -0.671559, -0.653173,	-0.634393, -0.615232, -0.595699, -0.575808, -0.555570, -0.534998,	-0.514103, -0.492898, -0.471397, -0.449611, -0.427555, -0.405241,	-0.382683, -0.359895, -0.336890, -0.313682, -0.290285, -0.266713,	-0.242980, -0.219101, -0.195090, -0.170962, -0.146730, -0.122411,	-0.098017, -0.073565, -0.049068, -0.024541};FILE	*Fpi, *Fpo;/* The program */main(argc, argv)int	argc;char	**argv;{	if (argc != 2)  {		fprintf(stderr, "Usage: fft input_file\n");		exit(1);	}	setcbrk(1);		/* Allow Control-C checking */	ctrlbrk(quit);		/* Execute quit() if Control-C detected */	/* open the data file */	if ((Fpi = fopen(*++argv, "rb")) == NULL)  {		fprintf(stderr,"fft: Unable to open data input file\n");		exit(1);	}	/* read in the data */	fread(Real, sizeof(double), SAMPLES, Fpi);	fread(Imag, sizeof(double), SAMPLES, Fpi);	fclose(Fpi);	build_window();	scale();	fft();	display();	/* wait for keyboard commands */	commands();}void scale(){	register int	loop;	for (loop = 0; loop < SAMPLES; loop++)  {		Real[loop] /= SAMPLES;		Imag[loop] /= SAMPLES;	}}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. */void display(){	register int	loop, offset;	int		n, x;	n = SAMPLES >> 1;	max_amp();	/*	 *	Graphics screen horizontal configuration:	 *	 *	| 64 pixels | 512 pixels | 64 pixels |	 *	|   margin  |   picture  |   margin  |	 *	 *	Spectral lines are two bits wide	 */	for (loop = n, offset = LEFT; loop < SAMPLES; loop++, offset++)  {		x = (int)(magnitude(loop) * Length / Maxn);		bar((offset + loop - n), Bottom - x, (offset + loop - n) + 1, Bottom);	}	for (loop = 0, offset = MIDDLE; loop < n; loop++, offset++)  {		x = (int)(magnitude(loop) * Length / Maxn);		bar((offset + loop + n), Bottom - x, (offset + loop + n) + 1, Bottom);	}}/* *	Find maximum amplitude */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 */double magnitude(n)int	n;{	n = permute(n);	return hypot(Real[n], Imag[n]);}void build_window(){	unsigned i_size;