home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume1 / 8708 / 5 < prev    next >
Encoding:
Internet Message Format  |  1990-07-13  |  11.3 KB
  1. From: mh@killer.UUCP (Mike Hobgood)
  2. Newsgroups: comp.sources.misc
  3. Subject: Original fft listing
  4. Keywords: Oops...
  5. Message-ID: <4080@ncoast.UUCP>
  6. Date: 7 Aug 87 01:07:27 GMT
  7. Sender: allbery@ncoast.UUCP
  8. Organization: The Unix(tm) Connection, Dallas, Texas
  9. Lines: 456
  10. Approved: allbery@ncoast.UUCP
  11. X-Archive: comp.sources.misc/8708/5
  12.  
  13.  
  14. ptsfa!dmt  says:
  15.  
  16. > I had a problem unshar'ing the fft program you posted to netnews.
  17. > complex.c has no } at the end, and I believe that there is another
  18. > problem in polar().
  19.  
  20. > Could you check your source and repost it?
  21. > Thanks (and thanks for the original posting).
  22.  
  23. Egad!  I uploaded a "messed with" version.  Here's the original:
  24. (and I don't know if it compiles, I didn't try).  Have fun.
  25.  
  26. ----------------------------- fft.c ------------------------------------
  27.  
  28. /*         Performs a Cooley-Tukey type fast fourier transform.
  29.  
  30.     original version Pascal written in "Simple Calculations with 
  31.     Complex Numbers by David Clark DDJ 10/84
  32.     
  33.     translated to C by R. Hellman
  34.     02/21/86
  35.     released to public domain
  36.  
  37. ---------------------------------------------------------------------------
  38.     Differences from original version:
  39.  
  40.     as in fortran versions that I have worked with instead of specifying total
  41.     array size in decimal instead specify number of array members as power of
  42.     2. Forward or reverse fft is specified by the sign of this argument.
  43.  
  44.     ****IMPORTANT****
  45.     for convenience sake members of an array are taken to begin at [1] rather
  46.     than [0] thus for an fft on a 64 member array of pointers to complex
  47.     numbers there must be  65 members.
  48. ----------------------------------------------------------------------------
  49.  
  50. f is a one dimensional array of pointers to complex numbers with length
  51. (2**abs(ln)).  The sign of ln determines in which direction the transform
  52. will be performed.  If ln is positive, then isinverse is set to -1 and
  53. a forward transform is carried out.  If ln is negative then isinverse 
  54. isinverse is set to 1 and a reverse transform is calculated.
  55.  
  56. for numpoints=(2**abs(ln)),transform[j]=sum(f[i]*w^((i-1)*(j-1))), where i,
  57. j are 1 to numpoints and w=exp(isinverse*2*pi*sqrt(-1)/numpoints).  The
  58. program also computes the inverse transform, for which the defining
  59. expression is: inverse DFT=(1/numpoints)*sum(f[i]*w^((i-1)*(j-1))).
  60.  
  61.      Run time is proportional to numpoints*log2(numpoints) rather than
  62. to numpoints**2 for the classical discrete Fourier transform. */
  63.  
  64. #include "stdio.h"
  65. #include "math.h"
  66.  
  67. typedef struct
  68.     {
  69.     double re;
  70.     double im;
  71.     }
  72. COMPLEX;
  73.  
  74. fft(f,ln)
  75. COMPLEX *f[];
  76. int ln;
  77. {
  78.     int i,isinverse,numpoints;
  79.     
  80.         if (ln>0) {
  81.             printf("Calculating forward Fourier transform\n");
  82.             isinverse = -1;
  83.             numpoints = power(2,ln);
  84.             }
  85.         else {
  86.             printf("Calculating reverse Fourier transform\n");
  87.             isinverse = 1;
  88.             numpoints = power(2,-ln);
  89.             }
  90.     scramble(numpoints,f);
  91.     butterflies(numpoints,isinverse,f);
  92. }
  93.  
  94. static scramble(numpoints,f)
  95. int numpoints;
  96. COMPLEX *f[];
  97.  
  98. /* put the data in bit reversed order */
  99.  
  100. {
  101.     int i,j,m;
  102.     COMPLEX *temp;
  103.  
  104.     j = 1;
  105.     for (i=1;i<=numpoints;i++) {
  106.         if (i<j) {
  107.             temp = f[j];
  108.             f[j] = f[i];
  109.             f[i] = temp;
  110.             }
  111.         m = numpoints>>1;
  112.         if (m<j)
  113.             while (m<j) {
  114.                 j -= m;
  115.                 m = (m + 1)>>1;
  116.                 }
  117.         j += m;
  118.         }
  119. }
  120.  
  121. static butterflies(numpoints,isinverse,f)
  122. int numpoints,isinverse;
  123. COMPLEX *f[];
  124.  
  125. /* calculate the butterflies for the bit reversed data of an fft -
  126.   normalize if a reverse fft is performed */
  127.  
  128. {
  129.    int mmax,step,index,m,i,j;
  130.    double theta;
  131.    COMPLEX w,cn,temp;
  132.  
  133.     mmax = 1;
  134.     while (numpoints>mmax) {
  135.         step = mmax<<1;
  136.         for (m=1;m<=mmax;m++) {
  137.             theta = PI*((double)(isinverse)*(double)(m-1))/(double)(mmax);
  138.             w.re = cos(theta);
  139.             w.im = sin(theta);
  140.             for (i=1;i<=((numpoints-m)/step)+1;i++) {
  141.                 index = m + ((i-1)*step);
  142.                 j = index + mmax;
  143.                 cmult(&temp,&w,f[j]);
  144.                 csub(f[j],f[index],&temp);
  145.                 cadd(f[index],f[index],&temp);
  146.                 }
  147.             }
  148.         mmax = step;
  149.         }
  150.     if (isinverse==1) {
  151.         cn.re = numpoints;
  152.         cn.im = 0.0;
  153.         for (i=1;i<=numpoints;i++)
  154.             cdiv(f[i],f[i],&cn);
  155.         }
  156. }
  157.  
  158. static power(x,y)
  159. int x,y;
  160. {
  161.     int i,j;
  162.     
  163.     j = 1;
  164.     for (i=0;i<y;i++)
  165.         j *= x;
  166.     return(j);
  167. }
  168.  
  169. ------------------------------ complex.c ---------------------------------
  170.  
  171. /*  Simple Complex number functions V1.0
  172.     Author: R. Hellman
  173.     Date: 02/21/86
  174.     Released for public domain use only */
  175.  
  176. #include "stdio.h"
  177. #include "math.h"
  178.  
  179. /* defines done in math.h (Lattice C)
  180.  
  181. #define PI   3.14159265358979323846
  182. #define PID2 1.57079632679489661923    /* PI divided by 2 */
  183. #define TINY 2.2e-308            /* tiny value */
  184. #define LOGHUGE 709.778            /* natural log of huge value */
  185.  
  186. math.h also identifies appropriate functions as double
  187.  
  188. */
  189.  
  190. typedef struct
  191.     {
  192.     double re;
  193.     double im;
  194.     }
  195.     COMPLEX;
  196.  
  197. initcmplx(f,arraysz)
  198. COMPLEX *f[];
  199. int arraysz;
  200.  
  201. /* initializes storage for each address of an array of complex pointers *f[]
  202.    with 'arraysz' number of elements */
  203.  
  204. {
  205.     int i;
  206.     char *getmem();
  207.  
  208.     for (i=0;i<arraysz;i++)
  209. #ifdef LATTICEC
  210.         if ((f[i]=(COMPLEX *)getmem(sizeof(COMPLEX)))==NULL) {
  211.             write(stderr,"Error: *COMPLEX number allocation\n",35);
  212.             exit(1);
  213.             }
  214. #else
  215.         if ((f[i]=(COMPLEX *)malloc(sizeof(COMPLEX)))==NULL) {
  216.             write(stderr,"Error: *COMPLEX number allocation\n",35);
  217.             exit(1);
  218.             }
  219. #endif
  220. }
  221.  
  222. /* COMPLEX number library translated from Pascal source listing in Dr. DOBBS
  223.    October, 1984 "Simple Calculations with Complex Numbers" by David Clark */
  224.  
  225. /* to insure that the same structure could be used as both an 'argument' and 
  226.    'result' immediate tranfer of structure members is done to a local variable
  227.    - in many of the functions this may not be necessary, but it was done in
  228.    all of the following functions for consistency */
  229.  
  230. cadd(result,arg1,arg2)
  231. COMPLEX *result,*arg1,*arg2;
  232.  
  233. /* adds arg1 and arg2 and returns sum in result */
  234.  
  235. {
  236.     COMPLEX targ1,targ2;
  237.  
  238.     targ1.re = arg1->re;
  239.     targ1.im = arg1->im;
  240.     targ2.re = arg2->re;
  241.     targ2.im = arg2->im;
  242.     result->re = targ1.re + targ2.re;
  243.     result->im = targ1.im + targ2.im;
  244. }
  245.  
  246. csub(result,arg1,arg2)
  247. COMPLEX *result,*arg1,*arg2;
  248.  
  249. /* subtracts arg1 and arg2 and returns sum in result */
  250.  
  251. {
  252.     COMPLEX targ1,targ2;
  253.  
  254.     targ1.re = arg1->re;
  255.     targ1.im = arg1->im;
  256.     targ2.re = arg2->re;
  257.     targ2.im = arg2->im;
  258.     result->re = targ1.re - targ2.re;
  259.     result->im = targ1.im - targ2.im;
  260. }
  261.  
  262. cmult(result,arg1,arg2)
  263. COMPLEX *result,*arg1,*arg2;
  264.  
  265. /* multiplies arg1 and arg2 and returns product in result */
  266.  
  267. {
  268.     COMPLEX targ1,targ2;
  269.  
  270.     targ1.re = arg1->re;
  271.     targ1.im = arg1->im;
  272.     targ2.re = arg2->re;
  273.     targ2.im = arg2->im;
  274.     result->re = (targ1.re * targ2.re) - (targ1.im * targ2.im);
  275.     result->im = (targ1.im * targ2.re) + (targ1.re * targ2.im);
  276. }
  277.  
  278. cdiv(result,arg1,arg2)
  279. COMPLEX *result,*arg1,*arg2;
  280.  
  281. /* divides arg1 and arg2 and returns quotient in result */
  282.  
  283. {
  284.     double denom;
  285.     COMPLEX targ1,targ2;
  286.  
  287.     targ1.re = arg1->re;
  288.     targ1.im = arg1->im;
  289.     targ2.re = arg2->re;
  290.     targ2.im = arg2->im;
  291.     denom = (targ2.re*targ2.re) + (targ2.im*targ2.im);
  292.     result->re = (targ1.re*targ2.re + targ1.im*targ2.im)/denom;
  293.     result->im = (targ1.im*targ2.re - targ1.re*targ2.im)/denom;
  294. }
  295.  
  296. polar(arg,modulus,amplitude)
  297. COMPLEX *arg;
  298. double *modulus,*amplitude;
  299.  
  300. /* converts a complex number arg in rectangular form to polar form */
  301.  
  302. {
  303.     COMPLEX targ;
  304.  
  305.     targ.re = arg->re;
  306.     targ.im = arg->im;
  307.     if (fabs(targ.re) < TINY)
  308.         targ.re = 0.0;
  309.     if (fabs(targ.im) < TINY)
  310.         targ.im = 0.0;
  311.     *modulus = sqrt((targ.re*targ.re) + (targ.im*targ.im));
  312.     if (targ.im == 0.0)
  313.         *amplitude = 0.0;
  314.     else {
  315.         if (targ.re == 0.0) {
  316.             if (targ.im > 0.0)
  317.                 *amplitude = PID2;
  318.             else
  319.                 *amplitude = -PID2;
  320.             }
  321.         else {
  322.             if ((log(fabs(targ.im)) - log(fabs(targ.re))) > LOGHUGE) {
  323.                 if (targ.re > 0.0) {
  324.                     if (targ.im > 0.0)
  325.                         *amplitude = PID2;
  326.                     else
  327.                         *amplitude = -PID2;
  328.                     }
  329.                 else {
  330.                     if (targ.im > 0.0)
  331.                         *amplitude = -PID2;
  332.                     else
  333.                         *amplitude = PID2;
  334.                     }
  335.                 }
  336.             else
  337.                 *amplitude = atan(targ.im/targ.re);
  338.             }
  339.         }
  340. }
  341.  
  342. ctopower(result,arg,power)
  343. COMPLEX *result,*arg;
  344. int power;
  345.  
  346. /* raises arg to the positive integral power and returns answer in result */
  347.  
  348. {
  349.     int i;
  350.     double modulus,amplitude,newmod,powamp;
  351.     COMPLEX targ;
  352.     
  353.     targ.re = arg->re;
  354.     targ.im = arg->im;
  355.     if (power==0) {
  356.         result->re = 1.0;
  357.         result->im = 0.0;
  358.         }
  359.     else {
  360.         polar(&targ,&modulus,&litude);
  361.         newmod = 1.0;
  362.         if (power > 0)
  363.             for (i=0;i<power;i++)
  364.                 newmod = newmod*modulus;
  365.         else
  366.             for (i=0;i<abs(power);i++)
  367.                 newmod = newmod/modulus;
  368.         powamp = power*amplitude;
  369.         result->re = newmod*cos(powamp);
  370.         result->im = newmod*sin(powamp);
  371.         }
  372. }
  373.  
  374. cexp(result,arg)
  375. COMPLEX *result,*arg;
  376.  
  377. /* raises e to the arg and returns the answer in result */
  378.  
  379. {
  380.     double expo;
  381.     COMPLEX targ;
  382.     
  383.     targ.re = arg->re;
  384.     targ.im = arg->im;
  385.     expo = exp(targ.re);
  386.     result->re = expo*cos(targ.im);
  387.     result->im = expo*sin(targ.im);
  388. }
  389.  
  390. cln(result,arg)
  391. COMPLEX *result,*arg;
  392.  
  393. /* takes the natural logarithm of arg and returns the answer in result */
  394.  
  395. {
  396.     double modulus,amplitude;
  397.     COMPLEX targ;
  398.     
  399.     targ.re = arg->re;
  400.     targ.im = arg->im;
  401.     polar(&targ,&modulus,&litude);
  402.     result->re = log(modulus);
  403.     result->im = amplitude;
  404. }
  405.  
  406. ctoc(result,arg1,arg2)
  407. COMPLEX *result,*arg1,*arg2;
  408.  
  409. /* raise complex number arg1 to complex power arg2 and return answer in result */
  410.  
  411. {
  412.     COMPLEX logpart,expo;
  413.     COMPLEX targ1,targ2;
  414.  
  415.     targ1.re = arg1->re;
  416.     targ1.im = arg1->im;
  417.     targ2.re = arg2->re;
  418.     targ2.im = arg2->im;
  419.     cln(&logpart,&targ1);
  420.     cmult(&expo,&targ2,&logpart);
  421.     cexp(result,&expo);
  422. }
  423.  
  424. csin(result,arg)
  425. COMPLEX *result,*arg;
  426.  
  427. /* takes the sine of arg and returns it in result */
  428.  
  429. {
  430.     COMPLEX targ,exp1,exp2,part1,part2,sum,divisor;
  431.     
  432.     targ.re = arg->re;
  433.     targ.im = arg->im;
  434.     exp1.re = -targ.im;
  435.     exp1.im = targ.re;
  436.     cexp(&part1,&exp1);
  437.     exp2.re = targ.im;
  438.     exp2.im = -targ.re;
  439.     cexp(&part2,&exp2);
  440.     csub(&sum,&part1,&part2);
  441.     divisor.re=0.0;
  442.     divisor.im=2.0;
  443.     cdiv(result,&sum,&divisor);
  444. }
  445.  
  446. ccos(result,arg)
  447. COMPLEX *result,*arg;
  448.  
  449. /* takes the cosine of arg and returns it in result */
  450.  
  451. {
  452.     COMPLEX targ,exp1,exp2,part1,part2,sum,divisor;
  453.     
  454.     targ.re = arg->re;
  455.     targ.im = arg->im;
  456.     exp1.re = -targ.im;
  457.     exp1.im = targ.re;
  458.     cexp(&part1,&exp1);
  459.     exp2.re = targ.im;
  460.     exp2.im = -targ.re;
  461.     cexp(&part2,&exp2);
  462.     cadd(&sum,&part1,&part2);
  463.     divisor.re = 2.0;
  464.     divisor.im = 0.0;
  465.     cdiv(result,&sum,&divisor);
  466. }
  467.  
  468. ----------------------------------- EOF -------------------------------------
  469.  
  470.  
  471.