🚀 go-pugleaf

In comp.sys.transputer <1992Nov18.173159.5761@rmece49.upr.clu.edu>
Roque Pagan <rpagan@rmece02.upr.clu.edu> writes:

>	I'm trying to locate some FFT source code in OCCAM.  Any
>    information would be greatly appreciated.

This is from a response I posted a few weeks ago:

... The code consists of three parts: (1) generate the 'twiddle factors' in 
routine Omega,(2) do a fast in-place bit reversal in routine BitReverse,
and (3) do the fft operations.

To use it you need to declare four arrays: the real and imaginary part of your
data of size N; and the real and imaginary parts of the 'twiddle factors' of
size M where 2^M = N; for examples:

  [512]REAL32 DataReal,DataImaginary:
  [9]REAL32 OmegaReal,OmegaImaginary:
  SEQ

    DataReal      :=   ( define or acquire your real data )
    DataImaginary :=   ( define or acquire your imaginary data )

    Fft(DataReal,DataImaginary,OmegaReal,OmegaImaginary)

    -- DataReal and DataImaginary now correspond to the spectrum of the data.

Good Luck!

Jacob Hirbawi
JcbHrb@CERF.net

-- CUT HERE -------------------------------------------------------------------
PROC Omega([]REAL32 w1,w2)
  INT i1,j1 :
  REAL32 a1 :
  SEQ
    i1:=1
    SEQ j1 = 0 FOR SIZE w1
      SEQ
        a1:=3.1415926345 (REAL32) /(REAL32 TRUNC i1)
        w1[j1]:=  COS(a1)
        w2[j1]:= -SIN(a1)
        i1:=i1+i1
:
PROC BitReverse([]REAL32 x1,x2)
  INT i1,i2,i3,j1,j2,n :
  REAL32 t1,t2 :
  SEQ
    n := SIZE x1
    i2:=1
    SEQ i1 = 1 FOR (n-1)
      SEQ
        IF
          i1<i2
            SEQ
              t1:=x1[i1-1]
              t2:=x2[i1-1]
              x1[i1-1]:=x1[i2-1]
              x2[i1-1]:=x2[i2-1]
              x1[i2-1]:=t1
              x2[i2-1]:=t2
          TRUE
            SKIP
        i3 := n/2
        WHILE i3<i2
          SEQ
            i2:=i2-i3
            i3:=i3/2
        i2:=i2+i3
:

PROC Fft([]REAL32 x1,x2,w1,w2)
  INT i1,i2,i3,j1,j2,n,m,k1,k2:
  REAL32 u1,u2,q1,q2,t1,t2:
  SEQ
    Omega(w1,w2)
    BitReverse(x1,x2)
    m := SIZE w1
    n := SIZE x1
    i1:=1
    SEQ j1 = 0 FOR m
      SEQ
        i2:=0
        WHILE i2<n
          SEQ
            u1:= 1.0 (REAL32)
            u2:= 0.0 (REAL32)
            SEQ i3 = i2 FOR i1
              SEQ
                k2:=i3+i1
                q1:=(u1*x1[k2])-(u2*x2[k2])
                q2:=(u2*x1[k2])+(u1*x2[k2])
                x1[k2]:=x1[i3]-q1
                x2[k2]:=x2[i3]-q2
                x1[i3]:=x1[i3]+q1
                x2[i3]:=x2[i3]+q2
                t1:=u1
                t2:=u2
                u1:=(t1*w1[j1])-(t2*w2[j1])
                u2:=(t2*w1[j1])+(t1*w2[j1])
            i2:=i2+(i1+i1)
        i1:=i1+i1
:
-- CUT HERE -------------------------------------------------------------------