fft.c

快速傅利叶变换 C 语言源代码

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