/*
FFTEASY consists of the four C functions fftc1, fftcn, fftr1, and
fftrn. FFTEASY is free. I am not in any way, shape, or form expecting
to make money off of these routines. I wrote them because I needed
them for some work I was doing and I'm putting them out on the
Internet in case other people might find them useful. Feel free to
download them, incorporate them into your code, modify them, translate
the comment lines into Swahili, or whatever else you want. What I do
want is the following:
1) Leave this notice (i.e. this entire paragraph beginning with
``FFTEASY consists of...'' and ending with my email address) in with
the code wherever you put it. Even if you're just using it in-house in
your department, business, or wherever else I would like these credits
to remain with it. This is partly so that people can...
2) Give me feedback. Did FFTEASY work great for you and help your
work? Did you hate it? Did you find a way to improve it, or translate
it into another programming language? Whatever the case might be, I
would love to hear about it. Please let me know at the email address
below.
3) Finally, insofar as I have the legal right to do so I forbid you
to make money off of this code without my consent. In other words if
you want to publish these functions in a book or bundle them into
commercial software or anything like that contact me about it
first. I'll probably say yes, but I would like to reserve that right.
For any comments or questions you can reach me at
gfelder@physics.stanford.edu.
*/
/* These declarations are put here so you can easily cut and paste them into your program. */
void fftc1(float f[], int N, int skip, int forward);
void fftcn(float f[], int ndims, int size[], int forward);
void fftr1(float f[], int N, int forward);
void fftrn(float f[], float fnyquist[], int ndims, int size[], int forward);
#include
#include
#include
struct complex
{
float real;
float imag;
};
/*
Do a Fourier transform of an array of N complex numbers separated by steps of (complex) size skip.
The array f should be of length 2N*skip and N must be a power of 2.
Forward determines whether to do a forward transform (1) or an inverse one(-1)
*/
void fftc1(float f[], int N, int skip, int forward)
{
int b,index1,index2,trans_size,trans;
float pi2 = 4.*asin(1.);
float pi2n,cospi2n,sinpi2n; /* Used in recursive formulas for Re(W^b) and Im(W^b) */
struct complex wb; /* wk = W^k = e^(2 pi i b/N) in the Danielson-Lanczos formula for a transform of length N */
struct complex temp1,temp2; /* Buffers for implementing recursive formulas */
struct complex *c = (struct complex *)f; /* Treat f as an array of N complex numbers */
/* Place the elements of the array c in bit-reversed order */
for(index1=1,index2=0;index1=b;b/=2) /* To find the next bit reversed array index subtract leading 1's from index2 */
index2-=b;
index2+=b; /* Next replace the first 0 in index2 with a 1 and this gives the correct next value */
if(index2>index1) /* Swap each pair only the first time it is found */
{
temp1 = c[index2*skip];
c[index2*skip] = c[index1*skip];
c[index1*skip] = temp1;
}
}
/* Next perform successive transforms of length 2,4,...,N using the Danielson-Lanczos formula */
for(trans_size=2;trans_size<=N;trans_size*=2) /* trans_size = size of transform being computed */
{
pi2n = forward*pi2/(float)trans_size; /* +- 2 pi/trans_size */
cospi2n = cos(pi2n); /* Used to calculate W^k in D-L formula */
sinpi2n = sin(pi2n);
wb.real = 1.; /* Initialize W^b for b=0 */
wb.imag = 0.;
for(b=0;b=0;dim--) /* Loop over dimensions */
{
planesize *= size[dim]; /* Planesize = Product of all sizes up to and including size[dim] */
for(i=0;i=0) or an inverse one(<0)
*/
void fftr1(float f[], int N, int forward)
{
int b;
float pi2n = 4.*asin(1.)/N,cospi2n=cos(pi2n),sinpi2n=sin(pi2n); /* pi2n = 2 Pi/N */
struct complex wb; /* wb = W^b = e^(2 pi i b/N) in the Danielson-Lanczos formula for a transform of length N */
struct complex temp1,temp2; /* Buffers for implementing recursive formulas */
struct complex *c = (struct complex *)f; /* Treat f as an array of N/2 complex numbers */
if(forward==1)
fftc1(f,N/2,1,1); /* Do a transform of f as if it were N/2 complex points */
wb.real = 1.; /* Initialize W^b for b=0 */
wb.imag = 0.;
for(b=1;b=0;j--) /* If the rightmost indices are maximal reset them to 0. Indexneg goes from 1 to 0 in these dimensions. */
{
indices[j]=0;
indexneg -= stepsize;
stepsize *= size[j];
}
if(indices[j]==0) /* If index[j] goes from 0 to 1 indexneg[j] goes from 0 to size[j]-1 */
indexneg += stepsize*(size[j]-1);
else /* Otherwise increasing index[j] decreases indexneg by one unit. */
indexneg -= stepsize;
if(j>=0) /* This avoids writing outside the array bounds on the last pass through the array loop */
indices[j]++;
} /* End of i loop (over total array) */
if(forward==-1) /* Inverse transform */
fftcn(f,ndims,size,-1);
size[ndims-1] *= 2; /* Give the user back the array size[] in its original condition */
free(indices); /* Free up memory allocated for indices array */
}