Bobby@cup.portal.com (Robert Jules Shaughnessy) (06/11/90)
I recieved the fallowing FFT Pascal code in the mail a while
back. It came without any documentation and I dont know how to use
it. I know, this is kinda a stupid question, but does anyone know
what the variables NN and Isign are for? Also, in what form does
the output come, I noticed it uses only one array and not 2 (One
for real and imaginary.)
One might ask, why use this code? This is the smallest FFT I have,
and I want to make a 68000 assembler version of it. If you have a
well documented smaller FFT, PLEASE upload it!
Thanks for your time!
Bobby Shaughnessy
------------------------------ Pascal FFT ------------------------
Program FFT;
Var
array [1000] :doubled;
Procedure FFT( var data : gldarray; nn, isign : integer);
VAR
ii, jj : integer;
mmax : integer;
m, n : integer;
i, j : integer;
istep : integer;
wtemp : double;
wr, wpr : double;
wi, wpi : double;
theta : double;
tempr : double;
tempi : double;
BEGIN
n := 2 * nn;
j := 1;
For ii := 1 to nn do
Begin
i := 2 * ii - 1;
If j > i then
Begin
tempr := data[j];
tempi := data[j+1];
data[j] := data[i];
data[j+1] := data[i+1];
data[i] := tempr;
data[i+1] := tempi;
End;
m := n DIV 2;
While ((m >= 2) and (j > m)) do
Begin
j := j-m;
m := m DIV 2;
End;
j := j + m;
End;
mmax := 2;
While( n > mmax) do
Begin
istep := 2 * mmax;
theta := 2 * Pi / (isign * mmax);
wpr := -2.0 * sqr( sin( 0.5 * theta));
wpi := sin( theta);
wr := 1.0;
wi := 0.0;
For ii := 1 to (mmax DIV 2) do
Begin
m := 2 * ii -1;
For jj := 0 to ((n - m) DIV istep) do
Begin
i := m + jj * istep;
j := i + mmax;
tempr := wr * data[j] - wi * data[j+1];
tempi := wr * data[j+1] + wi * data[j];
data[j] := data[i] - tempr;
data[j+1] := data[i+1] - tempi;
data[i] := data[i] + tempr;
data[i+1] := data[i+1] + tempi;
End;
wtemp := wr;
wr := wr * wpr - wi * wpi + wr;
wi := wi * wpr + wtemp * wpi + wi;
End;
mmax := istep;
End;
End;
Begin (Main)
FFT
End.
paulr@syma.sussex.ac.uk (Paul T Russell) (06/11/90)
This looks very much like the code from "Numerical Recipes in Pascal".
(Vetterling et al ?). The arrays consist of alternate real, complex
values and start at an indexs of 1, so array[1] = x[0].re, array[2] =
x[0].im, array[3] = x[1].re, etc...
Incidentally, I think they have the signs for the FT and inverse FT
the wrong way round - can anyone confirm this ?
//Paul
--
Paul Russell, Department of Experimental Psychology
University of Sussex, Falmer, Brighton BN1 9QG, England
Janet: paulr@uk.ac.sussex.syma Nsfnet: paulr@syma.sussex.ac.uk
Bitnet: paulr%sussex.syma@ukacrl.bitnet Usenet: ...ukc!syma!paulrvlo@xydeco.siemens.com (John Vlontzos) (06/11/90)
Looking at the code: nn is the log(N) (base 2) that is for 1024 points, nn=10 isign controls the sign of exp(2kp/N) (forward/inverse fft) ifft=(1/N)fft(X*) where X*=complex conjugate. I didn't look too hard but I think the routine doesn't divide by N when computing the inverse fft. John Vlontzos Siemens
moyarzun@bbn.com (Miguel Oyarzun) (06/13/90)
In article <2846@syma.sussex.ac.uk> paulr@syma.sussex.ac.uk (Paul T Russell) writes: >This looks very much like the code from "Numerical Recipes in Pascal". >(Vetterling et al ?). The arrays consist of alternate real, complex >values and start at an indexs of 1, so array[1] = x[0].re, array[2] = >x[0].im, array[3] = x[1].re, etc... > >Incidentally, I think they have the signs for the FT and inverse FT >the wrong way round - can anyone confirm this ? > If you look at the Numerical Recipes' definition of the Fourier integral, you will see that their algorithm is consitent with their definition. Unfortunately, they selected the opposite sign convention on the complex exponential from that commonly used in engineering work :-(. I guess it's the old scientist-vs-engineer's way of looking at the world. Miguel