Category : C++ Source Code
Archive : ZPP10.ZIP
Filename : FFT.CPP
Output of file : FFT.CPP contained in archive : ZPP10.ZIP
/* fft.c by Carl W. Moreland Version 2.3 04/03/92 */
/* -------------------------------------------------------------------- */
/* */
/* Fast Fourier Transform */
/* */
/* The following program is intended for Fourier analysis of circuit */
/* simulations. It can be invoked from the command line with no inter- */
/* action so it is ideal for batch files. The syntax is: */
/* */
/* fft infile.ext
/* */
/* infile.ext is the name of the file containing data. There must be */
/* 2^n data points and must represent an integer number of cycles. */
/* */
/* outfile.ext is an optional file to store the harmonic information. */
/* */
/* /c allows the user to enter comments for each analysis. This is in */
/* addition to the command line comment. Each comment line is limi- */
/* ted to 80 characters; there is no line limit. */
/* */
/* /n=x specifies the highest order harmonic to be printed. By default */
/* the program will print up to the 10th, or numpts/2-1, whatever */
/* is larger. If x is larger than numpts/2-1 it will be ignored. */
/* */
/* /e specifies exponential format instead of decimal. This is for */
/* those working on 14+ bit projects. */
/* */
/* /s suppresses screen output */
/* */
/* Version 2 now supports comments from the source file. They must */
/* begin with the ""#"" character. DC offset and fundamental amplitude */
/* are now printed. */
/* */
/* Version 2.1 adds dynamic memory allocation for the data array so max */
/* number of points is now limited by RAM. (05/04/90) */
/* */
/* Version 2.2 fixes a bug that gives incorrect results when the DC */
/* offset is greater than the fundamental amplitude. Harmonic bins are */
/* now sorted wrt the fundamental. (04/17/91) */
/* */
/* Version 2.3 uses the C++ complex class. (04/03/92) */
/* */
/* Sun compiler: cc fft.c -lm -offt (Version 2.2 or before) */
/* Note: Sun C does not recognize # conditionals or modern parameter */
/* passing. */
/* -------------------------------------------------------------------- */
#include
#include
#include
#include
#include
#include
#include "complex.h"
void get_args(int argc, char *argv[]);
int get_data(void);
void put_data(complex x[]);
void fft(int N, complex *x);
struct arg_struct /* Holds the command line arguments */
{
char input[20];
char output[20];
int num_harm;
int com_flag;
int exp_flag;
int sup_flag;
int out_flag;
char comment[81];
} args;
char *message[] =
{
"+----------------------------------------------------------------------+",
"| Fast Fourier Transform Version 2.3 04/03/92 Carl W. Moreland |",
"| |",
"| Syntax: |",
"| |",
"| fft infile.ext
"| |",
"| infile.ext is the name of the file containing data. There must be |",
"| 2^n data points and must represent an integer number of cycles. |",
"| |",
"| outfile.ext is an optional file to store the harmonic information. |",
"| |",
"| /c allows interactive comments, 80 characters max. |",
"| /n=x specifies the highest order harmonic to be printed. By default |",
"| the program will print up to the 10th or numpts/2-1. |",
"| /e specifies exponential format instead of decimal for the output. |",
"| /s suppresses screen output |",
"+----------------------------------------------------------------------+"
};
complex *x;
/* ----- main --------------------------------------------------------- */
main(int argc, char *argv[])
{
register int i;
int numpts;
if(argc < 2) /* Check for proper syntax */
{
for(i=0; i<18; i++)
printf("%s\n", message[i]);
exit(0);
}
get_args(argc, argv); /* Parse the command line arguments */
x = (complex *)malloc(8*sizeof(complex));
if(x == NULL)
{
printf("Insufficient memory for 8 data points\n");
exit(0);
}
numpts = get_data(); /* Read the input file */
if(numpts/2 < args.num_harm) /* See if there's less than 10 harmonics*/
args.num_harm = numpts/2 - 1;
fft(numpts, x); /* Perform FFT on data, return Mag & Ph */
put_data(x); /* Write the output to screen/file */
free(x);
return 0;
}
/* ----- get_args ----------------------------------------------------- */
void get_args(int argc, char *argv[])
{
register int i, j, k;
char temp[81]; /* Holds the comment string */
args.output[0] = '\0'; /* Initialize all parameters */
args.num_harm = 10;
args.com_flag = 0;
args.exp_flag = 0;
args.sup_flag = 0;
args.out_flag = 0;
args.comment[0] = '\0';
strcpy(args.input, argv[1]); /* store the input filename */
for(i=2; i
/* Check for "/" or "-" switch */
if(argv[i][0] == '/' || argv[i][0] == '-')
{
if(argv[i][1] == 'c' || argv[i][1] == 'C') args.com_flag = 1;
if(argv[i][1] == 'e' || argv[i][1] == 'E') args.exp_flag = 1;
if(argv[i][1] == 's' || argv[i][1] == 'S') args.sup_flag = 1;
if(argv[i][1] == 'n' || argv[i][1] == 'N')
{
k = 0;
for(j=2;;j++)
{
if(argv[i][j] == '\0')
{
temp[k++] = argv[i][j];
break;
}
if(isdigit(argv[i][j])) temp[k++] = argv[i][j];
}
args.num_harm = atoi(temp);
}
continue;
}
if(argv[i][0] == 0x27) /* Check for "'" for comment */
{
j = 1;
for(k=0; k<81; k++) /* up to 80 characters */
{
if(argv[i][j] == 0x27) /* Check for "'" end of comment */
{
temp[k] = '\0';
break;
}
if(argv[i][j] == '\0') /* Check for new word */
{
temp[k] = ' ';
i++;
j=0;
continue;
}
if(k == 80) /* 80 is the max so terminate */
{
temp[k] = '\0';
break;
}
temp[k] = argv[i][j++]; /* Copy character to temp string*/
}
strcpy(args.comment, temp);
continue;
}
/* If it's nothing else then it */
strcpy(args.output, argv[i]); /* must be the output file name */
args.out_flag = 1;
}
}
/* ----- get_data ----------------------------------------------------- */
int get_data()
{
register int i;
int points = 8; /* start with 8 data points */
char temp[81];
FILE *infile, *outfile;
if((infile = fopen(args.input, "r")) == NULL) /* Open input file */
{
puts(" File not found\n");
free(x);
exit(0);
}
if(args.out_flag) outfile = fopen(args.output, "a");
i = 0;
/* Read file until EOF */
while(fgets(temp, 81, infile) != NULL)
{
if(temp[0] == '#')
{
if(args.out_flag)
{
temp[0] = 0x20;
fprintf(outfile, "%s", temp);
}
continue;
}
if(i == points)
{
points *= 2;
x = (complex *)realloc(x, points*sizeof(complex));
if(x == NULL)
{
printf("Insufficient memory for %d data points\n", points);
exit(0);
}
}
x[i].re = atof(temp);
x[i++].im = 0;
}
fclose(infile);
if(args.out_flag) fclose(outfile);
return(i);
}
/* ----- put_data ----------------------------------------------------- */
void put_data(complex x[])
{
register int i;
char comment[81];
double db, bits;
double sum = 0;
FILE *outfile;
if(args.out_flag) /* Write data to an output file */
{
outfile = fopen(args.output, "a"); /* Open file for appending */
if(args.com_flag) /* Allow interactive comments */
{
printf("Enter comment(s),
for(;;)
{
gets(comment);
if(comment[0] == NULL) break;
fprintf(outfile, " %s\n", comment);
}
}
}
if(!args.sup_flag) /* Write data to the screen */
{
printf(" DC offset is %.2f\n", x[0].re);
printf(" Fundamental amplitude is %.4f\n\n", x[1].re);
printf(" # Magnitude Phase dB Bits \n");
printf(" --- --------- -------- -------- ------\n");
}
if(args.out_flag)
{
fprintf(outfile, "\n");
fprintf(outfile, " DC offset is %.2f\n", x[0].re);
fprintf(outfile, " Fundamental amplitude is %.4f\n\n", x[1].re);
fprintf(outfile, " # Magnitude Phase dB Bits \n");
fprintf(outfile, " --- --------- -------- -------- ------\n");
}
for(i=2; i<=args.num_harm; i++) /* Loop thru each harmonic */
{
if(x[i].re == 0)
{
db = -999;
bits = 999;
}
else
{
db = 20*log10(x[i].re); /* Convert to dB */
bits = -db/6.0206; /* Calculate number of bits */
}
if(args.exp_flag) /* Check for /e switch */
{
if(!args.sup_flag)
printf(" %2d %.4e %8.3f %8.3f %6.3f\n", i, x[i].re, x[i].im, db, bits);
if(args.out_flag)
fprintf(outfile, " %2d %.4e %8.3f %8.3f %6.3f\n", i, x[i].re, x[i].im, db, bits);
}
else
{
if(!args.sup_flag)
printf(" %2d %9.7f %8.3f %8.3f %6.3f\n", i, x[i].re, x[i].im, db, bits);
if(args.out_flag)
fprintf(outfile, " %2d %9.7f %8.3f %8.3f %6.3f\n", i, x[i].re, x[i].im, db, bits);
}
sum += x[i].re*x[i].re; /* Keep a running total */
}
sum = sqrt(sum); /* Calculate THD */
if(sum == 0)
{
bits = 99;
}
else
{
db = 20*log10(sum); /* Convert to dB */
bits = -db/6.0206; /* Calculate number of bits */
}
if(!args.sup_flag)
printf(" ----------------------------------------------\n");
if(args.out_flag)
fprintf(outfile, " ----------------------------------------------\n");
if(args.exp_flag)
{
if(!args.sup_flag)
printf(" THD = %.4e%% or %6.3f bits\n\n", 100*sum, bits);
if(args.out_flag)
fprintf(outfile, " THD = %.4e%% or %6.3f bits\n\n\n", 100*sum, bits);
}
else
{
if(!args.sup_flag)
printf(" THD = %7.4f%% or %6.3f bits\n\n", 100*sum, bits);
if(args.out_flag)
fprintf(outfile, " THD = %7.4f%% or %6.3f bits\n\n\n", 100*sum, bits);
}
if(args.out_flag) fclose(outfile);
return;
}
/* -------------------------------------------------------------------- */
/* */
/* >>>>>>>>>>>>>>>>> Decimation in Time FFT Algorithm <<<<<<<<<<<<<<<<< */
/* */
/* The following subroutine is taken from Digital Signal Processing by */
/* A. Oppenheim and R. Schafer, p.331. The original FORTRAN code was */
/* translated into C by Carl Moreland. Since FORTRAN arrays are indexed */
/* from 1 to N several modifications were made to accomodate C's 0 to */
/* N-1 index. Usage: */
/* */
/* void fft(N, x) */
/* */
/* integer N = # of samples */
/* complex x is a pointer to an array of complex of size N, where */
/* complex is defined as a structure of double re, double im. */
/* */
/* -------------------------------------------------------------------- */
void fft(int N, complex *x)
{
register int i, j, k;
int M, LE, LE1, ip, fund=1, num;
complex u, w, t, tm;
double PI = 3.14159265358;
M = log((double)N)/log(2.0); /* Convert numpts to log base 2 */
j = -N/2; /* Initial value of j for bit reversal */
for(i=0; i
k = N/2; /* The next few lines calculate the */
/* proper j value for bit reversing */
while(k<=j) /* x[i] and x[j]. Some modifications */
{ /* were made here to deal with C's */
j = j - k; /* array indices such as the initial */
k = k/2; /* j value above. A Polish way do it */
} /* but it works... */
j = j + k;
if(i
t = x[j]; /* porary holder. */
x[j] = x[i];
x[i] = t;
}
}
for(k=1; k<=M; k++) /* Run through M stages of computations */
{
LE = pow(2.0,(double)k);
LE1 = LE/2;
u = Z1;
w = complex(cos(PI/(double)LE1), -sin(PI/(double)LE1));
for(j=0; j
for(i=j; i
ip = i+LE1; /* x[ip] is the other term... */
/* Butterfly algorithm calculates x[i]' and x[ip]' from x[i] */
/* and x[ip] as follows: */
/* x[i]' = x[i] + x[ip]*w */
/* x[ip]' = x[i] - x[ip]*w */
/* where w = exp(j2cr/LE) = W(LE)^r */
t = x[ip] * u; /* Temporary holder */
x[ip] = x[i] - t;
x[i] += t;
}
u *= w;
}
}
Complex.SetArgMode(Z_DEGREES);
x[0].topolar();
for(i=1; i
x[i].topolar(); /* Convert to magnitude & phase */
if(x[i].re > x[fund].re) /* Find the fundamental */
fund = i;
}
x[N/2] = x[0]; /* Sort the bins. Use the upper half of */
/* of the fft array to temporarily */
for(i=1; i
num = fund*i;
while(num > N/2)
num = abs(N - num);
x[N/2+i] = x[num];
}
x[0] = x[N/2];
x[1] = x[N/2+1]; /* Now copy the newly sorted harmonics */
/* back down to the lower half of the */
for(i=2; i
x[i] = x[N/2+i]; /* magnitudes. */
x[i].re /= x[1].re;
}
x[0].re /= N;
x[1].re /= N/2;
}
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/