/*
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.
*/
/*principal changes by J.Janacek (www.biomed.cas.cz/~janacek):
a)complex was renamed to ecomplex to avoid collision with VC++ complex from ,
b)the order of dimension was changed to be the same as in IRIS Explorer, leftmost dimension is now first,
c)real FFT was changed to work at place
d)added the frmult procedure for calculation of convolution or correlation in frequency domain
e)the nd procedures were renamed, to avoid possible collisions caused by b),c)
November 6: bug in fftrnd: initialization of pi2n
*/
/* 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 fftcnd(float f[], int ndims, int size[], int forward);
void fftr1(float f[], int N, int forward);
void fftrnd(float f[], int ndims, int size[], int forward);
void frmult(float i1[], float i2[], float o[], int ndims, int size[], int conj);
#include
#include
#include
struct ecomplex
{
float real;
float imag;
};
/*
Do a Fourier transform of an array of N ecomplex numbers separated by steps of (ecomplex) 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 ecomplex wb; /* wk = W^k = e^(2 pi i b/N) in the Danielson-Lanczos formula for a transform of length N */
struct ecomplex temp1,temp2; /* Buffers for implementing recursive formulas */
struct ecomplex *c = (struct ecomplex *)f; /* Treat f as an array of N ecomplex 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) 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 ecomplex wb; /* wb = W^b = e^(2 pi i b/N) in the Danielson-Lanczos formula for a transform of length N */
struct ecomplex temp1,temp2; /* Buffers for implementing recursive formulas */
struct ecomplex *c = (struct ecomplex *)f; /* Treat f as an array of N/2 ecomplex numbers */
if(forward==1)
fftc1(f,N/2,1,1); /* Do a transform of f as if it were N/2 ecomplex points */
wb.real = 1.; /* Initialize W^b for b=0 */
wb.imag = 0.;
for(b=1;b