mh@killer.UUCP (Mike Hobgood) (08/01/87)
[]
I found this file on one of my disks, maybe useful.
#!/bin/sh
# This is a shell archive. Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by mh on Fri Jul 31 00:44:01 CDT 1987
# Contents: fft.c complex.c
echo x - fft.c
sed 's/^XX//' > "fft.c" <<'@//E*O*F fft.c//'
XX/*
XX * fft.c
XX *
XX * Performs a Cooley-Tukey Fast Fourier Transform.
XX *
XX * Original version Pascal written in "Simple Calculations with
XX * Complex Numbers by David Clark Doctor Dobbs Journal 10/84
XX *
XX * Translated to C by R. Hellman 02/21/86
XX * Released to Public Domain
XX *
XX * Difference from original version:
XX *
XX * As in fortran versions that I have worked with, instead of
XX * specifying total array size in decimal - instead specify number
XX * of array members as power of two. Forward or reverse fft is
XX * specified by the sign of this argument.
XX *
XX * ****IMPORTANT****
XX * For convenience sake members of an array are taken to begin at [1]
XX * rather than [0] thus for an fft on a 64 member array of pointers
XX * to complex numbers, there must be 65 members.
XX *
XX * f is a one dimensional array of pointers to complex numbers with
XX * length (2**abs(ln)). The sign of ln determines in which direction
XX * the transform will be performed. If ln is positive, then isinverse
XX * is set to -1 and a forward transform is carried out. If ln is
XX * negative then isinverse is set to 1 and a reverse transform is
XX * calculated.
XX *
XX * For numpoints = (2**abs(ln)),
XX * transform[j] = sum(f[i] * w ^ ((i-1) * (j-1))),
XX * where i, j are 1 to numpoints and
XX * w = exp(isinverse * 2 * pi * sqrt(-1)/numpoints).
XX *
XX * The program also computes the inverse transform, for which the
XX * defining expression is:
XX * inverse DFT = (1/numpoints) * sum(f[i] * w ^ ((i-1) * (j-1))).
XX *
XX * Run time is proportional to numpoints * log2(numpoints) rather than
XX * to numpoints ** 2 for the classical discrete Fourier transform.
XX */
XX#include <stdio.h>
XX#include "complex.c"
XXfft(f, ln)
XXcomplex *f[];
XXint ln;
XX{
XX int i, isinverse, numpoints;
XX if (ln > 0) {
XX puts("Calculating forward Fourier transform\n");
XX isinverse = -1;
XX numpoints = pow(2.0, ln);
XX }
XX else {
XX puts("Calculating reverse Fourier transform\n");
XX isinverse = 1;
XX numpoints = (int) pow(2.0, -ln);
XX }
XX bitreverse(numpoints, f);
XX butterfly(numpoints, isinverse, f);
XX}
XX/*
XX * Reverse the bits
XX *
XX */
XXstatic bitreverse(numpoints, f)
XXint numpoints;
XXcomplex *f[];
XX{
XX int i, j, m;
XX complex *temp;
XX j = 1;
XX for (i = 1; i <= numpoints; i++) {
XX if (i < j) {
XX temp = f[j];
XX f[j] = f[i];
XX f[i] = temp;
XX }
XX m = numpoints >> 1;
XX if (m < j)
XX while (m < j) {
XX j -= m;
XX m = (m + 1) >> 1;
XX }
XX j += m;
XX }
XX}
XX/*
XX * Calculate the butterflies for the bit reversed data
XX * Normalize if a reverse fft is performed
XX *
XX */
XXstatic butterfly(numpoints, isinverse, f)
XXint numpoints, isinverse;
XXcomplex *f[];
XX{
XX int mmax, step, index, m, i, j;
XX double theta;
XX complex w, cn, temp;
XX mmax = 1;
XX while (numpoints > mmax) {
XX step = mmax << 1;
XX for (m = 1; m <= mmax; m++) {
XX theta = M_PI * ((double)(isinverse) *
XX (double)(m - 1)) / (double)(mmax);
XX w.re = cos(theta);
XX w.im = sin(theta);
XX for (i = 1; i <= ((numpoints-m) / step) + 1; i++) {
XX index = m + ((i-1) * step);
XX j = index + mmax;
XX cmult(&temp, &w, f[j]);
XX csub(f[j], f[index], &temp);
XX cadd(f[index], f[index], &temp);
XX }
XX }
XX mmax = step;
XX }
XX if (isinverse == 1) {
XX cn.re = numpoints;
XX cn.im = 0.0;
XX for (i = 1; i <= numpoints; i++)
XX cdiv(f[i], f[i], &cn);
XX }
XX}
@//E*O*F fft.c//
chmod u=rw,g=r,o=r fft.c
echo x - complex.c
sed 's/^XX//' > "complex.c" <<'@//E*O*F complex.c//'
XX/*
XX * Complex number library translated from Pascal source listing
XX * in Dr. Dobbs October, 1984 "Simple Calculations with Complex Numbers"
XX * by David Clark.
XX *
XX * To insure that the same structure could be used as both an 'argument'
XX * and 'result' immediate tranfer of structure members is done to a
XX * local variable - in many of the functions this may not be necessary,
XX * but it was done in all of the following functions for consistency
XX *
XX * Author: R. Hellman
XX * Date: 02/21/86
XX *
XX * Released to the Public Domain
XX */
XX#include <malloc.h>
XX#include <math.h>
XX#define TINY 2.2e-308 /* tiny value */
XX#define LOGHUGE 709.778 /* natural log of huge value */
XXtypedef struct {
XX double re;
XX double im;
XX} complex;
XX/*
XX * Initializes storage for each address of an array of complex pointers
XX * *f[] with 'arraysz' number of elements
XX */
XXinitcmplx(f, arraysz)
XXcomplex *f[];
XXint arraysz;
XX{
XX int i;
XX for (i = 0; i < arraysz; i++)
XX if ((f[i] = (complex *)malloc(sizeof(complex))) == NULL) {
XX puts("Error: complex number allocation\n");
XX exit(1);
XX }
XX}
XX/* adds arg1 and arg2 and returns sum in result */
XXcadd(result, arg1, arg2)
XXcomplex *result, *arg1, *arg2;
XX{
XX complex targ1, targ2;
XX targ1.re = arg1->re;
XX targ1.im = arg1->im;
XX targ2.re = arg2->re;
XX targ2.im = arg2->im;
XX result->re = targ1.re + targ2.re;
XX result->im = targ1.im + targ2.im;
XX}
XX/* subtracts arg1 and arg2 and returns sum in result */
XXcsub(result, arg1, arg2)
XXcomplex *result, *arg1, *arg2;
XX{
XX complex targ1, targ2;
XX targ1.re = arg1->re;
XX targ1.im = arg1->im;
XX targ2.re = arg2->re;
XX targ2.im = arg2->im;
XX result->re = targ1.re - targ2.re;
XX result->im = targ1.im - targ2.im;
XX}
XX/* multiplies arg1 and arg2 and returns product in result */
XXcmult(result, arg1, arg2)
XXcomplex *result, *arg1, *arg2;
XX{
XX complex targ1, targ2;
XX targ1.re = arg1->re;
XX targ1.im = arg1->im;
XX targ2.re = arg2->re;
XX targ2.im = arg2->im;
XX result->re = (targ1.re * targ2.re) - (targ1.im * targ2.im);
XX result->im = (targ1.im * targ2.re) + (targ1.re * targ2.im);
XX}
XX/* divides arg1 and arg2 and returns quotient in result */
XXcdiv(result, arg1, arg2)
XXcomplex *result, *arg1, *arg2;
XX{
XX double denom;
XX complex targ1,targ2;
XX targ1.re = arg1->re;
XX targ1.im = arg1->im;
XX targ2.re = arg2->re;
XX targ2.im = arg2->im;
XX denom = (targ2.re * targ2.re) + (targ2.im * targ2.im);
XX result->re = (targ1.re * targ2.re + targ1.im * targ2.im) / denom;
XX result->im = (targ1.im * targ2.re - targ1.re * targ2.im) / denom;
XX}
XX/* converts a complex number arg in rectangular form to polar form */
XXpolar(arg, modulus, amplitude)
XXcomplex *arg;
XXdouble *modulus, *amplitude;
XX{
XX complex targ;
XX targ.re = arg->re;
XX targ.im = arg->im;
XX if (fabs(targ.re) < TINY)
XX targ.re = 0.0;
XX if (fabs(targ.im) < TINY)
XX targ.im = 0.0;
XX *modulus = sqrt((targ.re * targ.re) + (targ.im * targ.im));
XX if (targ.im == 0.0)
XX *amplitude = 0.0;
XX else if (targ.re == 0.0) {
XX if (targ.im > 0.0)
XX *amplitude = M_PI_2;
XX else
XX *amplitude = -M_PI_2;
XX }
XX else if ((log(fabs(targ.im)) - log(fabs(targ.re))) > LOGHUGE) {
XX if (targ.re > 0.0) {
XX if (targ.im > 0.0)
XX *amplitude = M_PI_2;
XX else
XX *amplitude = -M_PI_2;
XX }
XX else if (targ.im > 0.0)
XX *amplitude = -M_PI_2;
XX else
XX *amplitude = M_PI_2;
XX }
XX }
XX else
XX *amplitude = atan(targ.im / targ.re);
XX}
XX/* raises arg to the positive integral power and returns answer in result */
XXctopower(result, arg, power)
XXcomplex *result, *arg;
XXint power;
XX{
XX int i;
XX double modulus, amplitude, newmod, powamp;
XX complex targ;
XX targ.re = arg->re;
XX targ.im = arg->im;
XX if (power == 0) {
XX result->re = 1.0;
XX result->im = 0.0;
XX }
XX else {
XX polar(&targ, &modulus, &litude);
XX newmod = 1.0;
XX if (power > 0)
XX for (i = 0; i < power; i++)
XX newmod = newmod * modulus;
XX else
XX for (i = 0; i < abs(power); i++)
XX newmod = newmod / modulus;
XX powamp = power * amplitude;
XX result->re = newmod * cos(powamp);
XX result->im = newmod * sin(powamp);
XX }
XX}
XX/* raises e to the arg and returns the answer in result */
XXcexp(result, arg)
XXcomplex *result, *arg;
XX{
XX double expo;
XX complex targ;
XX targ.re = arg->re;
XX targ.im = arg->im;
XX expo = exp(targ.re);
XX result->re = expo * cos(targ.im);
XX result->im = expo * sin(targ.im);
XX}
XX/* takes the natural logarithm of arg and returns the answer in result */
XXcln(result, arg)
XXcomplex *result, *arg;
XX{
XX double modulus, amplitude;
XX complex targ;
XX
XX targ.re = arg->re;
XX targ.im = arg->im;
XX polar(&targ, &modulus, &litude);
XX result->re = log(modulus);
XX result->im = amplitude;
XX}
XX/* raise complex number arg1 to complex power arg2 and return answer in result */
XXctoc(result, arg1, arg2)
XXcomplex *result, *arg1, *arg2;
XX{
XX complex logpart, expo, targ1, targ2;
XX targ1.re = arg1->re;
XX targ1.im = arg1->im;
XX targ2.re = arg2->re;
XX targ2.im = arg2->im;
XX cln(&logpart, &targ1);
XX cmult(&expo, &targ2, &logpart);
XX cexp(result, &expo);
XX}
XX/* takes the sine of arg and returns it in result */
XXcsin(result, arg)
XXcomplex *result, *arg;
XX{
XX complex targ, exp1, exp2, part1,
XX part2, sum, divisor;
XX targ.re = arg->re;
XX targ.im = arg->im;
XX exp1.re = -targ.im;
XX exp1.im = targ.re;
XX cexp(&part1, &exp1);
XX exp2.re = targ.im;
XX exp2.im = -targ.re;
XX cexp(&part2, &exp2);
XX csub(&sum, &part1, &part2);
XX divisor.re = 0.0;
XX divisor.im = 2.0;
XX cdiv(result, &sum, &divisor);
XX}
XX/* takes the cosine of arg and returns it in result */
XXccos(result, arg)
XXcomplex *result, *arg;
XX{
XX complex targ, exp1, exp2, part1,
XX part2, sum, divisor;
XX targ.re = arg->re;
XX targ.im = arg->im;
XX exp1.re = -targ.im;
XX exp1.im = targ.re;
XX cexp(&part1, &exp1);
XX exp2.re = targ.im;
XX exp2.im = -targ.re;
XX cexp(&part2, &exp2);
XX cadd(&sum, &part1, &part2);
XX divisor.re = 2.0;
XX divisor.im = 0.0;
XX cdiv(result, &sum, &divisor);
@//E*O*F complex.c//
chmod u=rw,g=r,o=r complex.c
echo Inspecting for damage in transit...
temp=/tmp/sharin$$; dtemp=/tmp/sharout$$
trap "rm -f $temp $dtemp; exit" 0 1 2 3 15
cat > $temp <<\!!!
152 569 3300 fft.c
306 872 5702 complex.c
458 1441 9002 total
!!!
wc fft.c complex.c | sed 's=[^ ]*/==' | diff -b $temp - >$dtemp
if test -s $dtemp
then echo "Ouch [diff of wc output]:" ; cat $dtemp
else echo "No problems found."
fi
exit 0