🚀 go-pugleaf

RetroBBS NetNews Server

Inspired by RockSolid Light RIP Retro Guy

4 total messages Started by Robert Jules Shaughnessy Sun, 10 Jun 1990 23:09
Mystery FFT Code...
#813
Author: Robert Jules Shaughnessy
Date: Sun, 10 Jun 1990 23:09
109 lines
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...
#820
Author: Paul T Russell
Date: Mon, 11 Jun 1990 14:09
15 lines
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...
#815
Author: John Vlontzos
Date: Mon, 11 Jun 1990 14:24
9 lines
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...
#826
Author: Miguel Oyarzun
Date: Tue, 12 Jun 1990 21:05
17 lines
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