Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!mnetor!uunet!seismo!rutgers!labrea!decwrl!pyramid!uccba!hal!ncoast!allbery
From: mh@killer.UUCP (Mike Hobgood)
Newsgroups: comp.sources.misc
Subject: Original fft listing
Message-ID: <4080@ncoast.UUCP>
Date: Thu, 6-Aug-87 21:07:27 EDT
Article-I.D.: ncoast.4080
Posted: Thu Aug  6 21:07:27 1987
Date-Received: Sat, 8-Aug-87 19:05:32 EDT
Sender: allbery@ncoast.UUCP
Organization: The Unix(tm) Connection, Dallas, Texas
Lines: 456
Keywords: Oops...
Approved: allbery@ncoast.UUCP
X-Archive: comp.sources.misc/8708/5



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,&amplitude);
        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,&amplitude);
    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 -------------------------------------