<F0O@psuvm.psu.edu> (05/13/91)
I'd like to thank everyone who send a reply on my question of
finding an fft program in C!
Below are several fft programs. Enjoy!
[Tim]
=========================================================================
Date: Thu, 9 May 91 09:33:48 GMT-0800
From: duchow@watnxt3.ucr.edu (John Duchowski)
To: F0O@psuvm.psu.edu
Subject: Re: Would like C source for an FFT program!
Newsgroups: comp.sys.ibm.pc.misc,comp.lang.c
Organization: University of California, Riverside
In article <91129.104628F0O@psuvm.psu.edu> you write:
>
> Could anyone send me C source code, or know an FTP site where
>there is source code for an FFT program/function?
>
> Thanks!
>
> [Tim]
Send the following one-line command to listserv@ucf1vm:
index turboc-l
I believe they have something called fft there which comes xxencoded. However,
they also have a whole bunch of other stuff which you might find useful. If
you have never dealt with a listserv before, send one-line command:
help
to the same address. Some of the files may be Turbo-C specific but may only
require little readjustement to work on other systems. Good Luck !
- John
=========================================================================
Date: Fri, 10 May 91 13:12:24 EDT
From: johns@calvin.ee.cornell.edu (John Sahr)
To: F0O@psuvm.psu.edu
Subject: Would like C source for an FFT program!
Could anyone send me C source code, or know an FTP site where
there is source code for an FFT program/function?
Thanks!
[Tim]
Try these if you wish. These are some 1d and multidimensional
in-place complex ffts, single and double precision. You may copy
the code, and freely distribute it, but you may not sell it. You
may not restrict the distribution of the source, a la GNU.
These algorithms are based on Numerical Recipes in C code; I have
recoded them, because NRiC code is horrid.
There are 3 files here; complex.h, fatal-error.h, and fft.c. The two
include files are needed by fft.c
-john
--------- complex.h (abridged to FFT stuff) -------------------------------
/* jdsahr 10 jan 1989 */
/* by John Sahr, Cornell Space Plasma Physics Group */
/* 607 255 8298 johns@magneto.ee.cornell.edu */
#ifndef __complex_h__
#define __complex_h__
#include "math.h"
#ifndef M_PI
#define M_PI 3.14159265358979
#endif
typedef struct {double re, im;} complex;
typedef struct {float re, im;} complex_float;
/* from fft.c */
extern void cfft_nD(); /* multidimensional complex fft */
extern void cfft_stride(); /* 1D fft with stride */
extern void cfft(); /* 1D fft on close-packed data */
/* single precision versions of the above */
extern void cfft_nD_f(); /* multidimensional complex fft */
extern void cfft_stride_f(); /* 1D fft with stride */
extern void cfft_f(); /* 1D fft on close-packed data */
#endif __complex_h_
--------- fatal-error.h --------------------------------------------------
/* fatal-error.h, jdsahr 10jan1989, 6feb1989 */
#ifndef fatal_error
#include "stdio.h"
extern void exit();
#define fatal_error(ex,note) {if (ex){fprintf(stderr,"fatal error in %s:%s at
line %d\n--> %s\n",__FILE__, __SUBNAME__, __LINE__, note);exit(1);}}
#define fatal(note) {fprintf(stderr,"fatal error in %s:%s at line %d\n-->
%s\n",__FILE__, __SUBNAME__, __LINE__, note);exit(1);}
#ifndef __SUBNAME__
#define __SUBNAME__ " \0"
#endif __SUBNAME__
#endif
/***************************************************************************
These macro are nice error handlers built on the model of assert(3X).
If "ex" evaluates to nonzero, then a message is printed to stderr, and
exit(2) is called. The __FILE__ and __LINE__ symbols are defined at
compile time, and the __SUBNAME__ symbol may be set to the subroutine
name (or any other string desired, but its purpose is to hold the
subroutine name).
The "ifndef lint" switch is to keep lint from complaining that the
integer value returned by fprintf() is never used. It is pointless
and difficult to handle it, as the program is about to exit().
***************************************************************************/
#ifdef DEBUG
#define Dprintf(a) fprintf a
#else
#define Dprintf(a)
#endif
/***************************************************************************
The Dprintf(a) macro gets turned on if DEBUG is defined. Usage:
Dprintf((stderr,"the time is %s",now));
Note that the double parens are necessary.
***************************************************************************/
--------- fft.c ----------------------------------------------------
/* $Id$ */
/* by John Sahr, Cornell Space Plasma Physics Group */
/* 607 255 8298 johns@magneto.ee.cornell.edu */
#include "stdio.h"
#include "fatal-error.h"
#include "complex.h"
/* This function provides a pointer to a table for bit reversing. The
integer passed must be a power of 2, as these subroutines are intended
strictly for radix 2 FFTs. */
static int *get_bit_reverse_table(old_table,n)
int *old_table, n;
{
static int last_n = -1;
static int *last_old_table = NULL;
int i, j, top;
if(n == last_n && old_table == last_old_table) return old_table;
if(old_table)
{
old_table = (int *) realloc((char *) old_table,n*(sizeof (int)));
fatal_error(old_table == NULL,"realloc, malloc(3C) failed");
}
else
{
old_table = (int *) malloc(n*(sizeof (int)));
fatal_error(old_table == NULL,"malloc(3C) failed");
}
last_old_table = old_table;
last_n = n;
old_table[0] = 0;
for(top = 1; top < n; top <<= 1)
{
for(i = 0, j = top; i < top; i++, j++)
{
old_table[i] <<= 1;
old_table[j] = old_table[i]|1;
}
}
fatal_error(top != n,"n must be a power of 2");
return(old_table);
}
/* cfft() is a routine for calculating 1D ffts. The data must be compact
in memory.
data points to the complex vector
n length of vector. must be power of 2 (error in get_bit_rev...)
sense direction of transform: (sense > 0) ? forward : reverse
sense == 0 is indeterminate, and considered to be an error
If the data is not closely packed, use cfft_stride(). */
void cfft(data,n,sense)
complex *data;
int n, sense;
{
complex z, w, wp, *ai, *aj;
double theta;
int step, bigstep, i, ii, j, jj, m;
static int *swap_table = (int *) NULL;
fatal_error(data == NULL,"storage must be allocated for data");
fatal_error(n <= 0,"the length of the vector must be greater than zero");
fatal_error(sense == 0,"the sense parameter must be positive or negative");
swap_table = get_bit_reverse_table(swap_table,n);
for(i = 0; i < n; i++)
{
if(i < (j = swap_table[i]))
{
z = data[i];
data[i] = data[j];
data[j] = z;
}
}
step = 1;
bigstep = 2;
theta = (sense > 0) ? M_PI : -M_PI;
while(step < n)
{
wp.re = cos(theta); wp.im = sin(theta);
w.re = 1.0; w.im = 0.0;
for(m = 0; m < step; m++)
{
for(i = m; i < n; i += bigstep)
{
ai = data + i;
aj = ai + step;
z.re = w.re*aj->re - w.im*aj->im;
z.im = w.re*aj->im + w.im*aj->re;
aj->re = ai->re - z.re;
aj->im = ai->im - z.im;
ai->re = ai->re + z.re;
ai->im = ai->im + z.im;
}
z.re = w.re;
z.im = w.im;
w.re = z.re*wp.re - w.im*wp.im;
w.im = z.im*wp.re + z.re*wp.im;
}
step = bigstep;
bigstep *= 2;
theta *= 0.5;
}
}
/* cfft_stride() is a routine for calculating 1D ffts. By supplying a
"stride" parameter, you can do an fft on "sparse" data, when the vector
is not close packed.
data points to the complex vector
n length of vector. must be power of 2 (error in get_bit_rev...)
stride memory spacing of adjacent vector elements. May be negative if
you want to scan the vector in the opposite direction (why?)
sense direction of transform: (sense > 0) ? forward : reverse
sense == 0 is indeterminate, and considered to be an error */
void cfft_stride(data,n,stride,sense)
complex *data;
int n, stride, sense;
{
complex z, w, wp, *ai, *aj;
double theta;
int step, bigstep, i, ii, j, jj, m;
static int *swap_table = (int *) NULL;
fatal_error(data == NULL,"storage must be allocated for data");
fatal_error(n <= 0,"the length of the vector must be greater than zero");
fatal_error(sense == 0,"the sense parameter must be positive or negative");
swap_table = get_bit_reverse_table(swap_table,n);
for(i = 0; i < n; i++)
{
if((ii = i) < (jj = swap_table[i]))
{
ii *= stride;
jj *= stride;
z = data[ii];
data[ii] = data[jj];
data[jj] = z;
}
}
step = 1;
bigstep = 2;
theta = (sense > 0) ? M_PI : -M_PI;
while(step < n)
{
wp.re = cos(theta); wp.im = sin(theta);
w.re = 1.0; w.im = 0.0;
for(m = 0; m < step; m++)
{
for(i = m, j = m + step; i < n; i += bigstep, j += bigstep)
{
ai = data + i*stride;
aj = ai + step*stride;
z.re = w.re*aj->re - w.im*aj->im;
z.im = w.re*aj->im + w.im*aj->re;
aj->re = ai->re - z.re;
aj->im = ai->im - z.im;
ai->re = ai->re + z.re;
ai->im = ai->im + z.im;
}
z.re = w.re;
z.im = w.im;
w.re = z.re*wp.re - w.im*wp.im;
w.im = z.im*wp.re + z.re*wp.im;
}
step = bigstep;
bigstep *= 2;
theta *= 0.5;
}
}
/* The following procedure calculates the forward or backward in-place
multidimensional FFT, using the factoring method of Numerical
Recipes in C. The ghost of NRiC may be noted in the code, but
I have renamed (usefully) almost all of the variables, and
have implemented a different bit reverse scheme.
data pointer to complex array
nn integer array containing dimensions
ndim (integer) number of dimension
sense flag indicating direction of FFT
(sense > 0) ? forward transform : inverse transform
sense == 0 is considered to be indeterminate, a fatal error.
Note that the inverse transform is not scaled. This can be
done by calling cfft_scale() or cfft_scale_nD().
nn[0] contains the dimension which is scanned fastest, followed
by nn[1], nn[2], ... In other words, (fake C syntax)
data[i][j][k] = data[i + nn[0]*(j + nn[1]*k)]
Note that it is allowed to compute 1-D FFTs with this routine,
if ndim = 1 and nn[0] = the length of the array. */
void cfft_nD(data, nn, ndim, sense)
complex *data;
int *nn, ndim, sense;
{
static int *swap_table = (int *) NULL;
int lower, end_lower, step_lower;
int upper, end_upper, step_upper;
int stride_index, stride, total_elements, which_dim, dimension;
int i, m, ilow, irev, bigstep, step;
complex *pt, *ptrev, *p1, *p2;
complex temp, w, wp;
double signed_pi, theta;
fatal_error(data == NULL,"no storage supplied for data");
fatal_error(nn == NULL,"dimensions must be supplied");
fatal_error(ndim <= 0,"the number of dimensions must be positive");
fatal_error(sense == 0,"the sense paramater must be positive or negative");
for (total_elements = 1, which_dim = 0; which_dim < ndim; which_dim++)
{
total_elements *= nn[which_dim];
}
signed_pi = (sense > 0) ? M_PI : -M_PI;
end_lower = 1;
step_lower = 1;
end_upper = total_elements;
step_upper = 1;
for(which_dim = 0; which_dim < ndim; which_dim++)
{
dimension = nn[which_dim];
swap_table = get_bit_reverse_table(swap_table,dimension);
fatal_error(swap_table == NULL,"get_bit_reverse_table() failed");
stride = end_lower;
step_upper *= dimension;
for(upper = 0; upper < end_upper; upper += step_upper)
{
for(lower = 0; lower < end_lower; lower += step_lower)
{
pt = data + upper + lower;
for(i = 0; i < dimension; i++)
{
irev = swap_table[i];
if(irev < i)
{
irev *= stride;
ilow = i*stride;
temp = pt[ilow];
pt[ilow] = pt[irev];
pt[irev] = temp;
}
}
theta = signed_pi;
step = 1;
bigstep = 2;
while(step < dimension)
{
wp.re = cos(theta); wp.im = sin(theta);
w.re = 1.0; w.im = 0.0;
for(m = 0; m < step; m++)
{
for(i = m; i < dimension; i += bigstep)
{
p1 = pt + stride*i;
p2 = p1 + stride*step;
temp.re = w.re*p2->re - w.im*p2->im;
temp.im = w.re*p2->im + w.im*p2->re;
p2->re = p1->re - temp.re;
p2->im = p1->im - temp.im;
p1->re = p1->re + temp.re;
p1->im = p1->im + temp.im;
}
temp.re = w.re;
temp.im = w.im;
w.re = temp.re*wp.re - temp.im*wp.im;
w.im = temp.im*wp.re + temp.re*wp.im;
}
step = bigstep;
bigstep *= 2;
theta *= 0.5;
}
}
}
end_lower = step_upper;
}
}
/**************************************************************************
Single precision complex routines (expect pointers to *fcomplex);
**************************************************************************/
/* cfft_f() is a routine for calculating 1D ffts. The data must be compact
in memory.
data points to the fcomplex vector
n length of vector. must be power of 2 (error in get_bit_rev...)
sense direction of transform: (sense > 0) ? forward : reverse
sense == 0 is indeterminate, and considered to be an error
If the data is not closely packed, use cfft_stride_f(). */
void cfft_f(data,n,sense)
fcomplex *data;
int n, sense;
{
fcomplex z, w, wp, *ai, *aj;
double theta;
int step, bigstep, i, ii, j, jj, m;
static int *swap_table = (int *) NULL;
fatal_error(data == NULL,"storage must be allocated for data");
fatal_error(n <= 0,"the length of the vector must be greater than zero");
fatal_error(sense == 0,"the sense parameter must be positive or negative");
swap_table = get_bit_reverse_table(swap_table,n);
for(i = 0; i < n; i++)
{
if(i < (j = swap_table[i]))
{
z = data[i];
data[i] = data[j];
data[j] = z;
}
}
step = 1;
bigstep = 2;
theta = (sense > 0) ? M_PI : -M_PI;
while(step < n)
{
wp.re = cos(theta); wp.im = sin(theta);
w.re = 1.0; w.im = 0.0;
for(m = 0; m < step; m++)
{
for(i = m; i < n; i += bigstep)
{
ai = data + i;
aj = ai + step;
z.re = w.re*aj->re - w.im*aj->im;
z.im = w.re*aj->im + w.im*aj->re;
aj->re = ai->re - z.re;
aj->im = ai->im - z.im;
ai->re = ai->re + z.re;
ai->im = ai->im + z.im;
}
z.re = w.re;
z.im = w.im;
w.re = z.re*wp.re - w.im*wp.im;
w.im = z.im*wp.re + z.re*wp.im;
}
step = bigstep;
bigstep *= 2;
theta *= 0.5;
}
}
/* cfft_stride() is a routine for calculating 1D ffts. By supplying a
"stride" parameter, you can do an fft on "sparse" data, when the vector
is not close packed.
data points to the fcomplex vector
n length of vector. must be power of 2 (error in get_bit_rev...)
stride memory spacing of adjacent vector elements. May be negative if
you want to scan the vector in the opposite direction (why?)
sense direction of transform: (sense > 0) ? forward : reverse
sense == 0 is indeterminate, and considered to be an error */
void cfft_stride_f(data,n,stride,sense)
fcomplex *data;
int n, stride, sense;
{
fcomplex z, w, wp, *ai, *aj;
double theta;
int step, bigstep, i, ii, j, jj, m;
static int *swap_table = (int *) NULL;
fatal_error(data == NULL,"storage must be allocated for data");
fatal_error(n <= 0,"the length of the vector must be greater than zero");
fatal_error(sense == 0,"the sense parameter must be positive or negative");
swap_table = get_bit_reverse_table(swap_table,n);
for(i = 0; i < n; i++)
{
if((ii = i) < (jj = swap_table[i]))
{
ii *= stride;
jj *= stride;
z = data[ii];
data[ii] = data[jj];
data[jj] = z;
}
}
step = 1;
bigstep = 2;
theta = (sense > 0) ? M_PI : -M_PI;
while(step < n)
{
wp.re = cos(theta); wp.im = sin(theta);
w.re = 1.0; w.im = 0.0;
for(m = 0; m < step; m++)
{
for(i = m, j = m + step; i < n; i += bigstep, j += bigstep)
{
ai = data + i*stride;
aj = ai + step*stride;
z.re = w.re*aj->re - w.im*aj->im;
z.im = w.re*aj->im + w.im*aj->re;
aj->re = ai->re - z.re;
aj->im = ai->im - z.im;
ai->re = ai->re + z.re;
ai->im = ai->im + z.im;
}
z.re = w.re;
z.im = w.im;
w.re = z.re*wp.re - w.im*wp.im;
w.im = z.im*wp.re + z.re*wp.im;
}
step = bigstep;
bigstep *= 2;
theta *= 0.5;
}
}
/* The following procedure calculates the forward or backward in-place
multidimensional FFT, using the factoring method of Numerical
Recipes in C. The ghost of NRiC may be noted in the code, but
I have renamed (usefully) almost all of the variables, and
have implemented a different bit reverse scheme.
data pointer to fcomplex array
nn integer array containing dimensions
ndim (integer) number of dimension
sense flag indicating direction of FFT
(sense > 0) ? forward transform : inverse transform
sense == 0 is considered to be indeterminate, a fatal error.
Note that the inverse transform is not scaled. This can be
done by calling cfft_scale() or cfft_scale_nD().
nn[0] contains the dimension which is scanned fastest, followed
by nn[1], nn[2], ... In other words, (fake C syntax)
data[i][j][k] = data[i + nn[0]*(j + nn[1]*k)]
Note that it is allowed to compute 1-D FFTs with this routine,
if ndim = 1 and nn[0] = the length of the array. */
void cfft_nD_f(data, nn, ndim, sense)
fcomplex *data;
int *nn, ndim, sense;
{
static int *swap_table = (int *) NULL;
int lower, end_lower, step_lower;
int upper, end_upper, step_upper;
int stride_index, stride, total_elements, which_dim, dimension;
int i, m, ilow, irev, bigstep, step;
fcomplex *pt, *ptrev, *p1, *p2;
fcomplex temp, w, wp;
double signed_pi, theta;
fatal_error(data == NULL,"no storage supplied for data");
fatal_error(nn == NULL,"dimensions must be supplied");
fatal_error(ndim <= 0,"the number of dimensions must be positive");
fatal_error(sense == 0,"the sense paramater must be positive or negative");
for (total_elements = 1, which_dim = 0; which_dim < ndim; which_dim++)
{
total_elements *= nn[which_dim];
}
signed_pi = (sense > 0) ? M_PI : -M_PI;
end_lower = 1;
step_lower = 1;
end_upper = total_elements;
step_upper = 1;
for(which_dim = 0; which_dim < ndim; which_dim++)
{
dimension = nn[which_dim];
swap_table = get_bit_reverse_table(swap_table,dimension);
fatal_error(swap_table == NULL,"get_bit_reverse_table() failed");
stride = end_lower;
step_upper *= dimension;
for(upper = 0; upper < end_upper; upper += step_upper)
{
for(lower = 0; lower < end_lower; lower += step_lower)
{
pt = data + upper + lower;
for(i = 0; i < dimension; i++)
{
irev = swap_table[i];
if(irev < i)
{
irev *= stride;
ilow = i*stride;
temp = pt[ilow];
pt[ilow] = pt[irev];
pt[irev] = temp;
}
}
theta = signed_pi;
step = 1;
bigstep = 2;
while(step < dimension)
{
wp.re = cos(theta); wp.im = sin(theta);
w.re = 1.0; w.im = 0.0;
for(m = 0; m < step; m++)
{
for(i = m; i < dimension; i += bigstep)
{
p1 = pt + stride*i;
p2 = p1 + stride*step;
temp.re = w.re*p2->re - w.im*p2->im;
temp.im = w.re*p2->im + w.im*p2->re;
p2->re = p1->re - temp.re;
p2->im = p1->im - temp.im;
p1->re = p1->re + temp.re;
p1->im = p1->im + temp.im;
}
temp.re = w.re;
temp.im = w.im;
w.re = temp.re*wp.re - temp.im*wp.im;
w.im = temp.im*wp.re + temp.re*wp.im;
}
step = bigstep;
bigstep *= 2;
theta *= 0.5;
}
}
}
end_lower = step_upper;
}
}
=========================================================================
Date: Fri, 10 May 91 16:44:28 EST
From: milind@agni (Milind Saraph)
To: F0O@psuvm.psu.edu
Subject: Re: Would like C source for an FFT program!
Here it is:
/*
fft - calculate FFT
Carl Crawford
Purdue University
W. Lafayette, IN. 47907
Modified by J. P. Allebach 11 May 1983
Calling Sequence....fft(real,im,m,iopt)
Where real and im are the real and imaginary
parts of the input data. The result is
returned in place. M is the log base 2
of the number of elements in the array.
Call with iopt = 0 for DFT and iopt = 1 for
inverse DFT.
*/
#include <math.h>
fft(a,b,m,iopt)
double a[]; /* real part of data */
double b[]; /* imaginary part of data */
int m; /* size of data = 2**m */
int iopt; /* 0=dft, 1=idft */
{
int nv2,nm1,n,le,le1,ip;
double pi,pile1,tmp;
double ua,ub,wa,wb,ta,tb,*ap,*bp;
register i,j,l;
extern double sin(), cos();
n = 1<<m;
if(iopt){
for(i=0,ap=a,bp=b;i<n;i++){
*ap++ /= n;
*bp++ /= -n;
}
}
nv2 = n/2;
nm1 = n - 1;
j = 0;
for(i=0;i<nm1;i++){
if(i<j){
ta = a[j]; tb = b[j];
a[j] = a[i]; b[j] = b[i];
a[i] = ta; b[i] = tb;
}
l = nv2;
while(l < (j+1) ){
j = j - l;
l = l / 2;
}
j = j + l;
}
pi = 3.1415926535;
for(l=1;l<=m;l++){
le = 1<<l;
le1 = le>>1;
ua = 1.0; ub = 0.0;
pile1 = pi / le1;
wa = cos(pile1); wb = -sin(pile1);
for(j=0;j<le1;j++){
for(i=j;i<n;i += le){
ip = i + le1;
ta = a[ip] * ua - b[ip] * ub;
tb = a[ip] * ub + b[ip] * ua;
a[ip] = a[i] - ta;
b[ip] = b[i] - tb;
a[i] += ta;
b[i] += tb;
}
ua = (tmp = ua) * wa - ub * wb;
ub = tmp * wb + ub * wa;
}
}
if(iopt != 0){
for(i=0;i<n;i++)
b[i] = -b[i];
}
}
=========================================================================
Date: Fri, 10 May 91 22:07:28 edt
From: Chris Sherman <sherman@unx.sas.com>
To: F0O@psuvm.psu.edu
Subject: Re: Would like C source for an FFT program!
Newsgroups: comp.sys.ibm.pc.misc,comp.lang.c
In comp.lang.c you write:
> Could anyone send me C source code, or know an FTP site where
>there is source code for an FFT program/function?
> Thanks!
> [Tim]
This got posted to the net a while ago:
==================
>Hi,
>I am looking for C source for Fast Fourier Transform (any dimensions and
>power of 2). The netlib only has the Fortran version and I do not want
>to use f2c.
>If you have a good FFT program (fast!), that is excellent.
>Thank you.
>Joseph Yip
Maybe this will help.
tom
Tom Friedel JDyx Enterprises (404) 320-7624 tpf@jdyx.UUCP
Unix BBS: (404) 325-1719 <= 2400 ; (404) 321-5020 >= 2400
"Live simply, so that others may simply live."
Well, according to the amount of E-mail I received, and as I cannot
answer to all of them, I am posting the 2D FFT program I have written in C.
How to use it : first, it can be used only on matrix of the type
complexe, as defined at the begining of the program, with a size in
2^maxlevel (2 power maxlevel). One just have to call the program as following :
InvFFT(matrix,maxlevel);
As you have noticed, it is in fact an Inverse Fourier Transform
(sampling matrix in the frequency space -> matrix of heights) which
is performed, but just an initialisation line at the begining of the
procedure InvFFT has to be change to obtain the FFT :
Omega = 2 * M_PI / dimension; must be change into
Omega = - 2 * M_PI / dimension;
I am clearly aware that this program has not been optimised as it
could be, mainly in order to keep the the clarity of the algorithm, but I am
open to any suggestion that you could make in order to improve the process.
I have tested it in a fractal landscape project and it seems to work
correctly (when there were bugs in it such as forgetting the multiplication
of omega by two at the recursive calls, I had funny landscapes looking like
a micro-chip seen through an electronic microscop).
Well, I think that is all, for the moment. Here comes the program.
Note : you will have probably to cut the signature at the end of the
file as well.
Jean-Marc
----cut here -----cut here -----cut here -----cut here -----cut here -----
/* It is a pity I have to add this notice, but it must be done. */
/* NOTICE *
* Copyright (c) 1990, Jean-Marc de Lauwereyns *
* for the procedures mirror, RecFFT, InvFFT. *
* These procedures can be used and given freely to anybody but can not *
* be sold or used in any commercial program or in any book without the *
* permission of the original author. *
* These procedures can be modified at the condition that the nature of *
* modification is mentioned in comments at the place of the modification *
* with the original lines accompanied by the name of the original author.*
* This notice itself can not be removed or modified except by the *
* original author himself. */
/* you must have included the file /usr/include/math.h at the begining of *
* your program */
typedef struct {
double a,b;
} complexe;
/****************************************************/
/* function mirror just transform the binary */
/* representation of an integer into its mirror */
/* binary representation on maxlevel bits */
/* eg : with maxlevel = 6, 000101 becomes 101000 */
/****************************************************/
int mirror(index,maxlevel)
int index;
int maxlevel;
{
register int i;
int res;
res = 0;
for (i=0;i<maxlevel;i++)
{
res = (res << 1) | (index & (int)1);
index = (index >> 1);
}
return(res);
}
/*****************************************************/
/* recursive procedure to calculate the FFT on 2^N */
/*****************************************************/
void RecFFT(A,k,l,p,q,Omega)
complexe **A;
int k,l; /* indexes for the begin and the end on the lines taken */
int p,q; /* indexes for the begin and the end of the columns taken */
double Omega;
{
int i,j,m;
complexe a,b,c,d;
double argument[3];
if (k != l)
{
m = (l-k)/2;
for (i = k;i < k + m;i++)
{
argument[0] = Omega * (double)(i-k);
for (j = p;j < p + m;j++)
{
argument[1] = Omega * (double)(j-p);
argument[2] = Omega * (double)(i-k + j-p);
a.a = A[i][j].a;a.b = A[i][j].b;
b.a = A[i+m][j].a;b.b = A[i+m][j].b;
c.a = A[i][j+m].a;c.b = A[i][j+m].b;
d.a = A[i+m][j+m].a;d.b = A[i+m][j+m].b;
A[i][j].a = a.a + b.a + c.a + d.a;
A[i][j].b = a.b + b.b + c.b + d.b;
A[i+m][j].a = (a.a - b.a + c.a - d.a)*cos(argument[0]) -
(a.b - b.b + c.b - d.b)*sin(argument[0]);
A[i+m][j].b = (a.a - b.a + c.a - d.a)*sin(argument[0]) +
(a.b - b.b + c.b - d.b)*cos(argument[0]);
A[i][j+m].a = (a.a + b.a - c.a - d.a)*cos(argument[1]) -
(a.b + b.b - c.b - d.b)*sin(argument[1]);
A[i][j+m].b = (a.a + b.a - c.a - d.a)*sin(argument[1]) +
(a.b + b.b - c.b - d.b)*cos(argument[1]);
A[i+m][j+m].a = (a.a - b.a - c.a + d.a)*cos(argument[2]) -
(a.b - b.b - c.b + d.b)*sin(argument[2]);
B[i+m][j+m].b = (a.a - b.a - c.a + d.a)*sin(argument[2]) +
(a.b - b.b - c.b + d.b)*cos(argument[2]);
}
}
if (m != 0)
{
RecFFT(k,k+m,p,p+m,Omega*2);
RecFFT(k+m,l,p,p+m,Omega*2);
RecFFT(k,k+m,p+m,q,Omega*2);
RecFFT(k+m,l,p+m,q,Omega*2);
}
}
}
/********************************************************/
/* InvFFT is a procedure which produce the inverse of */
/* Fourier Transformation. */
/********************************************************/
void InvFFT(A,interlevel)
complexe **A;
int interlevel;
{
int dimension;
int i,j,i1,j1;
double coeff,Omega;
dimension = (1 << interlevel);
Omega = 2 * M_PI / dimension;
RecFFT(A,0,dimension,0,dimension,Omega);
for (i=0;i<=(dimension - 1);i++)
/* swap every rows i with its mirror image */
{
i1 = mirror(i,interlevel);
if (i1 > i) /* if i1 <= i, the swaaping have been already made */
for (j=0;j<dimension;j++)
{
coeff = A[i][j].a;
A[i][j].a = A[i1][j].a;
A[i1][j].a = coeff;
/* this may be remove if you are sure that the result is a *
* matrix of real which means that the original matrix must *
* verify the following conditions : *
* ______ ______ *
* A(N-i,N-j) = A(i,j) , A(0,N-j) = A(0,j) *
* ______ *
* A(N-i,0) = A(i,0) for i,j > 0 and A(0,0) is a pure *
* real */
coeff = A[i][j].b;
A[i][j].a = A[i1][j].a;
A[i1][j].a = coeff;
}
}
for (j=0;j<=(dimension - 1);j++) /* idem for the columns */
{
j1 = mirror(j,interlevel);
if (j1 >j)
for (i=0;i<dimension;i++)
{
coeff = A[i][j].a;
A[i][j].a = A[i][j1].a;
A[i][j1].a = coeff;
/* part which may be removed if (see above) */
coeff = A[i][j].b;
A[i][j].b = A[i][j1].b;
A[i][j1].b = coeff;
}
}
}
----cut here -----cut here -----cut here -----cut here -----cut here ----
Jean-Marc de Lauwereyns | ____ | e-mail addresses :
Heriot-Watt University | |\ /| | \ | JANET: lauwer@uk.ac.hw.cs
Computer Science Department | | \/ | |___/ | ARPA.: lauwer@cs.hw.ac.uk
Edinburgh | \___/ |___| \ |