Path: utzoo!utgpu!news-server.csri.toronto.edu!mailrus!uwm.edu!rpi!zaphod.mps.ohio-state.edu!usc!snorkelwacker!bloom-beacon!eru!luth!sunic!mcsun!ukc!edcastle!hwcs!lauwer From: lauwer@cs.hw.ac.uk (Jean-Marc de Lauwereyns) Newsgroups: comp.graphics Subject: 2D FFT in C : a non optimised source code. Message-ID: <881@odin.cs.hw.ac.uk> Date: 8 May 90 10:00:03 GMT Sender: news@cs.hw.ac.uk Reply-To: lauwer@cs.hw.ac.uk (Jean-Marc de Lauwereyns) Organization: Computer Science, Heriot-Watt U., Scotland Lines: 204 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> 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 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