Thread View: comp.dsp
4 messages
4 total messages
Started by Robert Jules Shaughnessy
Sun, 10 Jun 1990 23:09
Mystery FFT Code...
Author: Robert Jules Shaughnessy
Date: Sun, 10 Jun 1990 23:09
Date: Sun, 10 Jun 1990 23:09
109 lines
2877 bytes
2877 bytes
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.
Re: Mystery FFT Code...
Author: Paul T Russell
Date: Mon, 11 Jun 1990 14:09
Date: Mon, 11 Jun 1990 14:09
15 lines
970 bytes
970 bytes
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!paulr
Re: Mystery FFT Code...
Author: John Vlontzos
Date: Mon, 11 Jun 1990 14:24
Date: Mon, 11 Jun 1990 14:24
9 lines
739 bytes
739 bytes
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
Re: Mystery FFT Code...
Author: Miguel Oyarzun
Date: Tue, 12 Jun 1990 21:05
Date: Tue, 12 Jun 1990 21:05
17 lines
1336 bytes
1336 bytes
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
Thread Navigation
This is a paginated view of messages in the thread with full content displayed inline.
Messages are displayed in chronological order, with the original post highlighted in green.
Use pagination controls to navigate through all messages in large threads.
Back to All Threads