(*---------------------------------------------------------------------------
  :Program.     FFT.mod
  :Contents.    Der Fast-Fourier-Transformation Algorithmus
  :Author.      Christian Stiens
  :Address.     Heustiege 2, W-4710 Lüdinghausen
  :Copyright.   PD
  :Language.    Oberon, Assembler-INLINE's
  :Translator.  Amiga Oberon V2.12e [fbs]
  :History.     V1.0, 01-Nov-91
---------------------------------------------------------------------------*)

MODULE FFT;

  (* $ClearVars- $NilChk- $RangeChk- $OvflChk- $ReturnChk- $CaseChk- $StackChk+ *)

  IMPORT
    sys : SYSTEM;

  CONST
    TRAP0 = 04E40H;
    sinTabLen = 1024;

  VAR
    sinTab ["_sintab_1024x16"] : ARRAY sinTabLen OF INTEGER;


  PROCEDURE Sqr (i{0}: INTEGER): LONGINT; (* $EntryExitCode- *)
  BEGIN
    sys.INLINE(0C1C0H,04E75H);   (* MULS D0,D0 *)
  END Sqr;


  PROCEDURE Muls (i{0},j{1}: INTEGER): LONGINT; (* $EntryExitCode- *)
  BEGIN
    sys.INLINE(0C1C1H,04E75H);   (* MULS D1,D0 *)
  END Muls;


  PROCEDURE Mulu (i{0},j{1}: INTEGER): LONGINT; (* $EntryExitCode- *)
  BEGIN
    sys.INLINE(0C0C1H,04E75H);   (* MULU D1,D0 *)
  END Mulu;


  PROCEDURE Sqrt (i{1}: LONGINT): INTEGER; (* $EntryExitCode- *)
  BEGIN
    sys.INLINE(07000H,            (*      MOVEQ  #0,D0     *)
               02201H,            (*      MOVE.L D1,D1     *)
               0671CH,            (*      BEQ.S  ende      *)
               0303CH,08000H,     (*      MOVE.W #32768,D0 *)
               03400H,            (*      MOVE.W D0,D2     *)
               0E24AH,            (*      LSR.W  #1,D2     *)
               03600H,            (* loop MOVE.W D0,D3     *)
               0C6C3H,            (*      MULU   D3,D3     *)
               0B283H,            (*      CMP.L  D3,D1     *)
               06708H,            (*      BEQ.S  \2        *)
               06204H,            (*      BHI.S  \1        *)
               09042H,            (*      SUB.W  D2,D0     *)
               06002H,            (*      BRA.S  \2        *)
               0D042H,            (* \1   ADD.W  D2,D0     *)
               0E24AH,            (* \2   LSR.W  #1,D2     *)
               066ECH,            (*      BNE.S  loop      *)
               04E75H)            (* ende RTS              *)
  END Sqrt;


  PROCEDURE Swap (VAR x,y: INTEGER);
    VAR t: INTEGER;
  BEGIN
    t := x; x := y; y := t;
  END Swap;


  PROCEDURE BitRev (i{1},n{2}: INTEGER): INTEGER;  (* $EntryExitCode- *)
  BEGIN
    sys.INLINE(07000H,        (*      MOVEQ   #0,D0   *)
               05342H,        (*      SUBQ.W  #1,D2   *)
               06B08H,        (*      BMI.S   ende    *)
               0E289H,        (* loop LSR.L   #1,D1   *)
               0E390H,        (*      ROXL    #1,D0   *)
               051CAH,0FFFAH, (*      DBRA    D2,loop *)
               04E75H)        (* ende RTS             *)
  END BitRev;


  PROCEDURE FFT * (VAR xreal,ximag: ARRAY OF INTEGER; n: INTEGER);

    VAR
      log2n,n2,i,l,k,kn2    : INTEGER;
      treal,timag           : LONGINT;
      c,s                   : INTEGER;
      argCos,argSin,argInc  : INTEGER;
      shift                 : BOOLEAN;
      xrkn2,xikn2           : POINTER TO INTEGER;
      xrk,xik               : POINTER TO INTEGER;
      temp                  : LONGINT;

  BEGIN
    IF (n <= 0) OR (n > sinTabLen) THEN sys.INLINE(TRAP0) END;
    log2n := 0;
    WHILE ASH(LONG(1),log2n) < n DO INC(log2n) END;
    IF n # ASH(LONG(1),log2n) THEN sys.INLINE(TRAP0) END;
    IF (LEN(xreal) < n) OR (LEN(ximag) < n) THEN sys.INLINE(TRAP0) END;

    k := n;
    REPEAT
      DEC(k);
      i := BitRev(k,log2n);
      IF i < k THEN
        Swap(xreal[k],xreal[i]);
        Swap(ximag[k],ximag[i]);
      END;
    UNTIL k = 0;

    n2 := 1;
    argInc := sinTabLen;
    shift := TRUE;

    l := 0; WHILE l < log2n DO
      IF l + 1 = log2n THEN shift := FALSE END;
      k := 0;
      argSin := 0;
      argCos := sinTabLen DIV 4;
      argInc := argInc DIV 2;

      REPEAT
        i := n2;
        REPEAT
          s := sinTab[argSin MOD sinTabLen];
          c := sinTab[argCos MOD sinTabLen];

          kn2 := k + n2;

          xrkn2 := sys.ADR(xreal[kn2]);
          xikn2 := sys.ADR(ximag[kn2]);

          xrk := sys.ADR(xreal[k]);
          xik := sys.ADR(ximag[k]);

          treal := ASH(Muls(xrkn2^,c) + Muls(xikn2^,s),-15);
          timag := ASH(Muls(xikn2^,c) - Muls(xrkn2^,s),-15);

          IF shift THEN

            xrkn2^ := SHORT(ASH(LONG(xrk^) - treal,-1));
            xikn2^ := SHORT(ASH(LONG(xik^) - timag,-1));

            xrk^ := SHORT(ASH(LONG(xrk^) + treal,-1));
            xik^ := SHORT(ASH(LONG(xik^) + timag,-1));

          ELSE

            temp := LONG(xrk^) - treal;
            IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
            IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
            xrkn2^ := SHORT(temp);

            temp := LONG(xik^) - timag;
            IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
            IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
            xikn2^ := SHORT(temp);

            temp := LONG(xrk^) + treal;
            IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
            IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
            xrk^ := SHORT(temp);

            temp := LONG(xik^) + timag;
            IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
            IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
            xik^ := SHORT(temp);

          END;

          INC(argSin,argInc);
          INC(argCos,argInc);

          INC(k); DEC(i);
        UNTIL i = 0;

        INC(argSin,sinTabLen DIV 2);
        INC(argCos,sinTabLen DIV 2);

        INC(k,n2);
      UNTIL k >= n;
      INC(n2,n2);
    INC(l) END;
  END FFT;


  PROCEDURE Abs * (VAR xreal,ximag: ARRAY OF INTEGER; n: INTEGER);
  BEGIN
    IF (n > LEN(xreal)) OR (n > LEN(ximag)) THEN
      sys.INLINE(TRAP0)
    END;
    WHILE n > 0 DO
      DEC(n);
      xreal[n] := SHORT(sys.LSH(Mulu(Sqrt(Sqr(xreal[n])+Sqr(ximag[n])),181),-8));
    END;
  END Abs;


END FFT.