mh@killer.UUCP (Mike Hobgood) (08/07/87)
ptsfa!dmt says: > I had a problem unshar'ing the fft program you posted to netnews. > complex.c has no } at the end, and I believe that there is another > problem in polar(). > Could you check your source and repost it? > Thanks (and thanks for the original posting). Egad! I uploaded a "messed with" version. Here's the original: (and I don't know if it compiles, I didn't try). Have fun. ----------------------------- fft.c ------------------------------------ /* Performs a Cooley-Tukey type fast fourier transform. original version Pascal written in "Simple Calculations with Complex Numbers by David Clark DDJ 10/84 translated to C by R. Hellman 02/21/86 released to public domain --------------------------------------------------------------------------- Differences from original version: as in fortran versions that I have worked with instead of specifying total array size in decimal instead specify number of array members as power of 2. Forward or reverse fft is specified by the sign of this argument. ****IMPORTANT**** for convenience sake members of an array are taken to begin at [1] rather than [0] thus for an fft on a 64 member array of pointers to complex numbers there must be 65 members. ---------------------------------------------------------------------------- f is a one dimensional array of pointers to complex numbers with length (2**abs(ln)). The sign of ln determines in which direction the transform will be performed. If ln is positive, then isinverse is set to -1 and a forward transform is carried out. If ln is negative then isinverse isinverse is set to 1 and a reverse transform is calculated. for numpoints=(2**abs(ln)),transform[j]=sum(f[i]*w^((i-1)*(j-1))), where i, j are 1 to numpoints and w=exp(isinverse*2*pi*sqrt(-1)/numpoints). The program also computes the inverse transform, for which the defining expression is: inverse DFT=(1/numpoints)*sum(f[i]*w^((i-1)*(j-1))). Run time is proportional to numpoints*log2(numpoints) rather than to numpoints**2 for the classical discrete Fourier transform. */ #include "stdio.h" #include "math.h" typedef struct { double re; double im; } COMPLEX; fft(f,ln) COMPLEX *f[]; int ln; { int i,isinverse,numpoints; if (ln>0) { printf("Calculating forward Fourier transform\n"); isinverse = -1; numpoints = power(2,ln); } else { printf("Calculating reverse Fourier transform\n"); isinverse = 1; numpoints = power(2,-ln); } scramble(numpoints,f); butterflies(numpoints,isinverse,f); } static scramble(numpoints,f) int numpoints; COMPLEX *f[]; /* put the data in bit reversed order */ { int i,j,m; COMPLEX *temp; j = 1; for (i=1;i<=numpoints;i++) { if (i<j) { temp = f[j]; f[j] = f[i]; f[i] = temp; } m = numpoints>>1; if (m<j) while (m<j) { j -= m; m = (m + 1)>>1; } j += m; } } static butterflies(numpoints,isinverse,f) int numpoints,isinverse; COMPLEX *f[]; /* calculate the butterflies for the bit reversed data of an fft - normalize if a reverse fft is performed */ { int mmax,step,index,m,i,j; double theta; COMPLEX w,cn,temp; mmax = 1; while (numpoints>mmax) { step = mmax<<1; for (m=1;m<=mmax;m++) { theta = PI*((double)(isinverse)*(double)(m-1))/(double)(mmax); w.re = cos(theta); w.im = sin(theta); for (i=1;i<=((numpoints-m)/step)+1;i++) { index = m + ((i-1)*step); j = index + mmax; cmult(&temp,&w,f[j]); csub(f[j],f[index],&temp); cadd(f[index],f[index],&temp); } } mmax = step; } if (isinverse==1) { cn.re = numpoints; cn.im = 0.0; for (i=1;i<=numpoints;i++) cdiv(f[i],f[i],&cn); } } static power(x,y) int x,y; { int i,j; j = 1; for (i=0;i<y;i++) j *= x; return(j); } ------------------------------ complex.c --------------------------------- /* Simple Complex number functions V1.0 Author: R. Hellman Date: 02/21/86 Released for public domain use only */ #include "stdio.h" #include "math.h" /* defines done in math.h (Lattice C) #define PI 3.14159265358979323846 #define PID2 1.57079632679489661923 /* PI divided by 2 */ #define TINY 2.2e-308 /* tiny value */ #define LOGHUGE 709.778 /* natural log of huge value */ math.h also identifies appropriate functions as double */ typedef struct { double re; double im; } COMPLEX; initcmplx(f,arraysz) COMPLEX *f[]; int arraysz; /* initializes storage for each address of an array of complex pointers *f[] with 'arraysz' number of elements */ { int i; char *getmem(); for (i=0;i<arraysz;i++) #ifdef LATTICEC if ((f[i]=(COMPLEX *)getmem(sizeof(COMPLEX)))==NULL) { write(stderr,"Error: *COMPLEX number allocation\n",35); exit(1); } #else if ((f[i]=(COMPLEX *)malloc(sizeof(COMPLEX)))==NULL) { write(stderr,"Error: *COMPLEX number allocation\n",35); exit(1); } #endif } /* COMPLEX number library translated from Pascal source listing in Dr. DOBBS October, 1984 "Simple Calculations with Complex Numbers" by David Clark */ /* to insure that the same structure could be used as both an 'argument' and 'result' immediate tranfer of structure members is done to a local variable - in many of the functions this may not be necessary, but it was done in all of the following functions for consistency */ cadd(result,arg1,arg2) COMPLEX *result,*arg1,*arg2; /* adds arg1 and arg2 and returns sum in result */ { COMPLEX targ1,targ2; targ1.re = arg1->re; targ1.im = arg1->im; targ2.re = arg2->re; targ2.im = arg2->im; result->re = targ1.re + targ2.re; result->im = targ1.im + targ2.im; } csub(result,arg1,arg2) COMPLEX *result,*arg1,*arg2; /* subtracts arg1 and arg2 and returns sum in result */ { COMPLEX targ1,targ2; targ1.re = arg1->re; targ1.im = arg1->im; targ2.re = arg2->re; targ2.im = arg2->im; result->re = targ1.re - targ2.re; result->im = targ1.im - targ2.im; } cmult(result,arg1,arg2) COMPLEX *result,*arg1,*arg2; /* multiplies arg1 and arg2 and returns product in result */ { COMPLEX targ1,targ2; targ1.re = arg1->re; targ1.im = arg1->im; targ2.re = arg2->re; targ2.im = arg2->im; result->re = (targ1.re * targ2.re) - (targ1.im * targ2.im); result->im = (targ1.im * targ2.re) + (targ1.re * targ2.im); } cdiv(result,arg1,arg2) COMPLEX *result,*arg1,*arg2; /* divides arg1 and arg2 and returns quotient in result */ { double denom; COMPLEX targ1,targ2; targ1.re = arg1->re; targ1.im = arg1->im; targ2.re = arg2->re; targ2.im = arg2->im; denom = (targ2.re*targ2.re) + (targ2.im*targ2.im); result->re = (targ1.re*targ2.re + targ1.im*targ2.im)/denom; result->im = (targ1.im*targ2.re - targ1.re*targ2.im)/denom; } polar(arg,modulus,amplitude) COMPLEX *arg; double *modulus,*amplitude; /* converts a complex number arg in rectangular form to polar form */ { COMPLEX targ; targ.re = arg->re; targ.im = arg->im; if (fabs(targ.re) < TINY) targ.re = 0.0; if (fabs(targ.im) < TINY) targ.im = 0.0; *modulus = sqrt((targ.re*targ.re) + (targ.im*targ.im)); if (targ.im == 0.0) *amplitude = 0.0; else { if (targ.re == 0.0) { if (targ.im > 0.0) *amplitude = PID2; else *amplitude = -PID2; } else { if ((log(fabs(targ.im)) - log(fabs(targ.re))) > LOGHUGE) { if (targ.re > 0.0) { if (targ.im > 0.0) *amplitude = PID2; else *amplitude = -PID2; } else { if (targ.im > 0.0) *amplitude = -PID2; else *amplitude = PID2; } } else *amplitude = atan(targ.im/targ.re); } } } ctopower(result,arg,power) COMPLEX *result,*arg; int power; /* raises arg to the positive integral power and returns answer in result */ { int i; double modulus,amplitude,newmod,powamp; COMPLEX targ; targ.re = arg->re; targ.im = arg->im; if (power==0) { result->re = 1.0; result->im = 0.0; } else { polar(&targ,&modulus,&litude); newmod = 1.0; if (power > 0) for (i=0;i<power;i++) newmod = newmod*modulus; else for (i=0;i<abs(power);i++) newmod = newmod/modulus; powamp = power*amplitude; result->re = newmod*cos(powamp); result->im = newmod*sin(powamp); } } cexp(result,arg) COMPLEX *result,*arg; /* raises e to the arg and returns the answer in result */ { double expo; COMPLEX targ; targ.re = arg->re; targ.im = arg->im; expo = exp(targ.re); result->re = expo*cos(targ.im); result->im = expo*sin(targ.im); } cln(result,arg) COMPLEX *result,*arg; /* takes the natural logarithm of arg and returns the answer in result */ { double modulus,amplitude; COMPLEX targ; targ.re = arg->re; targ.im = arg->im; polar(&targ,&modulus,&litude); result->re = log(modulus); result->im = amplitude; } ctoc(result,arg1,arg2) COMPLEX *result,*arg1,*arg2; /* raise complex number arg1 to complex power arg2 and return answer in result */ { COMPLEX logpart,expo; COMPLEX targ1,targ2; targ1.re = arg1->re; targ1.im = arg1->im; targ2.re = arg2->re; targ2.im = arg2->im; cln(&logpart,&targ1); cmult(&expo,&targ2,&logpart); cexp(result,&expo); } csin(result,arg) COMPLEX *result,*arg; /* takes the sine of arg and returns it in result */ { COMPLEX targ,exp1,exp2,part1,part2,sum,divisor; targ.re = arg->re; targ.im = arg->im; exp1.re = -targ.im; exp1.im = targ.re; cexp(&part1,&exp1); exp2.re = targ.im; exp2.im = -targ.re; cexp(&part2,&exp2); csub(&sum,&part1,&part2); divisor.re=0.0; divisor.im=2.0; cdiv(result,&sum,&divisor); } ccos(result,arg) COMPLEX *result,*arg; /* takes the cosine of arg and returns it in result */ { COMPLEX targ,exp1,exp2,part1,part2,sum,divisor; targ.re = arg->re; targ.im = arg->im; exp1.re = -targ.im; exp1.im = targ.re; cexp(&part1,&exp1); exp2.re = targ.im; exp2.im = -targ.re; cexp(&part2,&exp2); cadd(&sum,&part1,&part2); divisor.re = 2.0; divisor.im = 0.0; cdiv(result,&sum,&divisor); } ----------------------------------- EOF -------------------------------------