/*
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.
*/

/* changes by J.Janacek (www.biomed.cas.cz/~janacek):
a)complex was renamed to ecomplex to avoid collision with VC++ complex from <math.h>,
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)
f)procedure frmultexp for shifting the origin to middle of the real array in frequency domain
*/

/* 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);
void frmultexp(float idat[], int ndims, int size[]);

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

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 < N; index1++) { /* Loop through all elements of c */
for(b= N / 2; index2 >= 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 < trans_size / 2; b++) { /* Step over half of the elements in the transform */
for(trans= 0; trans < N / trans_size; trans++) { /* Iterate over all transforms of size trans_size to be computed */
index1= (trans * trans_size + b) * skip; /* Index of element in first half of transform being computed */
index2= index1 + trans_size / 2 * skip; /* Index of element in second half of transform being computed */
temp1= c[index1];
temp2= c[index2];
c[index1].real= temp1.real + wb.real * temp2.real - wb.imag * temp2.imag; /* Implement D-L formula */
c[index1].imag= temp1.imag + wb.real * temp2.imag + wb.imag * temp2.real;
c[index2].real= temp1.real - wb.real * temp2.real + wb.imag * temp2.imag;
c[index2].imag= temp1.imag - wb.real * temp2.imag - wb.imag * temp2.real;
}
temp1= wb;
wb.real= cospi2n * temp1.real - sinpi2n * temp1.imag; /* Real part of eˆ(2 pi i b/trans_size) used in D-L formula */
wb.imag= cospi2n * temp1.imag + sinpi2n * temp1.real; /* Imaginary part of eˆ(2 pi i b/trans_size) used in D-L formula */
}
}

/* For an inverse transform divide by the number of grid points */
if(forward<0.)
for(index1= 0; index1 < skip * N; index1+= skip)
{
c[index1].real/= N;
c[index1].imag/= N;
}
}

/*
Do a Fourier transform of an ndims dimensional array of ecomplex numbers
Array dimensions are given by size[0], ..., size[ndims-1]. Note that these are sizes of ecomplex arrays.
The array f should be of length 2*size[0]*...*size[ndims-1] and all sizes must be powers of 2.
Forward determines whether to do a forward transform (1) or an inverse one(-1)
*/
void fftcnd(float f[], int ndims, int size[], int forward) {

int i, j, dim;
int planesize= 1, skip= 1; /* These determine where to begin successive transforms and the skip between their elements (see below) */
int totalsize= 1; /* Total size of the ndims dimensional array */

for(dim= 0; dim < ndims; dim++) /* Determine total size of array */
totalsize*= size[dim];

for(dim= 0; dim < ndims; dim++) { /* Loop over dimensions */
planesize*= size[dim]; /* Planesize = Product of all sizes up to and including size[dim] */
for(i= 0; i < totalsize; i+= planesize) /* Take big steps to begin loops of transforms */
for(j= 0; j < skip; j++) /* Skip sets the number of transforms in between big steps as well as the skip between elements */
fftc1(f + 2 * (i + j), size[dim], skip, forward); /* 1-D Fourier transform. (Factor of two converts ecomplex index to float index.) */
skip *= size[dim]; /* Skip = Product of all sizes up to (but not including) size[dim] */
}
}

/*
Do a Fourier transform of an array of N real numbers
N must be a power of 2
Forward determines whether to do a forward transform (>= 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 < N / 4; b++) { /* Loop over elements of transform. See documentation for these formulas */
temp1= wb;
wb.real= cospi2n * temp1.real - sinpi2n * temp1.imag; /* Real part of eˆ(2 pi i b/N) used in D-L formula */
wb.imag= cospi2n * temp1.imag + sinpi2n * temp1.real; /* Imaginary part of eˆ(2 pi i b/N) used in D-L formula */
temp1= c[b];
temp2= c[N / 2 - b];
c[b].real= .5 * (temp1.real + temp2.real + forward * wb.real * (temp1.imag + temp2.imag) + wb.imag * (temp1.real - temp2.real));
c[b].imag= .5 * (temp1.imag - temp2.imag - forward * wb.real * (temp1.real - temp2.real) + wb.imag * (temp1.imag + temp2.imag));
c[N / 2 - b].real= .5*( temp1.real + temp2.real - forward * wb.real * (temp1.imag + temp2.imag) - wb.imag * (temp1.real - temp2.real));
c[N / 2 - b].imag= .5*(-temp1.imag + temp2.imag - forward * wb.real * (temp1.real - temp2.real) + wb.imag * (temp1.imag + temp2.imag));
}
temp1= c[0];
c[0].real= temp1.real + temp1.imag; /* Set b= 0 term in transform */
c[0].imag= temp1.real - temp1.imag; /* Put b= N/2 term in imaginary part of first term */

if(forward == -1)
{
c[0].real*= .5;
c[0].imag*= .5;
fftc1(f, N / 2, 1, -1);
}
}

void fftrnd(float f[], int ndims, int isize[], int forward) {

int i, j, b;
int index, indexneg= 0; /* Positions in the 1-d arrays of points labeled by indices (i0, i1, ..., i(ndims-1));
indexneg gives the position in the array of the corresponding negative frequency */
int stepsize; // Used in calculating indexneg
int N= isize[0] / 2; /* The size of the FIRST dimension is used often enough to merit its own name. */
double pi2n= 2. * 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 array of ecomplex numbers */
int totalsize= 1; /* Total number of ecomplex points in array */
int *indices= (int *) malloc(ndims * sizeof(int)); /* Indices for looping through array */
int *size= (int *) malloc(ndims * sizeof(int)); /* Indices for looping through array */
if(!indices || !size) { /* Make sure memory was correctly allocated */
printf("Error allocating memory in fftrnd routine. Exiting.\n");
exit(1);
}

for(i= 0; i<ndims; i++) {
size[i]= isize[i];
totalsize *= size[i];
indices[i] = 0;
}
totalsize/= 2;
size[0]/= 2; /* Set size[] to be the sizes of f viewed as a ecomplex array */
if(forward == 1) /* Forward transform */
fftcnd(f, ndims, size, 1); /* Do a transform of f as if it were N/2 ecomplex points */
for(index= 0; index < totalsize; index+= N) { /* Loop over all but first array index */
wb.real= 1.; /* Initialize Wˆb for b= 0 */
wb.imag= 0.;
for(b= 1; b < N / 2; b++) { /* Loop over elements of transform. See documentation for these formulas */
temp1= wb;
wb.real= cospi2n * temp1.real - sinpi2n * temp1.imag; /* Real part of eˆ(2 pi i b/N_real) used in D-L formula */
wb.imag= cospi2n * temp1.imag + sinpi2n * temp1.real; /* Imaginary part of eˆ(2 pi i b/N_real) used in D-L formula */
temp1= c[index + b];
temp2= c[indexneg + N - b]; /* Note that N-b is NOT the negative frequency for b. Only nonnegative b momenta are stored. */
c[index + b].real= .5 * (temp1.real + temp2.real + forward * wb.real * (temp1.imag + temp2.imag) + wb.imag * (temp1.real - temp2.real));
c[index + b].imag= .5 * (temp1.imag - temp2.imag - forward * wb.real * (temp1.real - temp2.real) + wb.imag * (temp1.imag + temp2.imag));
c[indexneg + N - b].real= .5*( temp1.real + temp2.real - forward * wb.real * (temp1.imag + temp2.imag) - wb.imag * (temp1.real - temp2.real));
c[indexneg + N - b].imag= .5*(-temp1.imag + temp2.imag - forward * wb.real * (temp1.real - temp2.real) + wb.imag * (temp1.imag + temp2.imag));
}
if (index < indexneg) {
temp1= c[index];
temp2= c[indexneg];
c[index].real= .5 * (temp1.real+temp2.real + forward * (temp1.imag + temp2.imag));
c[index].imag= .5 * (temp1.imag-temp2.imag - forward * (temp1.real - temp2.real));
c[indexneg].real= .5 * ( temp1.real + temp2.real - forward * (temp1.imag + temp2.imag));
c[indexneg].imag= .5 * (-temp1.imag + temp2.imag - forward * (temp1.real - temp2.real));
}
else if (index == indexneg) {
temp1= c[index];
c[index].real= temp1.real + temp1.imag;
c[index].imag= temp1.real - temp1.imag;
if(forward == -1) {
c[index].real*= .5;
c[index].imag*= .5;
}
}

/* Find indices for positive and single index for negative frequency.
In each dimension indexneg[j]= 0 if index[j]= 0, indexneg[j]= size[j]-index[j] otherwise. */
stepsize= N; /* Amount to increment indexneg by as each individual index is incremented */
for(j= 1; j < ndims && indices[j] == size[j]-1; j++) { /* If the leftmost 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(j < ndims) { /* This avoids writing outside the array bounds on the last pass through the array loop */
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;
indices[j]++;
}
} /* End of index loop (over total array) */

if(forward == -1) /* Inverse transform */
fftcnd(f, ndims, size, -1);

free(indices); /* Free up memory allocated for indices array */
free(size); /* Free up memory allocated for indices array */
}

void frmult(float i1[], float i2[], float o[], int ndims, int size[], int conj) {
/*multiplies two Fourrier images of nD real matrices to calculate convolution (conj= = 0)
or correlation (conj= = 1) in frquency domain*/
float *p1, *p2, *q, rp1, rp2;
long i, j, k, n, ok, N, M;
long indexneg= 0; /* Positions in the 1-d arrays of points labeled by indices (i0, i1, ..., i(ndims-1));
indexneg gives the position in the array of the corresponding negative frequency */
long stepsize; // Used in calculating indexneg
long *indices;

indices= (long *)malloc(ndims * sizeof(long)); /* Indices for looping through array */
if(!indices) { /* Make sure memory was correctly allocated */
printf("allocating memory\n");
return;
}
for (j= 0, n= 1; j < ndims; j++) {
indices[j]= 0;
n*= size[j];
}

#define CONJUG \
rp1= *p1++; rp2= *p2++; \
*(q++)= rp1 * rp2 + *p1 * *p2; \
*(q++)= *p1 * rp2 - rp1 * *p2; \
*p1++; *p2++;

#define SIMPLE \
rp1= *p1++; rp2= *p2++; \
*(q++)= rp1 * rp2 - *p1 * *p2; \
*(q++)= *p1 * rp2 + rp1 * *p2; \
*p1++; *p2++;

#define INDEX \
stepsize= N; \
for(j= 1; j < ndims && indices[j] == size[j] - 1; j++) \
{ \
indices[j]= 0; \
indexneg-= stepsize; \
stepsize*= size[j]; \
} \
if(j < ndims) {\
if(indices[j] == 0) \
indexneg+= stepsize*(size[j] - 1); \
else \
indexneg-= stepsize; \
indices[j]++; \
}

#define DOCYC(OPER) \
for (i= 0; i < n; i+= M) { \
p1= i1 + i; \
p2= i2 + i; \
q= o + i; \
if (i != indexneg) { \
OPER \
} \
else { \
*(q++)= *(p1++) * *(p2++); \
*(q++)= *(p1++) * *(p2++); \
} \
for (k= 1; k < M/2; k++) { \
OPER \
} \
INDEX \
}

q= o; p1= i1; p2= i2;
M= size[0];
N= 2 * M;

if (conj)
DOCYC(CONJUG)
else
DOCYC(SIMPLE)

free(indices);

#undef INDEX
#undef CONJUG
#undef SIMPLE
#undef DOCYC
}

void frmultexp(float idat[], int ndims, int size[]) {
/*multiplies Fourrier image of nD real matrice to move 0 position to center
in spatial domain*/
float *q, f;
long i, j, k, n, ok, N, M;
long indexneg= 0; /* Positions in the 1-d arrays of points labeled by indices (i0, i1, ..., i(ndims-1));
indexneg gives the position in the array of the corresponding negative frequency */
long stepsize; // Used in calculating indexneg
long *indices;

indices= (long *)malloc(ndims * sizeof(long)); /* Indices for looping through array */
if(!indices) { /* Make sure memory was correctly allocated */
printf("allocating memory\n");
return;
}
for (j= 0, n= 1; j < ndims; j++) {
indices[j]= 0;
n*= size[j];
}

M= size[0];
N= 2 * M;

f= 1.0;

for (i= 0; i < n; i+=M) {
q= idat + i;
*(q++)*= f; *(q++)*= f;
f= -f;
for (k= 1; k < M/2; k++, f= -f) {
*(q++)*= f; *(q++)*= f;
}
stepsize=N;
for(j=1;j<ndims && indices[j]==size[j]-1;j++) {
indices[j]=0;
stepsize *= size[j];
f= -f;
}
if(j<ndims) indices[j]++;
f= -f;
}

free(indices);
}