fft.c

傅立叶变化

	else 
*/
		epsilon=wt;
	p=exp(log(-log(epsilon))/6.);	/* (-ln(epsilon))**(1/6) */
/*	fprintf(stderr, "ratio C/A = %7.4f \n", p);  */
	dp=p/(m-1);
	j=m;
	while(j)
		{p2=p*p*p;  /* p**3 */
		y[j] *= exp(-p2*p2);	/* p**6 */
		p -= dp; --j;
		}
	return 1;
}

listt(y, i, j) FLOAT y[]; int i, j;	/* display y[i] through y[j]  */
{	while (i<=j) {printf("%4d %15.9f \n", i, y[i]); ++i;}
}


read_data(argc, argv) int argc; char **argv;
{	int i, j, length;
	double xx, yy, d, *pd, sum;
	char *s, *t, *stop;
	char *malloc();
	char *strchr();

	x=(FLOAT *)malloc(ENTRIES*sizeof(FLOAT));
	y=(FLOAT *)malloc(ENTRIES*sizeof(FLOAT));
	if(x==0 || y==0) {fprintf(stderr, "can\'t allocate buffer"); exit(1);}
	argc--; argv++;
	if(strchr(argv[0], '?')) help();
/*            open input file         */
	if(argc && *argv[0]!='-')
		{ifile=fopen(argv[0], "r");
		if(ifile==0) {fprintf(stderr, "file %s not found\n", argv[0]); exit(1);}
		strcpy(iname, argv[0]);
		argc--; argv++;
		}
/*            open output file         */
	if(argc && *argv[0]!='-')
		{strcpy(oname, argv[0]);
		argc--; argv++;
		unlink(oname);
		if((ofile=fopen(oname, "w"))==0)
			{fprintf(stderr, "can\'t open output file %s", oname);
			exit(1);
			}		
		}
	p_data[0]=-1;

	fprint_cmd(ofile, "; fft %s\n");

	i=1;		/* note x[0] and y[0] aren't defined */
	while(i<ENTRIES)
		{if(fgets(buf, BUFSIZE, ifile)==0)
			{if(standard_input) {fclose(ifile); ifile=fopen("/dev/con", "r");}
			break;
			}
		t=buf;
		while(*t && isspace(*t)) t++;
		if(*t == 0) continue;		/* skip blank lines */
		buf[strlen(buf)-1]=0;		/* zero out the line feed */
		if(buf[0]==';')				/* skip comment */
			{fprintf(ofile, "%s\n", buf); continue;
			}
		if(t = strchr(buf, ';')) *t=0;	/* zap same-line comment */
		if(automatic_abscissas)
			{xx=abscissa;
			abscissa+=abscissa_step;
			sscanf(buf, "%lf", &yy); 
			}
		else
			{
			sscanf(buf, "%lf %lf", &xx, &yy);
			}
		range_check(xx); range_check(yy);
		x[i]=xx; y[i]=yy;		/* convert doubles to floats here */
		s=buf;                      /* start looking for label */
		while(*s==' ')s++;			/* skip first number */
		while(*s && (*s!=' '))s++;
		if(!automatic_abscissas)	/* skip second number */
			{while(*s==' ')s++;
			while(*s && (*s!=' '))s++;
			}
		while(*s==' ')s++;
		if((length=strlen(s))&&(labels<MAXLABELS))
			{if(*s=='\"')           /* label in quotes */
				{t=s+1;
				while(*t && (*t!='\"')) t++;
				t++;
				}
			else                    /* label without quotes */
				{t=s;
				while(*t && (*t!=' '))t++;
				}
			*t=0;
			length=t-s;
			p_data[labels]=i;
			p_text[labels]=(char *)malloc(length+1);
			if(p_text[labels]) strcpy(p_text[labels++], s);
			}
		i++;
		}
	n=i-1;
/*	if(ifile==stdin) {close(ifile); ifile=fopen();} */
	if(breaking && (!labels || p_data[labels-1]!=n-1))
		{p_data[labels]=i-1;
		if(p_text[labels]=(char *)malloc(1)) *p_text[labels++]=0;
		}
}

/* check whether number is too big for a float */
range_check(x) double x;
{	if(fabs(x)>BIGFLOAT)
		{fprintf(stderr, "input number too large: %f\n", x);
		exit(1);
		}
}


int streq(a,b) char *a,*b;
{	while(*a)
		{if(*a!=*b)return 0;
		a++; b++;
		}
	return (*b==0);
}

/* get_parameter - process one command line option
		(return # parameters used) */
get_parameter(argc, argv) int argc; char **argv;
{	int i;
	if(streq(*argv, "-a"))
		{i=get_double(argc, argv, 1, &abscissa_step, &abscissa, &abscissa);
		abscissa_arguments=i-1;
		automatic_abscissas=1;
		return i;
		}
	else if(streq(*argv, "-b")) {breaking=1; return 1;}
	else if(streq(*argv, "-s"))
		{super=10; wt=DEFAULT_EDGE_WEIGHT;
		if((argc>1) && numeric(argv[1])) 
			{super=atoi(argv[1]);
			if((super!=16) && (super!=10) && (super!=6))
				{fprintf(stderr, "windowing only of degrees 6, 10, & 16 available");
				exit(1);
				}
			if((argc>2) && numeric(argv[2]))
				{wt=atof(argv[2]);
				if((wt<=0.)||(wt>=1.))
					{fprintf(stderr, "wt=%f, need  0 < wt < 1", wt); exit(1);
					}
				return 3;
				}
			return 2;
			}
		return 1;
		}
	else if(streq(*argv, "-n"))
		{if((argc>1) && numeric(argv[1])) keeping=atoi(argv[1]);
		if(keeping<2 || keeping>ENTRIES)
			{fprintf(stderr, "switch -n option out of range 2...%d\n",ENTRIES);
			exit(1);
			}
		return 2;
		}
	else if(streq(*argv, "-o"))
		{if((argc>1) && numeric(argv[1])) output=atof(argv[1]);
		if(output<1 || output>ENTRIES)
			{fprintf(stderr, "switch -o option out of range 1...%d\n",ENTRIES);
			exit(1);
			}
		return 2;
		}
	else if(streq(*argv, "-p"))
		{if((argc>1) && numeric(argv[1])) pad_factor=atoi(argv[1]);
		return 2;
		}
	else if(streq(*argv, "-f"))
		{i=get_double(argc, argv, 2, &fmin, &fmax, &fmax);
		f_arguments=i-1;
		return i;
		}
	else if(streq(*argv, "-t"))
		{i=get_double(argc, argv, 2, &xmin, &xmax, &xmax);
		x_arguments=i-1;
		return i;
		}
	else if(streq(*argv, "-m"))
		{magnitudes=1;
		return 1;
		}
	else if(streq(*argv, "-z"))
		{origin=0.;
		i=get_double(argc, argv, 1, &origin, &fmax, &fmax);
		shifting=1;
		return i;
		}
	else gripe(argv);
}

char *message[]={
"fft   version ", VERSION,
" - calculate fast Fourier transform of time domain curve\n",
"usage: fft  [infile  [outfile]]  [options]\n",
"options are:\n",
"  -a  [step]       automatic abscissas \n",
/*
"  -b               break input after each label,  \n",
"                   find separate transforms\n",
"  -f  min [max]    minimum and maximum frequencies\n",
"  -t  min [max]    minimum and maximum times\n",
*/
"  -m               print only frequencies & magnitudes\n",
"  -p  num          pad data by factor of num\n",
"  -n  num          keep only first  num  input data points\n",
"  -o  num          print out num values\n",
"  -s  [deg [wt]]   perform supergaussian windowing of\n",
"             degree deg (6, 10, or 16, default 10) with\n",
"             weight wt on last point (default .9)\n",
"  -z  origin       subtract  origin  from each abscissa value\n",
0
};

/*	name...
		fft

	purpose...
		perform forward or inverse FFT

	history...
		May 84  translated from FORTRAN
		12 Jun 84  diagnostic printouts shortened
*/

	int ip, n2, i, j, k, nu, nn;
	int k1n2, l, nu1;
	double arg, p, q, s, c;
	double ninv, time, dt, freq, df;
	double xplus, xminus, yplus, yminus;
	double uplus, uminus, vplus, vminus;

	/* format time domain function for calculation
		of its fast Fourier transform */
/*
on input:
n = a power of 2
x[1...n] has  t(0)  t(1)  t(2)  ...  t(n-1)
y[1...n] has  q(0)  q(1)  q(2)  ...  q(n-1)

on output:
n = (n/2+1), one more than a power of 2
x[1...n]  has  f(0)  f(1)  f(2)  ...  f(n-1)
y[1...2n] has  Qr(0) Qi(0)   Qr(1) Qi(1)   Qr(2) Qi(2)  ...  Qr(n-1)  Qi(n-1)
*/
fft(nnn, x, y) int *nnn; FLOAT x[], y[];
{	n= *nnn;
	if(n<=0){puts("fft: illegal # points"); return;}
	dt=x[2]-x[1];
	df=1./(n*dt);
	n2=4; nu=0;
	while(++nu<=20)
		{if(n2==n)break;
		n2=n2+n2;
		}
	if(nu>20){puts("fft: n not a power of 2"); return;}
/*	puts("fft: n is ok, nu= %d\n", nu);	*/
	n=n/2; i=0;
	while(++i<=n) {x[i]=y[2*i]; y[i]=y[2*i-1];}
/*
	printf("fft: data reformatted\n");
	i=0; 
	while(++i<=10) {printf("%6d %15.6f %15.6f\n", i, y[i], x[i]);}
*/
	fft2(y, x, n, nu);
/*
	puts("fft: transform calculated\n"); 
	i=0; 
	while(++i<=10) {printf("%6d %15.6f %15.6f\n", i, y[i], x[i]);}
*/
	ip=n+1; y[ip]=y[1]; x[ip]=x[1];
	ninv=3.14159265358979/n;
	n=n/2; n2=n+1; i=0;
	while(++i<=n2)
		{arg=ninv*(i-1); s=sin(arg); c=cos(arg);
		uplus=y[i]+y[ip]; uminus=y[i]-y[ip];
		vplus=x[i]+x[ip]; vminus=x[i]-x[ip];
		p= c*vplus-s*uminus; q= -s*vplus-c*uminus;
		y[i]=.5*(uplus+p); y[ip]=.5*(uplus-p);
		x[i]=.5*(q+vminus); x[ip]=.5*(q-vminus);
		--ip;
		}
	n=n+n+1; ip=n;
	while(ip) {y[ip+ip-1]=dt*y[ip]; --ip;}
	freq=0.; i=0;
	while(++i<=n) {y[i+i]=x[i]*dt; x[i]=freq; freq=freq+df;}
	*nnn=n;
}

fft2(xreal, ximag, n, nu)
FLOAT xreal[], ximag[];
int n, nu;
{	double treal, timag;
	FLOAT *xr1, *xr2, *xi1, *xi2;

	n2=n/2; nu1=nu-1; k=0; l=0;
	while(++l<=nu)
		{while(k<n)
			{p=ibitr(k>>nu1, nu);
			arg=3.14159265358979*2*p/n;
			c=cos(arg); s=sin(arg); i=0;
			while(++i<=n2)
			    {++k; k1n2=k+n2;
			    /* initialize four pointers */
			    xr1=xreal+k; xr2=xreal+k1n2;
			    xi1=ximag+k; xi2=ximag+k1n2;
			    treal= *xr2*c+ *xi2*s;
			    timag= *xi2*c- *xr2*s;
			    *xr2= *xr1-treal;
			    *xi2= *xi1-timag;
			    *xr1= *xr1+treal;
			    *xi1= *xi1+timag;
			    }
			k=k+n2;
			}
		k=0; --nu1; n2=n2/2;
		}
	k=0;
	while(++k<=n)
		{i=ibitr(k-1, nu)+1;
		if(i>k)
			{treal=xreal[k];   timag=ximag[k];
			xreal[k]=xreal[i]; ximag[k]=ximag[i];
			xreal[i]=treal;    ximag[i]=timag;
			}
		}
}

/* reverse the last nu bits of j */
ibitr(j, nu) int j, nu;
{	int ib;
	ib=0;
	while(nu--){ib=(ib<<1)+(j&1); j=j>>1;}
	return ib;
}