lauwer@cs.hw.ac.uk (Jean-Marc de Lauwereyns) (05/09/90)
There is an error, that I haven't noticed. Mea culpa : I just added
the increminated part especially for the posting I made. It is in the lines
interlacing in the procedure InvFFT (exchanges between line i and line i1),
when it comes to process the imaginary part of the matrix. Those lines have
been added especially for the posting I made and I haven't verified them
carefully enough.
It was just two lines to coorect, but the bug led to desastrous
consequences for the result.
Once again, excuse me : nobody is perfect.
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;
/* erroneous lines : *
* A[i][j].a = A[i1][j].a; *
* A[i1][j].a = coeff; *
* it should be *.b and not *.a */
A[i][j].b = A[i1][j].b;
A[i1][j].b = 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 | \___/ |___| \ |bcs86318@zach.fit.edu ( CRAIG) (05/10/90)
Just wanted to thank you for posting the code for the 2 dimensional FFT. It might come in handy for a project I'm working on. Thanks again.