Source Code 1992 March

home *** CD-ROM | disk | FTP | other *** search

/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / compsrcs / misc / volume05 / fft2 < prev next >

Wrap

Internet Message Format | 1991-08-27 | 43.8 KB

From decwrl!purdue!mailrus!tut.cis.ohio-state.edu!cwjcc!hal!ncoast!allbery Sat Dec 17 21:43:57 PST 1988 Article 756 of comp.sources.misc: Path: granite!decwrl!purdue!mailrus!tut.cis.ohio-state.edu!cwjcc!hal!ncoast!allbery From: sampson@killer.DALLAS.TX.US (Steve Sampson) Newsgroups: comp.sources.misc Subject: v05i080: Unix and Turbo-C General Purpose FFT Message-ID: <6370@killer.DALLAS.TX.US> Date: 14 Dec 88 02:53:36 GMT Sender: allbery@ncoast.UUCP Reply-To: sampson@killer.DALLAS.TX.US (Steve Sampson) Organization: The Unix(R) Connection, Dallas, Texas Lines: 1824 Approved: allbery@ncoast.UUCP Posting-number: Volume 5, Issue 80 Submitted-by: "Steve Sampson" <sampson@killer.DALLAS.TX.US> Archive-name: fft2 This archive contains source for both a Unix and MS-DOS General Purpose FFT program. I've posted earlier versions before, but consider them obsolete. This all can be improved on I'm sure. Have fun with it. #! /bin/sh # Contents: uread.doc ufft.c usine.c upulse.c read.doc makefile fft.c # sine.c pulse.c # Wrapped by sampson@killer.DALLAS.TX on Fri Dec 09 23:04:22 1988 PATH=/bin:/usr/bin:/usr/ucb ; export PATH if test -f uread.doc -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"uread.doc\" else echo shar: Extracting \"uread.doc\" $5220 characters$ sed "s/^X//" >uread.doc <<'END_OF_uread.doc' X5 September 1988 X XFFT Program for Unix C compiler XSteve Sampson, Box 45668, Tinker AFB, OK 73145 Xsampson@killer (TAC Trained) X X X This Unix version of the FFT program uses a general purpose display routine, Xand the number of samples is based on the amount of memory. This version can Xalso be used in MS-DOS (if you change the file modes to "rb" and "wb"). X XThe original Byte Magazine program (see references below) was designed for real Xdata only. In my experiments I needed to preserve both real and imaginary Xdata. If you feed the FFT real data only, then the output will be a mirror Ximage, and you can ignore the left side. Two signal generators are included. XOne generates sine waves (sine) and the other generates pulses (pulse). Some Xpapers I found on the subject of FFTs are included at the end. There are Xseveral books devoted to the subject also. X XFor the Unix example try: X X sine 16 in X 1000 X 3000 X XWhich will sample the 1 Khz data every 333 microseconds (1 / 3 Khz). Note: The Xsample frequency should be greater than 2 times the input frequency (Nyquist Xand all that...). X XThen run fft: X X fft 16 in out X XAnd you should see a display like so: X X0 |======= (-1500.0 Hz) X1 |===== (-1312.5 Hz) X2 |==== (-1125.0 Hz) X3 |==== (-937.0 Hz) X4 |=== (-750.0 Hz) X5 |=== (-562.5 Hz) X6 |=== (-375.0 Hz) X7 |=== (-187.5 Hz) X8 |==== <------- DC (000.0 Hz) X9 |==== <------- Fundamental (187.5 Hz) X10 |====== <------- Second Harmonic (375.0 Hz) X11 |======== (562.5 Hz) X12 |============== (750.0 Hz) X13 |======================================================== X14 |============================ (1125.0 Hz) ^ X15 |=========== (1312.5 Hz) | X | X [13 - 8 (center)] * 187.5 = 937.0 Hz X XThe fundamental display frequency is: X X T = Time Increment Between Samples X N = Number Of Samples X Tp = N * T X X Then F = 1 / Tp X X In the example above, the time increment between samples is X 1 / 3000 or 333 microseconds. N = 16, so Tp = 5333 microseconds X and 1 / .005333 is 187.5 Hz. X X Therefore each filter is a multiple of 187.5 Hertz. Filter 8 in this X example is center, so that would be zero, 9 would be one, etc. X XIn this case the sample interval didn't work so good for the frequency and Xshows the limitation of the Discrete Fourier Transform in representing a Xcontinuous signal. A better sample rate for 1000 Hz would be 4000 Hz, Xin which case T = 250 ms, Tp = 4 ms, and F = 250 Hz. 1000 / 250 = 4. The Xpower should all be in filter 12 (8 + 4) in this example. X XLet's run it and see: X X sine 16 in X 1000 X 4000 X X fft 16 in out X X0 | X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | X9 | X10 | X11 | X12 |======================================================== X13 | X14 | X15 | X XWell what do you know... X XThe output data file can be used by other programs as needed. X XBy using negative frequencies in 'sine' you can generate opening targets: X X sine 16 in X -1000 X 3000 X X fft 16 in out X XProduces: X X0 |======= X1 |=========== X2 |============================ X3 |======================================================= X4 |============== X5 |======== X6 |====== X7 |==== X8 |==== <-------- Zero Hertz (DC) X9 |=== X10 |=== X11 |=== X12 |=== X13 |==== X14 |==== X15 |===== X XYou can see in these examples where weighting functions would be nice. X XFor an example of what happens when the imaginary data is not input X(ie, zeros put in) for a 1000 Hz frequency at 3000 Hz sample rate: X X0 |=============== X1 |================== X2 |=================================== X3 |======================================================== X4 |=========== X5 |==== X6 |== X7 |= Delete this part X--------------------------------------------------------------------- X8 | X9 |= X10 |== X11 |==== X12 |=========== X13 |======================================================= X14 |=================================== X15 |================== X XThe left side is redundant and can be deleted. This is what the original XByte Magazine article did. X XFor generating pulses, a second program 'pulse' is provided. It pre-loads Ximaginary data with zeros. For example: X X pulse 16 in X .000006 X .0000008 X XIs a radar with a 6 microsecond pulse and 800 nanosecond range gates. X X fft 16 in out X XWill produce: X X0 | X1 |======= X2 | X3 |======== X4 | X5 |============ X6 | X7 |=================================== X8 |======================================================== X9 |=================================== X10 | X11 |============ X12 | X13 |======== X14 | X15 |======= X XThe more filters you use, the prettier it looks. X X XFFT References X-------------- X X1. Fast Fourier Transforms On Your Home Computer, William D. Stanley, X Steven J. Peterson, BYTE Magazine, December 1978. Basic idea comes X from this program. X X2. 8052 Microcomputer simplifies FFT Design, Arnold Rosenberg, X Electronics, May 5, 1983. Used a bit reverse table idea based on the X routine in this program. X X3. A Fast Fourier Transform for the 8080, Robert D. Fusfeld, X Dr. Dobbs, Number 44. Gave me some ideas. X X4. A Guided Tour of the Fast Fourier Transform, G. D. Bergland, X IEEE Spectrum, July 1969. Math! X X5. FFT - Shortcut to Fourier Analysis, Richard Klahn, Richard R. Shively, X Electronics, April 15 1968. Math! X X/* EOF */ END_OF_uread.doc if test 5220 -ne `wc -c <uread.doc`; then echo shar: \"uread.doc\" unpacked with wrong size! fi # end of overwriting check fi if test -f ufft.c -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"ufft.c\" else echo shar: Extracting \"ufft.c\" $5897 characters$ sed "s/^X//" >ufft.c <<'END_OF_ufft.c' X/* X * fft.c X * X * Unix Version 2.4 by Steve Sampson, Public Domain, September 1988 X * X * This program produces a Frequency Domain display from the Time Domain X * data input; using the Fast Fourier Transform. X * X * The Real data is generated by the in-phase (I) channel, and the X * Imaginary data is produced by the quadrature-phase (Q) channel of X * a Doppler Radar receiver. The middle filter is zero Hz. Closing X * targets are displayed to the right, and Opening targets to the left. X * X * Note: With Imaginary data set to zero the output is a mirror image. X * X * Usage: fft samples input output X * 1. samples is a variable power of two. X * 2. Input is (samples * sizeof(double)) characters. X * 3. Output is (samples * sizeof(double)) characters. X * 4. Standard error is help or debugging messages. X * X * See also: readme.doc, pulse.c, and sine.c. X */ X X/* Includes */ X X#include <stdio.h> X#include <malloc.h> X#include <math.h> X X/* Defines */ X X#define TWO_PI (2.0 * 3.14159265358979324) /* alias 360 deg */ X#define Chunk (Samples * sizeof(double)) X X#define sine(x) Sine[(x)] X#define cosine(x) Sine[((x) + (Samples >> 2)) % Samples] X X/* Globals, Forward declarations */ X Xstatic int Samples, Power; Xstatic double *Real, *Imag, Maxn, magnitude(); Xstatic void scale(), fft(), max_amp(), display(), err_report(); X Xstatic int permute(); Xstatic double *Sine; Xstatic void build_trig(); X Xstatic FILE *Fpi, *Fpo; X X X/* The program */ X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X if (argc != 4) X err_report(0); X X Samples = abs(atoi(*++argv)); X Power = log10((double)Samples) / log10(2.0); X X /* Allocate memory for the variable arrays */ X X if ((Real = (double *)malloc(Chunk)) == NULL) X err_report(1); X X if ((Imag = (double *)malloc(Chunk)) == NULL) X err_report(2); X X if ((Sine = (double *)malloc(Chunk)) == NULL) X err_report(3); X X /* open the data files */ X X if ((Fpi = fopen(*++argv, "r")) == NULL) X err_report(4); X X if ((Fpo = fopen(*++argv, "w")) == NULL) X err_report(5); X X /* read in the data from the input file */ X X fread(Real, sizeof(double), Samples, Fpi); X fread(Imag, sizeof(double), Samples, Fpi); X fclose(Fpi); X X build_trig(); X scale(); X fft(); X display(); X X /* write the frequency domain data to the standard output */ X X fwrite(Real, sizeof(double), Samples, Fpo); X fwrite(Imag, sizeof(double), Samples, Fpo); X fclose(Fpo); X X /* free up memory and let's get back to our favorite shell */ X X free((char *)Real); X free((char *)Imag); X free((char *)Sine); X X exit(0); X} X X Xstatic void err_report(n) Xint n; X{ X switch (n) { X case 0: X fprintf(stderr, "Usage: fft samples in_file out_file\n"); X fprintf(stderr, "Where samples is a power of two\n"); X break; X case 1: X fprintf(stderr, "fft: Out of memory getting real space\n"); X break; X case 2: X fprintf(stderr, "fft: Out of memory getting imag space\n"); X free((char *)Real); X break; X case 3: X fprintf(stderr, "fft: Out of memory getting sine space\n"); X free((char *)Real); X free((char *)Imag); X break; X case 4: X fprintf(stderr,"fft: Unable to open data input file\n"); X free((char *)Real); X free((char *)Imag); X free((char *)Sine); X break; X case 5: X fprintf(stderr,"fft: Unable to open data output file\n"); X fclose(Fpi); X free((char *)Real); X free((char *)Imag); X free((char *)Sine); X break; X } X X exit(1); X} X X Xstatic void scale() X{ X register int loop; X X for (loop = 0; loop < Samples; loop++) { X Real[loop] /= Samples; X Imag[loop] /= Samples; X } X} X X Xstatic void fft() X{ X register int loop, loop1, loop2; X unsigned i1; /* going to right shift this */ X int i2, i3, i4, y; X double a1, a2, b1, b2, z1, z2; X X i1 = Samples >> 1; X i2 = 1; X X /* perform the butterfly's */ X X for (loop = 0; loop < Power; loop++) { X i3 = 0; X i4 = i1; X X for (loop1 = 0; loop1 < i2; loop1++) { X y = permute(i3 / (int)i1); X z1 = cosine(y); X z2 = -sine(y); X X for (loop2 = i3; loop2 < i4; loop2++) { X X a1 = Real[loop2]; X a2 = Imag[loop2]; X X b1 = z1*Real[loop2+i1] - z2*Imag[loop2+i1]; X b2 = z2*Real[loop2+i1] + z1*Imag[loop2+i1]; X X Real[loop2] = a1 + b1; X Imag[loop2] = a2 + b2; X X Real[loop2+i1] = a1 - b1; X Imag[loop2+i1] = a2 - b2; X } X X i3 += (i1 << 1); X i4 += (i1 << 1); X } X X i1 >>= 1; X i2 <<= 1; X } X} X X/* X * Display the frequency domain. X * X * The filters are aranged so that DC is in the middle filter. X * Thus -Doppler is on the left, +Doppler on the right. X */ X Xstatic void display() X{ X register int loop, offset; X int n, x; X X max_amp(); X n = Samples >> 1; X X for (loop = n; loop < Samples; loop++) { X x = (int)(magnitude(loop) * 59.0 / Maxn); X printf("%d\t|", loop - n); X offset = 0; X while (++offset <= x) X putchar('='); X X putchar('\n'); X } X X for (loop = 0; loop < n; loop++) { X x = (int)(magnitude(loop) * 59.0 / Maxn); X printf("%d\t|", loop + n); X offset = 0; X while (++offset <= x) X putchar('='); X X putchar('\n'); X } X} X X/* X * Find maximum amplitude X */ X Xstatic void max_amp() X{ X register int loop; X double mag; X X Maxn = 0.0; X for (loop = 0; loop < Samples; loop++) { X if ((mag = magnitude(loop)) > Maxn) X Maxn = mag; X } X} X X/* X * Calculate Power Magnitude X */ X Xstatic double magnitude(n) Xint n; X{ X n = permute(n); X return hypot(Real[n], Imag[n]); X} X X/* X * Bit reverse the number X * X * Change 11100000b to 00000111b or vice-versa X */ X Xstatic int permute(index) Xint index; X{ X register int loop; X unsigned n1; X int result; X X n1 = Samples; X result = 0; X X for (loop = 0; loop < Power; loop++) { X n1 >>= 1; X if (index < n1) X continue; X X /* Unix has a truncation problem with pow() */ X X result += (int)(pow(2.0, (double)loop) + .05); X index -= n1; X } X X return result; X} X X/* X * Pre-compute the sine and cosine tables X */ X Xstatic void build_trig() X{ X register int loop; X double angle, increment; X X angle = 0.0; X increment = TWO_PI / (double)Samples; X X for (loop = 0; loop < Samples; loop++) { X Sine[loop] = sin(angle); X angle += increment; X } X} X X/* EOF */ END_OF_ufft.c if test 5897 -ne `wc -c <ufft.c`; then echo shar: \"ufft.c\" unpacked with wrong size! fi # end of overwriting check fi if test -f usine.c -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"usine.c\" else echo shar: Extracting \"usine.c\" $1763 characters$ sed "s/^X//" >usine.c <<'END_OF_usine.c' X/* X * sine.c X * X * Unix Version 2.4 by Steve Sampson, Public Domain, September 1988 X * X * This program is used to generate time domain sinewave data X * for fft.c. If you want an opening target - negate the X * test frequency. X */ X X#include <stdio.h> X#include <malloc.h> X#include <math.h> X X#define TWO_PI (2.0 * 3.14159265358979324) X#define Chunk (Samples * sizeof(double)) X Xstatic double Sample, Freq, Time, *Real, *Imag; Xstatic int Samples; Xstatic void err_report(); Xstatic FILE *Fp; X X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X register int loop; X X if (argc != 3) X err_report(0); X X Samples = abs(atoi(*++argv)); X X if ((Real = (double *)malloc(Chunk)) == NULL) X err_report(1); X X if ((Imag = (double *)malloc(Chunk)) == NULL) X err_report(2); X X printf("Input The Test Frequency (Hz) ? "); X scanf("%lf", &Freq); X printf("Input The Sampling Frequency (Hz) ? "); X scanf("%lf", &Sample); X Sample = 1.0 / Sample; X X Time = 0.0; X for (loop = 0; loop < Samples; loop++) { X Real[loop] = sin(TWO_PI * Freq * Time); X Imag[loop] = -cos(TWO_PI * Freq * Time); X Time += Sample; X } X X if ((Fp = fopen(*++argv, "w")) == NULL) X err_report(3); X X fwrite(Real, sizeof(double), Samples, Fp); X fwrite(Imag, sizeof(double), Samples, Fp); X fclose(Fp); X X putchar('\n'); X X free((char *)Real); X free((char *)Imag); X X exit(0); X} X X Xstatic void err_report(n) Xint n; X{ X switch (n) { X case 0: X fprintf(stderr,"Usage: sine samples output_file\n"); X fprintf(stderr,"Where samples is a power of two\n"); X break; X case 1: X fprintf(stderr,"sine: Out of memory getting real space\n"); X break; X case 2: X fprintf(stderr,"sine: Out of memory getting imag space\n"); X free((char *)Real); X break; X case 3: X fprintf(stderr,"sine: Unable to create write file\n"); X } X X exit(1); X} X X/* EOF */ END_OF_usine.c if test 1763 -ne `wc -c <usine.c`; then echo shar: \"usine.c\" unpacked with wrong size! fi # end of overwriting check fi if test -f upulse.c -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"upulse.c\" else echo shar: Extracting \"upulse.c\" $1613 characters$ sed "s/^X//" >upulse.c <<'END_OF_upulse.c' X/* X * pulse.c X * X * Unix Version 2.4 by Steve Sampson, Public Domain, September 1988 X * X * This program is used to generate time domain pulse data for fft.c. X */ X X#include <stdio.h> X#include <malloc.h> X X#define Chunk (Samples * sizeof(double)) X Xstatic double Sample, Pw, Time, *Real, *Imag; Xstatic int Samples; Xstatic void err_report(); Xstatic FILE *Fp; X X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X register int loop; X X if (argc != 3) X err_report(0); X X Samples = abs(atoi(*++argv)); X X if ((Real = (double *)malloc(Chunk)) == NULL) X err_report(1); X X if ((Imag = (double *)malloc(Chunk)) == NULL) X err_report(2); X X printf("Input The Pulse Width (Seconds) ? "); X scanf("%lf", &Pw); X printf("Input The Sampling Time (Seconds) ? "); X scanf("%lf", &Sample); X X Time = 0.0; X for (loop = 0; loop < Samples; loop++) { X if (Time < Pw) X Real[loop] = 1.0; X else X Real[loop] = 0.0; X X Imag[loop] = 0.0; X Time += Sample; X } X X if ((Fp = fopen(*++argv, "w")) == NULL) X err_report(3); X X fwrite(Real, sizeof(double), Samples, Fp); X fwrite(Imag, sizeof(double), Samples, Fp); X fclose(Fp); X X putchar('\n'); X X free((char *)Real); X free((char *)Imag); X X exit(0); X} X X Xstatic void err_report(n) Xint n; X{ X switch (n) { X case 0: X fprintf(stderr,"Usage: pulse samples output_file\n"); X fprintf(stderr,"Where samples is a power of two\n"); X break; X case 1: X fprintf(stderr,"pulse: Out of memory getting real space\n"); X break; X case 2: X fprintf(stderr,"pulse: Out of memory getting imag space\n"); X free((char *)Real); X break; X case 3: X fprintf(stderr,"pulse: Unable to create write file\n"); X } X X exit(1); X} X X/* EOF */ END_OF_upulse.c if test 1613 -ne `wc -c <upulse.c`; then echo shar: \"upulse.c\" unpacked with wrong size! fi # end of overwriting check fi if test -f read.doc -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"read.doc\" else echo shar: Extracting \"read.doc\" $5665 characters$ sed "s/^X//" >read.doc <<'END_OF_read.doc' X X X FFT Program for Turbo-C X X By X X Steve Sampson X X X Version 2.6, November 1988 X X X The FFT program and test signal generators included in the archive, can be Xused to perform signal analysis in the frequency domain, using samples in the Xtime domain. Historically the applications have ranged from music to radar. XI recently have been doing some research in the radar field using this FFT to Xperform relative velocity measurements. This program can be further refined Xto meet your needs. X X I initially uploaded a fairly basic program, and through feedback have Xmade some improvements. It's pretty complete now as far as a start for making Xyour own specific version. Earlier Unix compatability has been removed in Xfavor of IBM graphics adaptors. X X This program uses graphics to present a 256 filter window. My current Xapplication required complex data and resulted in 256 complex points. You may Xuse this with real data which results in 128 significant points. If you feed Xthe FFT real data only (Imaginary data set to zero), then the output will be a Xmirror image, and you can ignore the left side. X XSome papers I found on the subject of FFTs are included at the end. There are Xseveral books devoted to the subject also. X XFor an example try: X X sine in X 1000 X 3000 X XWhich will sample the 1 Khz data every 333 microseconds (1 / 3 Khz). Note: The Xsample frequency should be greater than 2 times the input frequency (Nyquist Xand all that...). X XThen run fft like so: X X fft in X XAnd you should see a graphics display which has filters 0 through 255 on the XX axis, and power on the Y axis. All DC power in the time samples will appear Xin filter 128. The fundamental display frequency is shown in filter 129 and Xcan be computed as follows: X X T = Time Increment Between Samples X N = Number Of Samples (256) X Tp = N * T X X Then F = 1 / Tp X X In the example above, the time increment between samples is X 1 / 3000 or 333 microseconds. N = 256, so Tp = 85 milliseconds X and 1 / .005333 is 11.7 Hz per filter. X X Therefore each filter is a multiple of 11.7 Hertz. X X Filter 126 -23.4 Hz X Filter 127 -11.7 Hz X Filter 128 DC 00.0 Hz X Filter 129 Fundamental 11.7 Hz X Filter 130 Second Harmonic 23.4 Hz X XIn this case you'll find the sampling interval didn't work very well for the Xinput frequency. The display shows the power spread out over several filters. XThis is a limitation of the Discrete Fourier Transform in representing a Xcontinuous signal. X XA better sample rate for 1000 Hz would be 4000 Hz, in which case T = 250 ms, XTp = 64 milliseconds, and F = 15.625 Hz. 1000 / 15.625 = 64. All of the Xpower should then appear in filter 192 (128 + 64) in this example. X XRun it and see using: X X sine in X 1000 X 4000 X X fft in X XBy using negative frequencies in 'sine' you can generate opening targets: X X sine in X -1000 X 3000 X X fft in X XYou can see in these examples where weighting functions would be useful. For Xinstance using weighting you could greatly attenuate the sidebands. Usually Xthe main lobe becomes a little wider however. The current version does not Xperform any weighting. X XFor generating pulses, a second program 'pulse' is provided. It pre-loads Ximaginary data with zeros. For example: X X pulse in X .000006 X .0000008 X X fft in X XSimulates a radar with a 6 microsecond pulse and 800 nanosecond range gates, Xand produces the typical pulse spectrum display. X X XHow to run the FFT program X-------------------------- X XThe program auto-detects CGA, EGA, and VGA graphics adaptors. The CGA version Xcursor line is hard to see though. VGA is untested. X XWhen the program is run, a ruler line is drawn with tick marks every 5 filters. XAlso "Computing FFT" is displayed at the top of the screen. The program is now Xbusy computing the FFT and will not do anything useful until finished. X XWhen the FFT computation is complete, the power plot will be performed and a Xverticle cursor line will appear. "Computing FFT" will be replaced with X"Filter # 128". The program is then ready for input of commands. X XThe commands are: X X ESC Escapes back to MS-DOS X LEFT Moves the cursor left X RIGHT Moves the cursor right X CTRLLEFT Moves the cursor 10 filters left X CTRLRIGHT Moves the cursor 10 filters right X HOME Moves the cursor to filter 0 X END Moves the cursor to filter 255 X F1 Prints the display on an Epson/IBM Compatable printer X XA bee-bop when F1 is pressed means the printer has a problem (paper, power...). X X XHow to compile the programs X--------------------------- X XI use a Hard Disk configured per the Borland Manuals. If this is your settup Xalso; do this: X X 1. Copy the archive contents into any directory. X 2. Run \TURBOC\MAKE X 3. Three programs are made: fft.exe, sine.exe, and pulse.exe. X XIf you have any other system specific changes, then consult the makefile and XBorland manuals. I can also be reached at the address below. X X XFFT References X-------------- X X1. Fast Fourier Transforms On Your Home Computer, William D. Stanley, X Steven J. Peterson, BYTE Magazine, December 1978. Basic idea comes X from this program. X X2. 8052 Microcomputer simplifies FFT Design, Arnold Rosenberg, X Electronics, May 5, 1983. Used a bit reverse table based on the X routine in this program. X X3. A Fast Fourier Transform for the 8080, Robert D. Fusfeld, X Dr. Dobbs, Number 44. Gave me some ideas. X X4. A Guided Tour of the Fast Fourier Transform, G. D. Bergland, X IEEE Spectrum, July 1969. X X5. FFT - Shortcut to Fourier Analysis, Richard Klahn, Richard R. Shively, X Electronics, April 15 1968. X X--- XAll programs are entered into the Public Domain, No rights reserved. XSteve Sampson, Box 45668, Tinker AFB, OK, 73145 X75136,626 END_OF_read.doc if test 5665 -ne `wc -c <read.doc`; then echo shar: \"read.doc\" unpacked with wrong size! fi # end of overwriting check fi if test -f makefile -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"makefile\" else echo shar: Extracting \"makefile\" $874 characters$ sed "s/^X//" >makefile <<'END_OF_makefile' X# X# Makefile for FFT small model X# X# Version 2.6, Public Domain by Steve Sampson, November 1988 X# X# See readme.doc file X# X XMOD = s XLIB = \turboc\lib XINC = \turboc\include X Xall: X make fft X make sine X make pulse X make clean X Xfft: fft.obj X tcc -c -f -m$(MOD) -L$(LIB) -I$(INC) fft.c X \turboc\bgiobj \turboc\examples\cga X \turboc\bgiobj \turboc\examples\egavga X tlink /x /c $(LIB)\C0$(MOD) fft cga egavga, fft,, \ X $(LIB)\emu $(LIB)\math$(MOD) $(LIB)\C$(MOD) $(LIB)\graphics X Xsine: sine.obj X tcc -c -f -m$(MOD) -L$(LIB) -I$(INC) sine.c X tlink /x /c $(LIB)\C0$(MOD) sine, sine,, \ X $(LIB)\emu $(LIB)\math$(MOD) $(LIB)\C$(MOD) X Xpulse: pulse.obj X tcc -c -f -m$(MOD) -L$(LIB) -I$(INC) pulse.c X tlink /x /c $(LIB)\C0$(MOD) pulse, pulse,, \ X $(LIB)\emu $(LIB)\math$(MOD) $(LIB)\C$(MOD) X Xclean: X del *.obj X X# X# Dependancies X# X Xfft.obj: fft.c Xsine.obj: sine.c Xpulse.obj: pulse.c X X# EOF END_OF_makefile if test 874 -ne `wc -c <makefile`; then echo shar: \"makefile\" unpacked with wrong size! fi # end of overwriting check fi if test -f fft.c -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"fft.c\" else echo shar: Extracting \"fft.c\" $15702 characters$ sed "s/^X//" >fft.c <<'END_OF_fft.c' X/* X * fft.c X * X * Version 2.6 by Steve Sampson, Public Domain, November 1988 X * X * This program produces a Frequency Domain display from the Time Domain X * data input; using the Fast Fourier Transform. X * X * Runs in @ 30 seconds on a 5 MHz PC XT Clone without 8087. X * X * The Real data is generated by the in-phase (I) channel, and the X * Imaginary data is produced by the quadrature-phase (Q) channel of X * a Doppler Radar receiver. The middle filter is zero Hz. Closing X * targets are displayed to the right, and Opening targets to the left. X * X * Note: With Imaginary data set to zero the output is a mirror image. X * X * Usage: fft input X * 1. samples is 256. X * 2. Input is (samples * sizeof(double)) characters. X * 3. Standard error is help or debugging messages. X * X * Auto detects CGA, EGA, and VGA in Turbo-C graphics mode. X * X * See also: readme.doc, pulse.c, and sine.c. X */ X X/* Includes */ X X#include <stdio.h> X#include <alloc.h> X#include <math.h> X X#include <conio.h> X#include <dos.h> X#include <bios.h> X#include <graphics.h> X#include <string.h> X#include <stdlib.h> X X/* Defines */ X X#define SAMPLES 256 X#define POWER 8 /* 2 to the 8th power is 256 */ X X#define ESC 27 /* exit program */ X#define LEFTKEY 331 /* cursor left */ X#define RIGHTKEY 333 /* cursor right */ X#define CTRLLEFTKEY 371 /* cursor 10 left */ X#define CTRLRIGHTKEY 372 /* cursor 10 right */ X#define HOMEKEY 327 /* cursor at filter 0 */ X#define ENDKEY 335 /* cursor at filter 255 */ X#define F1 315 /* print screen */ X X#define TOP 50 X#define LEFT 64 X#define RIGHT 575 X#define MIDDLE 192 X#define TEXTX 160 X#define TEXTXN 304 X#define TEXTY 25 X X/* X * Fixed constants should be used in the macros for speed. X * X * A cosine wave leads a sine wave by 90 degrees, so offset into X * the lookup table 1/4 the way into it (256 / 4 = 64). Using X * modulo 256, lookups will wrap around to zero for numbers greater X * than 255. (cosine(200) = Sine[264 % 256] = Sine[8]). X */ X X#define permute(x) Br_table[(x)] X#define sine(x) Sine[(x)] X#define cosine(x) Sine[((x) + 64) % 256] X X/* Globals, Forward declarations */ X Xdouble Real[SAMPLES], Imag[SAMPLES], Maxn, magnitude(); Xint GraphDriver = DETECT, GraphMode, Primary, Cursor, Length; Xint Bottom, Left, PrintChar(), PrintScreen(), getkey(); Xvoid *Save, quit(), beep(), build_window(), commands(), bee_bop(); Xvoid scale(), fft(), max_amp(), display(); X X/* X * Bit Reverse Table for size 256 X * Lookup saves 20 seconds in Turbo-C over the pow() function. X * X * Br_table[x] = x inverted (eg. 00000001 flipped to 10000000) X */ X Xunsigned char Br_table[] = { X 0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0, 0x10, 0x90, 0x50, 0xd0, X 0x30, 0xb0, 0x70, 0xf0, 0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8, X 0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8, 0x04, 0x84, 0x44, 0xc4, X 0x24, 0xa4, 0x64, 0xe4, 0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4, X 0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec, 0x1c, 0x9c, 0x5c, 0xdc, X 0x3c, 0xbc, 0x7c, 0xfc, 0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2, X 0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2, 0x0a, 0x8a, 0x4a, 0xca, X 0x2a, 0xaa, 0x6a, 0xea, 0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa, X 0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6, 0x16, 0x96, 0x56, 0xd6, X 0x36, 0xb6, 0x76, 0xf6, 0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee, X 0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe, 0x01, 0x81, 0x41, 0xc1, X 0x21, 0xa1, 0x61, 0xe1, 0x11, 0x91, 0x51, 0xd1, 0x31, 0xb1, 0x71, 0xf1, X 0x09, 0x89, 0x49, 0xc9, 0x29, 0xa9, 0x69, 0xe9, 0x19, 0x99, 0x59, 0xd9, X 0x39, 0xb9, 0x79, 0xf9, 0x05, 0x85, 0x45, 0xc5, 0x25, 0xa5, 0x65, 0xe5, X 0x15, 0x95, 0x55, 0xd5, 0x35, 0xb5, 0x75, 0xf5, 0x0d, 0x8d, 0x4d, 0xcd, X 0x2d, 0xad, 0x6d, 0xed, 0x1d, 0x9d, 0x5d, 0xdd, 0x3d, 0xbd, 0x7d, 0xfd, X 0x03, 0x83, 0x43, 0xc3, 0x23, 0xa3, 0x63, 0xe3, 0x13, 0x93, 0x53, 0xd3, X 0x33, 0xb3, 0x73, 0xf3, 0x0b, 0x8b, 0x4b, 0xcb, 0x2b, 0xab, 0x6b, 0xeb, X 0x1b, 0x9b, 0x5b, 0xdb, 0x3b, 0xbb, 0x7b, 0xfb, 0x07, 0x87, 0x47, 0xc7, X 0x27, 0xa7, 0x67, 0xe7, 0x17, 0x97, 0x57, 0xd7, 0x37, 0xb7, 0x77, 0xf7, X 0x0f, 0x8f, 0x4f, 0xcf, 0x2f, 0xaf, 0x6f, 0xef, 0x1f, 0x9f, 0x5f, 0xdf, X 0x3f, 0xbf, 0x7f, 0xff X}; X X/* X * Sine/Cosine Table for size 256, Lookup saves 9 seconds in Turbo-C X * X * Sine[n] = sin(x), x = x + (2Pi / 256) X */ X Xfloat Sine[] = { X 0.000000, 0.024541, 0.049068, 0.073565, 0.098017, 0.122411, X 0.146730, 0.170962, 0.195090, 0.219101, 0.242980, 0.266713, X 0.290285, 0.313682, 0.336890, 0.359895, 0.382683, 0.405241, X 0.427555, 0.449611, 0.471397, 0.492898, 0.514103, 0.534998, X 0.555570, 0.575808, 0.595699, 0.615232, 0.634393, 0.653173, X 0.671559, 0.689541, 0.707107, 0.724247, 0.740951, 0.757209, X 0.773010, 0.788346, 0.803208, 0.817585, 0.831470, 0.844854, X 0.857729, 0.870087, 0.881921, 0.893224, 0.903989, 0.914210, X 0.923880, 0.932993, 0.941544, 0.949528, 0.956940, 0.963776, X 0.970031, 0.975702, 0.980785, 0.985278, 0.989177, 0.992480, X 0.995185, 0.997290, 0.998795, 0.999699, 1.000000, 0.999699, X 0.998795, 0.997290, 0.995185, 0.992480, 0.989177, 0.985278, X 0.980785, 0.975702, 0.970031, 0.963776, 0.956940, 0.949528, X 0.941544, 0.932993, 0.923880, 0.914210, 0.903989, 0.893224, X 0.881921, 0.870087, 0.857729, 0.844854, 0.831470, 0.817585, X 0.803208, 0.788346, 0.773010, 0.757209, 0.740951, 0.724247, X 0.707107, 0.689541, 0.671559, 0.653173, 0.634393, 0.615232, X 0.595699, 0.575808, 0.555570, 0.534998, 0.514103, 0.492898, X 0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895, X 0.336890, 0.313682, 0.290285, 0.266713, 0.242980, 0.219101, X 0.195090, 0.170962, 0.146730, 0.122411, 0.098017, 0.073565, X 0.049068, 0.024541, 0.000000, -0.024541, -0.049068, -0.073565, X -0.098017, -0.122411, -0.146730, -0.170962, -0.195090, -0.219101, X -0.242980, -0.266713, -0.290285, -0.313682, -0.336890, -0.359895, X -0.382683, -0.405241, -0.427555, -0.449611, -0.471397, -0.492898, X -0.514103, -0.534998, -0.555570, -0.575808, -0.595699, -0.615232, X -0.634393, -0.653173, -0.671559, -0.689541, -0.707107, -0.724247, X -0.740951, -0.757209, -0.773010, -0.788346, -0.803208, -0.817585, X -0.831470, -0.844854, -0.857729, -0.870087, -0.881921, -0.893224, X -0.903989, -0.914210, -0.923880, -0.932993, -0.941544, -0.949528, X -0.956940, -0.963776, -0.970031, -0.975702, -0.980785, -0.985278, X -0.989177, -0.992480, -0.995185, -0.997290, -0.998795, -0.999699, X -1.000000, -0.999699, -0.998795, -0.997290, -0.995185, -0.992480, X -0.989177, -0.985278, -0.980785, -0.975702, -0.970031, -0.963776, X -0.956940, -0.949528, -0.941544, -0.932993, -0.923880, -0.914210, X -0.903989, -0.893224, -0.881921, -0.870087, -0.857729, -0.844854, X -0.831470, -0.817585, -0.803208, -0.788346, -0.773010, -0.757209, X -0.740951, -0.724247, -0.707107, -0.689541, -0.671559, -0.653173, X -0.634393, -0.615232, -0.595699, -0.575808, -0.555570, -0.534998, X -0.514103, -0.492898, -0.471397, -0.449611, -0.427555, -0.405241, X -0.382683, -0.359895, -0.336890, -0.313682, -0.290285, -0.266713, X -0.242980, -0.219101, -0.195090, -0.170962, -0.146730, -0.122411, X -0.098017, -0.073565, -0.049068, -0.024541 X}; X XFILE *Fpi, *Fpo; X X X/* The program */ X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X if (argc != 2) { X fprintf(stderr, "Usage: fft input_file\n"); X exit(1); X } X X setcbrk(1); /* Allow Control-C checking */ X ctrlbrk(quit); /* Execute quit() if Control-C detected */ X X /* open the data file */ X X if ((Fpi = fopen(*++argv, "rb")) == NULL) { X fprintf(stderr,"fft: Unable to open data input file\n"); X exit(1); X } X X /* read in the data */ X X fread(Real, sizeof(double), SAMPLES, Fpi); X fread(Imag, sizeof(double), SAMPLES, Fpi); X fclose(Fpi); X X build_window(); X scale(); X fft(); X display(); X X /* wait for keyboard commands */ X X commands(); X} X X Xvoid scale() X{ X register int loop; X X for (loop = 0; loop < SAMPLES; loop++) { X Real[loop] /= SAMPLES; X Imag[loop] /= SAMPLES; X } X} X X Xvoid fft() X{ X register int loop, loop1, loop2; X unsigned i1; /* going to right shift this */ X int i2, i3, i4, y; X double a1, a2, b1, b2, z1, z2; X X i1 = SAMPLES >> 1; X i2 = 1; X X /* perform the butterfly's */ X X for (loop = 0; loop < POWER; loop++) { X i3 = 0; X i4 = i1; X X for (loop1 = 0; loop1 < i2; loop1++) { X y = permute(i3 / (int)i1); X z1 = cosine(y); X z2 = -sine(y); X X for (loop2 = i3; loop2 < i4; loop2++) { X X a1 = Real[loop2]; X a2 = Imag[loop2]; X X b1 = z1*Real[loop2+i1] - z2*Imag[loop2+i1]; X b2 = z2*Real[loop2+i1] + z1*Imag[loop2+i1]; X X Real[loop2] = a1 + b1; X Imag[loop2] = a2 + b2; X X Real[loop2+i1] = a1 - b1; X Imag[loop2+i1] = a2 - b2; X } X X i3 += (i1 << 1); X i4 += (i1 << 1); X } X X i1 >>= 1; X i2 <<= 1; X } X} X X/* X * Display the frequency domain. X * X * The filters are aranged so that DC is in the middle filter. X * Thus -Doppler is on the left, +Doppler on the right. X */ X Xvoid display() X{ X register int loop, offset; X int n, x; X X n = SAMPLES >> 1; X max_amp(); X X /* X * Graphics screen horizontal configuration: X * X * | 64 pixels | 512 pixels | 64 pixels | X * | margin | picture | margin | X * X * Spectral lines are two bits wide X */ X X for (loop = n, offset = LEFT; loop < SAMPLES; loop++, offset++) { X x = (int)(magnitude(loop) * Length / Maxn); X bar((offset + loop - n), Bottom - x, (offset + loop - n) + 1, Bottom); X } X X for (loop = 0, offset = MIDDLE; loop < n; loop++, offset++) { X x = (int)(magnitude(loop) * Length / Maxn); X bar((offset + loop + n), Bottom - x, (offset + loop + n) + 1, Bottom); X } X} X X/* X * Find maximum amplitude X */ X Xvoid max_amp() X{ X register int loop; X double mag; X X Maxn = 0.0; X for (loop = 0; loop < SAMPLES; loop++) { X if ((mag = magnitude(loop)) > Maxn) X Maxn = mag; X } X} X X/* X * Calculate Power Magnitude X */ X Xdouble magnitude(n) Xint n; X{ X n = permute(n); X return hypot(Real[n], Imag[n]); X} X Xvoid build_window() X{ X unsigned i_size; X X /* settup graphics mode */ X X registerbgidriver(CGA_driver); X registerbgidriver(EGAVGA_driver); X X initgraph(&GraphDriver, &GraphMode, ""); X X if (graphresult() != grOk) { X fprintf(stderr, "fft: %s\n", grapherrormsg(graphresult())); X exit(1); X } X X if ((i_size = getgraphmode()) != CGAHI) X setbkcolor(BLACK); X X switch (i_size) { X case CGAHI: X Primary = WHITE; X Cursor = WHITE; X Length = 100.0; X Bottom = 150; X break; X case EGAHI: X Primary = LIGHTGRAY; X Cursor = LIGHTGREEN; X Length = 250.0; X Bottom = 300; X break; X case VGAHI: X Primary = LIGHTGRAY; X Cursor = LIGHTGREEN; X Length = 380.0; X Bottom = 430; X } X X setcolor(Primary); X setfillstyle(SOLID_FILL, Primary); X settextjustify(LEFT_TEXT, TOP_TEXT); X settextstyle(DEFAULT_FONT, HORIZ_DIR, 2); X X cleardevice(); X X outtextxy(TEXTX, TEXTY, "Computing FFT"); X X /* Draw ruler line */ X X line(LEFT, Bottom + 5, RIGHT, Bottom + 5); X for (i_size = LEFT; i_size <= RIGHT ; i_size += 10) { X line(i_size, Bottom + 5, i_size, Bottom + 10); X line(i_size + 1, Bottom + 5, i_size + 1, Bottom + 10); X } X X /* create under-cursor memory */ X X i_size = imagesize(LEFT, TOP, LEFT + 1, Bottom); X Save = malloc(i_size); X} X Xvoid commands() X{ X int key, filter; X char string[4]; X X setcolor(BLACK); /* erase text */ X outtextxy(TEXTX, TEXTY, "Computing FFT"); X setcolor(Primary); X X Left = 320; X filter = 128; X strcpy(string, "128"); X outtextxy(TEXTX, TEXTY, "Filter # 128"); X getimage(Left, TOP, Left + 1, Bottom, Save); X setfillstyle(SOLID_FILL, Cursor); X bar(Left, TOP, Left + 1, Bottom); X X for (;;) { X while (bioskey(1) == 0) /* wait for key press */ X ; X key = getkey(); X X switch(key) { /* find out which key was pressed */ X case HOMEKEY: Xj1: putimage(Left, TOP, Save, COPY_PUT); X filter = 0; X Left = LEFT; X getimage(Left, TOP, Left + 1, Bottom, Save); X bar(Left, TOP, Left + 1, Bottom); X X setcolor(BLACK); /* erase old text */ X outtextxy(TEXTXN, TEXTY, string); X setcolor(Primary); X X sprintf(string, "%d", filter); X outtextxy(TEXTXN, TEXTY, string); X break; X case ENDKEY: Xj2: putimage(Left, TOP, Save, COPY_PUT); X filter = 255; X Left = RIGHT - 1; X getimage(Left, TOP, Left + 1, Bottom, Save); X bar(Left, TOP, Left + 1, Bottom); X X setcolor(BLACK); /* erase old text */ X outtextxy(TEXTXN, TEXTY, string); X setcolor(Primary); X X sprintf(string, "%d", filter); X outtextxy(TEXTXN, TEXTY, string); X break; X case CTRLLEFTKEY: X if (filter <= 9) X goto j1; X X putimage(Left, TOP, Save, COPY_PUT); X Left -= 20; X filter -= 10; X getimage(Left, TOP, Left + 1, Bottom, Save); X bar(Left, TOP, Left + 1, Bottom); X X setcolor(BLACK); /* erase old text */ X outtextxy(TEXTXN, TEXTY, string); X setcolor(Primary); X X sprintf(string, "%d", filter); X outtextxy(TEXTXN, TEXTY, string); X break; X case LEFTKEY: X if (filter <= 0) { X beep(); X break; X } X X putimage(Left, TOP, Save, COPY_PUT); X Left -= 2; X filter--; X getimage(Left, TOP, Left + 1, Bottom, Save); X bar(Left, TOP, Left + 1, Bottom); X X setcolor(BLACK); /* erase old text */ X outtextxy(TEXTXN, TEXTY, string); X setcolor(Primary); X X sprintf(string, "%d", filter); X outtextxy(TEXTXN, TEXTY, string); X break; X case CTRLRIGHTKEY: X if (filter >= 245) X goto j2; X X putimage(Left, TOP, Save, COPY_PUT); X Left += 20; X filter += 10; X getimage(Left, TOP, Left + 1, Bottom, Save); X bar(Left, TOP, Left + 1, Bottom); X X setcolor(BLACK); /* erase old text */ X outtextxy(TEXTXN, TEXTY, string); X setcolor(Primary); X X sprintf(string, "%d", filter); X outtextxy(TEXTXN, TEXTY, string); X break; X case RIGHTKEY: X if (filter >= 255) { X beep(); X break; X } X X putimage(Left, TOP, Save, COPY_PUT); X Left += 2; X filter++; X getimage(Left, TOP, Left + 1, Bottom, Save); X bar(Left, TOP, Left + 1, Bottom); X X setcolor(BLACK); /* erase old text */ X outtextxy(TEXTXN, TEXTY, string); X setcolor(Primary); X X sprintf(string, "%d", filter); X outtextxy(TEXTXN, TEXTY, string); X break; X case F1: X if (PrintScreen()) X bee_bop(); X break; X case ESC: X quit(); X default: X beep(); X } X } X} X Xvoid beep() X{ X sound(200); X delay(125); X nosound(); X} X Xvoid bee_bop() X{ X sound(1000); X delay(80); X sound(200); X delay(160); X sound(1000); X nosound(); X} X X/* X * Return to MS-DOS operating system X */ X Xvoid quit() X{ X sound(1000); X delay(250); X nosound(); X X closegraph(); X X free((char *)Save); X X exit(0); X} X X/* X * Dumps a screen to printer in Double Density landscape format. X * This routine is from the Borland Programming Forum on CompuServe. X * X * returns TRUE on error X */ X Xint PrintScreen() X{ X int X,Y,YOfs; X int BitData,MaxBits; X unsigned int Height,Width; X struct viewporttype Vport; X X getviewsettings(&Vport); X Width = Vport.bottom - Vport.top; X Height = (Vport.right + 1) - Vport.left; X X if (PrintChar(27)) X return (1); X X if (PrintChar('3')) X return (1); X X if (PrintChar(24)) /* 24/216 inch line spacing */ X return (1); X X Y = 0; X while (Y < Height) { X if (PrintChar(27)) X return (1); X X if (PrintChar('L')) /* Double Density */ X return (1); X X if (PrintChar(Width & 0xff)) X return (1); X X if (PrintChar(Width >> 8)) X return (1); X X for (X = Width - 1; X >= 0; X--) { X BitData = 0; X if ((Y + 7) <= Height) X MaxBits = 7; X else X MaxBits = Height - Y; X X for (YOfs = 0; YOfs <= MaxBits; YOfs++) { X BitData = BitData << 1; X if (getpixel(YOfs + Y, X) > 0) X BitData++; X } X X if (PrintChar(BitData)) X return (1); X } X X if (PrintChar('\r')) X return (1); X X if (PrintChar('\n')) X return (1); X X Y += 8; X } X X PrintChar(12); /* formfeed */ X X return (0); X} X X/* X * returns TRUE on error X */ X Xint PrintChar(out) Xint out; X{ X unsigned char ret; X X ret = biosprint(0, out, 0); X if ((ret & 0x29) || !(ret & 0x10)) X return 1; X else X return 0; X} X Xint getkey() X{ X int key, lo, hi; X X key = bioskey(0); X lo = key & 0x00FF; X hi = (key & 0xFF00) >> 8; X X return((lo == 0) ? hi + 256 : lo); X} X X/* EOF */ END_OF_fft.c if test 15702 -ne `wc -c <fft.c`; then echo shar: \"fft.c\" unpacked with wrong size! fi # end of overwriting check fi if test -f sine.c -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"sine.c\" else echo shar: Extracting \"sine.c\" $1139 characters$ sed "s/^X//" >sine.c <<'END_OF_sine.c' X/* X * sine.c X * X * Version 2.6 by Steve Sampson, Public Domain, November 1988 X * X * This program is used to generate time domain sinewave data X * for fft.c. If you want an opening target - negate the X * test frequency. X */ X X#include <stdio.h> X#include <alloc.h> X#include <math.h> X X#define SAMPLES 256 X#define TWO_PI (2.0 * 3.14159265358979324) X Xdouble Sample, Freq, Time, Real[SAMPLES], Imag[SAMPLES]; XFILE *Fp; X X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X register int loop; X X if (argc != 2) { X fprintf(stderr,"Usage: sine output_file\n"); X exit(1); X } X X printf("Input The Test Frequency (Hz) ? "); X scanf("%lf", &Freq); X printf("Input The Sampling Frequency (Hz) ? "); X scanf("%lf", &Sample); X X Sample = 1.0 / Sample; X Time = 0.0; X X for (loop = 0; loop < SAMPLES; loop++) { X Real[loop] = sin(TWO_PI * Freq * Time); X Imag[loop] = -cos(TWO_PI * Freq * Time); X Time += Sample; X } X X if ((Fp = fopen(*++argv, "wb")) == NULL) { X fprintf(stderr,"sine: Unable to create write file\n"); X exit(1); X } X X fwrite(Real, sizeof(double), SAMPLES, Fp); X fwrite(Imag, sizeof(double), SAMPLES, Fp); X fclose(Fp); X X putchar('\n'); X X exit(0); X} END_OF_sine.c if test 1139 -ne `wc -c <sine.c`; then echo shar: \"sine.c\" unpacked with wrong size! fi # end of overwriting check fi if test -f pulse.c -a "${1}" != "-c" ; then echo shar: Will not over-write existing file \"pulse.c\" else echo shar: Extracting \"pulse.c\" $985 characters$ sed "s/^X//" >pulse.c <<'END_OF_pulse.c' X/* X * pulse.c X * X * Version 2.6 by Steve Sampson, Public Domain, November 1988 X * X * This program is used to generate time domain pulse data for fft.c. X */ X X#include <stdio.h> X#include <alloc.h> X X#define SAMPLES 256 X Xdouble Sample, Pw, Time, Real[SAMPLES], Imag[SAMPLES]; XFILE *Fp; X X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X register int loop; X X if (argc != 2) { X fprintf(stderr,"Usage: pulse output_file\n"); X exit(1); X } X X printf("Input The Pulse Width (Seconds) ? "); X scanf("%lf", &Pw); X printf("Input The Sampling Time (Seconds) ? "); X scanf("%lf", &Sample); X X Time = 0.0; X X for (loop = 0; loop < SAMPLES; loop++) { X if (Time < Pw) X Real[loop] = 1.0; X else X Real[loop] = 0.0; X X Imag[loop] = 0.0; X Time += Sample; X } X X if ((Fp = fopen(*++argv, "wb")) == NULL) { X fprintf(stderr,"pulse: Unable to create write file\n"); X exit(1); X } X X fwrite(Real, sizeof(double), SAMPLES, Fp); X fwrite(Imag, sizeof(double), SAMPLES, Fp); X fclose(Fp); X X putchar('\n'); X X exit(0); X} END_OF_pulse.c if test 985 -ne `wc -c <pulse.c`; then echo shar: \"pulse.c\" unpacked with wrong size! fi # end of overwriting check fi echo shar: End of shell archive. exit 0