fft.c

傅立叶变化

/*	name...
		fft

	bugs...
		Should do transforms on several functions (-b option).
		use large model.
		eliminate use of extra few entries, permitting use of exactly
			64K arrays (4 bytes * 8K entries * 2 arrays)

	performance...
		4096 point transform in 1:13.71 with 6 MHz 80286 and 4 MHz 80287.

	bugs...
		check for unequally spaced data is too strict.


	history...
		6 Jun 91	Ver 1.42: can read data from stdin.  Uses args.c
					to read switches and echo command line.
		19 Apr 89	Ver 1.41: forcing space between columns in output.
		1 Jun 87	Checking for unequally spaced data.
		-- 1.40 --
		11 May 87	-o option installed, -p	option accepts total number
					as well as padding factor.
		8 May 87	Reading into doubles and performing range check
					before converting to floats.
		1 May 87	calculating with floats rather than doubles.
					performing 4096 point transforms.
					-a option suppresses frequencies on output
		-- 1.3 --
		29 Apr 87	Fixed padding factor calculation
					added -m option to print only freq and mag
					added -z option to shift origin
					enhanced help display
		6 Apr 87	Default edge weight .1 instead of .9,
					Default abscissa step 1 (instead of 0!).
		20 Mar 87	Conversion to/from dBs or to magnitude/phase removed.
					Default output to STDOUT.
					Removed prompts.
		8 Feb 87	Inverse FFT code split off into separate file.
		7 Nov 86	-o omits phase printout, -p specifies padding factor.
			padding factor defaults to 1/4 of maximum, # points defaults
			to all.
		6 Nov 86	Switch s sets up windowing.
		2 Nov 86	Ported to DeSmet C, input & output files in ASCII,
			options read from command line.
		12 Jul 84  Printing time between samples.
		11 Jul 84  One sided windows, degree 6, 10, or 16.
			Calculating energy before & after windowing.
		10 Jul 84  No longer printing out y values. Announcing
			program version.  One sided windowing, degree 16.
			Allows peak to be any value from 200 to 22000 without
			rescaling.  No longer discarding zero frequency point.
		26 Jun 84  Weight at window edge (epsilon) specified by
			user.  Number of input data points to be retained
			is specified by user.
		25 Jun 84  Supergaussian windowing (degree 6) is optional,
			beeping after finishing transform.
		22 Jun 84  Scaling real & imaginary parts to the range
			10900...22000.  Padding by a factor of 4, then 
			discarding 50% of the data (high frequency values) so
			512 magitude values are retained.  No "magnitude"
			values printed out (not even zeros).
		18 Jun 84  writing out first magnitude, then phase.
			Correctly calculating time step.
		15 Jun 84  setting phase to zero.
		12 Jun 84  scaling so largest element is 2048.
			Calculating magnitude & phase.
		May 84     written
*/
#include <stdio.h>
#include <math.h>

#define VERSION "1.42"



#define BIGFLOAT 6.e38	/* no floats larger than this  */
#define ENTRIES 4098
#define MAXLABELS 50
#define BUFSIZE 200
#define DEFAULT_EDGE_WEIGHT .1
#define FLOAT float		/* change to "double" to restore higher precision */

FILE *ifile=stdin, *ofile=stdout;

char buffer[128], iname[35], oname[35];
char buf[BUFSIZE];
char *label_array[MAXLABELS];
char **p_text=label_array;


int m,			/* number of input values to be retained */
n,				/* number of data points to be Fourier transformed */
*nnn,			/* points to n */
breaking=0,		/* nonzero if finding separate transforms for several functions */
labels,			/* number of labels stored */
super=0,		/* if nonzero, the degree of supergaussian window desired */
automatic_abscissas, /* nonzero if abscissas to be calculated */
abscissa_arguments=0,
shifting=0,	/* nonzero if data to be shifted */
pad_factor=0,	/* nonzero if specific padding factor was requested */
keeping=0,		/* if nonzero, the number of input values to be kept */
magnitudes=0,	/* nonzero if only magnitudes are to be written */
f_arguments=0,
x_arguments=0;

int standard_input=0;		/* nonzero if data is coming from standard input */
int index_array[MAXLABELS];	/* indexes into x and y */
int *p_data=index_array;

FLOAT *x,		/* time or frequency values */
*y;				/* voltages (on input) or magnitudes (on output) */
double abscissa,
abscissa_step=1.,
before, after,	/* energy in data before & after windowing */
origin,			/* abscissa value to be treated as time zero */
wt=0.,			/* weight on last data point kept */
output=0.,		/* if nonzero, the number of calculated values to be printed */
fmin, fmax,
xmin, xmax;		/* first & last data points to be used */

double energy(), atan2();

main(argc, argv) int argc; char **argv;
{	int i, nn, windowing;

	double amp, phase, factor, q, pi, time, dt;
	double this, big, sweep, scale, yr, yi;

	{double junk;
	sscanf("3.4","%lf",&junk);		/* workaround for DESMET sscanf bug */
	}

	argc = args(argc, argv);		/* fetch switches */

	read_data(argc, argv);			/* read input data */

/*	puts(" time function\n"); listt(y, 1, 10); */
	nn=n;
		{double dt;
		dt=x[2]-x[1];
		if(nn*fabs(x[nn]-x[nn-1]-dt) + fabs(x[nn]-x[1]-(nn-1)*dt) > .001*nn*dt)
			{printf("fft: times not equally spaced\n");
			printf("x[1]=%f, x[2]=%f, x[%d]=%f, x[%d]=%f, dt=%f",
				x[1], x[2], nn-1, x[nn-1], nn, x[nn], dt);
			exit();
			}
		}
	nn=pad();
/*	puts(" time function after padding\n"); listt(y, 1, 10); */
	before=energy();
	windowing = wind16() || wind10() || wind6();
	if(windowing)
		{after=energy();
		fprintf(stderr, "Energy loss due to windowing = %5.2f%% \n",
			100.*(before-after)/before);
		}
	else after=before;
	fprintf(stderr, "Energy (sum of squares of data%s) %s= %10.3g \n", 
	(!automatic_abscissas || abscissa_arguments) ? ", times time step" : "",
	windowing?"after windowing ":"", after);
/*	puts(" time function\n"); listt(y, 1, 10); listt(y, nn-10, nn); */
	if(shifting) shift(nn);
	
	fprintf(stderr, "computing FFT...\n");
	fft(&nn, x, y);

	if(output)
		{if(output<32) output = nn/output;
		if(output<nn) nn=output;
		}
	if(abscissa_arguments)
		fprintf(ofile, "; output frequency step is %15.8g", x[2]-x[1]);
	i=1;
	while(1)
		{if(!automatic_abscissas)  fprintf(ofile, "%15.6g ", x[i]);
		yr=y[2*i-1]; yi=y[2*i];
		if(magnitudes)  fprintf(ofile, "%15.6g", sqrt(yr*yr + yi*yi));
		else            fprintf(ofile, "%15.6g %15.6g", yr, yi);
		if(++i>nn) break;
		fprintf(ofile, "\n");
		}
	fprintf(ofile, " \"\"\n");
	if(ofile!=stdout)
		{fclose(ofile);
		}
}

double energy()
{	double v, e;
	int i;
	i=1;
	e=0.;
	while(i<=m)
		{v=y[i++]; e=e+v*v;
		}
	return e*(x[2]-x[1]);
}

pad()		/* pad the data with zeros */
{	int i, k, num, max;
	char c;
	if(keeping)
		m=keeping;
	else 
		m=n;
/*
		while(1)
			{fprintf(stderr,
			"Last point to include (2...%d, default %d) ? ", n, n);
			gets(buf);
			if(buf[0]) sscanf(buf, "%d", &m);
			else m=n;
			if((m>=2) && (m<=n)) break;
			}
*/
	if(pad_factor)
		{if(pad_factor>=32) num=(pad_factor/10)*7;
		else if(pad_factor>0) num=((m*pad_factor)/10)*7;
		else num=m*4;	/*  default is a padding factor between 4 and 8 */
		}
	else num=m*4;	/*  default is a padding factor between 4 and 8 */
	for (k=1; k<num; k <<= 1) {}	/* find next power of 2 */
	while (k>ENTRIES) {k >>= 1;}	/* back off if too many for buffer */
	fprintf(stderr, "transforming %d points, for actual padding factor of %3.1f\n", 
		k, k/(double)m);
	i=m;
	while(i<k)	/* will have k values after this */
		{y[++i]=0.;}
	return k;
}

shift(n) int n;		/* shift data so time "origin" is at the beginning of the array */
{	int k;
	k=(int)((origin-x[1])/(x[2]-x[1]) + 1.5);
	if(k<1 || k>n) 
		{fprintf(stderr, "specified origin isn't within abscissa range");
		exit(1);
		}
			/* need to swap y[1]...y[k-1] with y[k]...y[n] */
	reverse(1, k-1);
	reverse(k, n);
	reverse(1, n);
}

reverse(i, j) int i,j;  /* reverse y[i]...y[j]  */
{	double t;
	while(i<j) {t=y[i]; y[i]=y[j]; y[j]=t; i++; j--;}
}

wind16()
{	char c;
	int j;
	double p, p2, dp, scale, epsilon;

/*	if(!super)
		{while(1)
			{puts("Perform one-sided supergaussian windowing, degree 16? ");
			c=toupper(getchar()); putchar('\n');
			if(c=='N')return 0;
			if(c=='Y')break;
			}
		}
	else
*/
	     if(super!=16) return 0;		
/*
	if(wt==0.)
		{fprintf(stderr, "Weight at window edge (0 < epsilon < 1, default=%f) ? ",
			DEFAULT_EDGE_WEIGHT);
		gets(buffer);
		if(buffer[0]) epsilon=atof(buffer); else epsilon=DEFAULT_EDGE_WEIGHT;
		}
	else 
*/		epsilon=wt;
	p=exp(log(-log(epsilon))/16.);	/* (-ln(epsilon))**(1/16) */
/*	fprintf(stderr, "ratio C/A = %7.4f \n", p);  */
	dp=p/(m-1);
	j=m;
	while(j)
		{p2=p*p; p2=p2*p2; p2=p2*p2;	/* p**8 */
		y[j] *= exp(-p2*p2);	/* p**16 */
		p -= dp; --j;
		}
	return 1;
}

wind10()
{	char c;
	int j;
	double p, p2, p4, dp, scale, epsilon;

/*	if(!super)
		{while(1)
			{puts("Perform one-sided supergaussian windowing, degree 10? ");
			c=toupper(getchar()); putchar('\n');
			if(c=='N')return 0;
			if(c=='Y')break;
			}
		}
	else
*/
	     if(super!=10) return 0;		
/*
	if(wt==0.)
		{fprintf(stderr,"Weight at window edge (0 < epsilon < 1, default=%f) ? ",
			DEFAULT_EDGE_WEIGHT);
		gets(buffer);
		if(buffer[0]) epsilon=atof(buffer); else epsilon=DEFAULT_EDGE_WEIGHT;
		}
	else 
*/
		epsilon=wt;
	p=exp(log(-log(epsilon))/10.);	/* (-ln(epsilon))**(1/10) */
/*	fprintf(stderr, "ratio C/A = %7.4f \n", p);  */
	dp=p/(m-1);
	j=m;
	while(j)
		{p2=p*p; p4=p2*p2;	/* p**2, p**4 */
		y[j] *= exp(-p2*p4*p4);	/* p**10 */
		p -= dp; --j;
		}
	return 1;
}

wind6()
{	char c;
	int j;
	double p, p2, dp, scale, epsilon;

/*	if(!super)
		{while(1)
			{puts("Perform one-sided supergaussian windowing, degree 6? ");
			c=toupper(getchar()); putchar('\n');
			if(c=='N')return 0;
			if(c=='Y')break;
			}
		}
	else
*/
	     if(super!=6) return 0;		
/*
	if(wt==0.)
		{fprintf(stderr,"Weight at window edge (0 < epsilon < 1, default=%f) ? ",
			DEFAULT_EDGE_WEIGHT);
		gets(buffer);
		if(buffer[0]) epsilon=atof(buffer); else epsilon=DEFAULT_EDGE_WEIGHT;
		}