FFT.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;
/* settup graphics mode */
registerbgidriver(CGA_driver);
registerbgidriver(EGAVGA_driver);
initgraph(&GraphDriver, &GraphMode, "");
if (graphresult() != grOk) {
fprintf(stderr, "fft: %s\n", grapherrormsg(graphresult()));
exit(1);
}
if ((i_size = getgraphmode()) != CGAHI)
setbkcolor(BLACK);
switch (i_size) {
case CGAHI:
Primary = WHITE;
Cursor = WHITE;
Length = 100.0;
Bottom = 150;
break;
case EGAHI:
Primary = LIGHTGRAY;
Cursor = LIGHTGREEN;
Length = 250.0;
Bottom = 300;
break;
case VGAHI:
Primary = LIGHTGRAY;
Cursor = LIGHTGREEN;
Length = 380.0;
Bottom = 430;
}
setcolor(Primary);
setfillstyle(SOLID_FILL, Primary);
settextjustify(LEFT_TEXT, TOP_TEXT);
settextstyle(DEFAULT_FONT, HORIZ_DIR, 2);
cleardevice();
outtextxy(TEXTX, TEXTY, "Computing FFT");
/* Draw ruler line */
line(LEFT, Bottom + 5, RIGHT, Bottom + 5);
for (i_size = LEFT; i_size <= RIGHT ; i_size += 10) {
line(i_size, Bottom + 5, i_size, Bottom + 10);
line(i_size + 1, Bottom + 5, i_size + 1, Bottom + 10);
}
/* create under-cursor memory */
i_size = imagesize(LEFT, TOP, LEFT + 1, Bottom);
Save = malloc(i_size);
}
void commands()
{
int key, filter;
char string[4];
setcolor(BLACK); /* erase text */
outtextxy(TEXTX, TEXTY, "Computing FFT");
setcolor(Primary);
Left = 320;
filter = 128;
strcpy(string, "128");
outtextxy(TEXTX, TEXTY, "Filter # 128");
getimage(Left, TOP, Left + 1, Bottom, Save);
setfillstyle(SOLID_FILL, Cursor);
bar(Left, TOP, Left + 1, Bottom);
for (;;) {
while (bioskey(1) == 0) /* wait for key press */
;
key = getkey();
switch(key) { /* find out which key was pressed */
case HOMEKEY:
j1: putimage(Left, TOP, Save, COPY_PUT);
filter = 0;
Left = LEFT;
getimage(Left, TOP, Left + 1, Bottom, Save);
bar(Left, TOP, Left + 1, Bottom);
setcolor(BLACK); /* erase old text */
outtextxy(TEXTXN, TEXTY, string);
setcolor(Primary);
sprintf(string, "%d", filter);
outtextxy(TEXTXN, TEXTY, string);
break;
case ENDKEY:
j2: putimage(Left, TOP, Save, COPY_PUT);
filter = 255;
Left = RIGHT - 1;
getimage(Left, TOP, Left + 1, Bottom, Save);
bar(Left, TOP, Left + 1, Bottom);
setcolor(BLACK); /* erase old text */
outtextxy(TEXTXN, TEXTY, string);
setcolor(Primary);
sprintf(string, "%d", filter);
outtextxy(TEXTXN, TEXTY, string);
break;
case CTRLLEFTKEY:
if (filter <= 9)
goto j1;
putimage(Left, TOP, Save, COPY_PUT);
Left -= 20;
filter -= 10;
getimage(Left, TOP, Left + 1, Bottom, Save);
bar(Left, TOP, Left + 1, Bottom);
setcolor(BLACK); /* erase old text */
outtextxy(TEXTXN, TEXTY, string);
setcolor(Primary);
sprintf(string, "%d", filter);
outtextxy(TEXTXN, TEXTY, string);
break;
case LEFTKEY:
if (filter <= 0) {
beep();
break;
}
putimage(Left, TOP, Save, COPY_PUT);
Left -= 2;
filter--;
getimage(Left, TOP, Left + 1, Bottom, Save);
bar(Left, TOP, Left + 1, Bottom);
setcolor(BLACK); /* erase old text */
outtextxy(TEXTXN, TEXTY, string);
setcolor(Primary);
sprintf(string, "%d", filter);
outtextxy(TEXTXN, TEXTY, string);
break;
case CTRLRIGHTKEY:
if (filter >= 245)
goto j2;
putimage(Left, TOP, Save, COPY_PUT);
Left += 20;
filter += 10;
getimage(Left, TOP, Left + 1, Bottom, Save);
bar(Left, TOP, Left + 1, Bottom);
setcolor(BLACK); /* erase old text */
outtextxy(TEXTXN, TEXTY, string);
setcolor(Primary);
sprintf(string, "%d", filter);
outtextxy(TEXTXN, TEXTY, string);
break;
case RIGHTKEY:
if (filter >= 255) {
beep();
break;
}
putimage(Left, TOP, Save, COPY_PUT);
Left += 2;
filter++;
getimage(Left, TOP, Left + 1, Bottom, Save);
bar(Left, TOP, Left + 1, Bottom);
setcolor(BLACK); /* erase old text */
outtextxy(TEXTXN, TEXTY, string);
setcolor(Primary);
sprintf(string, "%d", filter);
outtextxy(TEXTXN, TEXTY, string);
break;
case F1:
if (PrintScreen())
bee_bop();
break;
case ESC:
quit();
default:
beep();
}
}
}
void beep()
{
sound(200);
delay(125);
nosound();
}
void bee_bop()
{
sound(1000);
delay(80);
sound(200);
delay(160);
sound(1000);
nosound();
}
/*
* Return to MS-DOS operating system
*/
void quit()
{
sound(1000);
delay(250);
nosound();
closegraph();
free((char *)Save);
exit(0);
}
/*
* Dumps a screen to printer in Double Density landscape format.
* This routine is from the Borland Programming Forum on CompuServe.
*
* returns TRUE on error
*/
int PrintScreen()
{
int X,Y,YOfs;
int BitData,MaxBits;
unsigned int Height,Width;
struct viewporttype Vport;
getviewsettings(&Vport);
Width = Vport.bottom - Vport.top;
Height = (Vport.right + 1) - Vport.left;
if (PrintChar(27))
return (1);
if (PrintChar('3'))
return (1);
if (PrintChar(24)) /* 24/216 inch line spacing */
return (1);
Y = 0;
while (Y < Height) {
if (PrintChar(27))
return (1);
if (PrintChar('L')) /* Double Density */
return (1);
if (PrintChar(Width & 0xff))
return (1);
if (PrintChar(Width >> 8))
return (1);
for (X = Width - 1; X >= 0; X--) {
BitData = 0;
if ((Y + 7) <= Height)
MaxBits = 7;
else
MaxBits = Height - Y;
for (YOfs = 0; YOfs <= MaxBits; YOfs++) {
BitData = BitData << 1;
if (getpixel(YOfs + Y, X) > 0)
BitData++;
}
if (PrintChar(BitData))
return (1);
}
if (PrintChar('\r'))
return (1);
if (PrintChar('\n'))
return (1);
Y += 8;
}
PrintChar(12); /* formfeed */
return (0);
}
/*
* returns TRUE on error
*/
int PrintChar(out)
int out;
{
unsigned char ret;
ret = biosprint(0, out, 0);
if ((ret & 0x29) || !(ret & 0x10))
return 1;
else
return 0;
}
int getkey()
{
int key, lo, hi;
key = bioskey(0);
lo = key & 0x00FF;
hi = (key & 0xFF00) >> 8;
return((lo == 0) ? hi + 256 : lo);
}
/* EOF */
类似资料：