home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d6xx / d666 / spectroscope.lha / Spectroscope / Source / FFT.mod < prev    next >
Text File  |  1992-05-21  |  6KB  |  206 lines
  1. (*---------------------------------------------------------------------------
  2.   :Program.     FFT.mod
  3.   :Contents.    The Fast-Fourier-Transformation Algorithm
  4.   :Author.      Christian Stiens
  5.   :Address.     Heustiege 2, W-4710 Lüdinghausen, GERMANY
  6.   :Copyright.   public domain
  7.   :Language.    Oberon, Assembler-INLINE's
  8.   :Translator.  Amiga Oberon V2.25 (inofficial beta version)
  9.   :History.     V1.0, 01-Nov-91
  10. ---------------------------------------------------------------------------*)
  11.  
  12. MODULE FFT;
  13.  
  14.   (* $ClearVars- $NilChk- $RangeChk- $OvflChk- $ReturnChk- $CaseChk- $StackChk+ *)
  15.  
  16.   IMPORT
  17.     sys : SYSTEM;
  18.  
  19.   CONST
  20.     TRAP0 = 04E40H;
  21.     sinTabLen = 1024;
  22.  
  23.   VAR
  24.     sinTab ["_sintab_1024x16"] : ARRAY sinTabLen OF INTEGER;
  25.  
  26.  
  27.   PROCEDURE Sqr (i{0}: INTEGER): LONGINT; (* $EntryExitCode- *)
  28.   BEGIN
  29.     sys.INLINE(0C1C0H,04E75H);   (* MULS D0,D0 *)
  30.   END Sqr;
  31.  
  32.  
  33.   PROCEDURE Muls (i{0},j{1}: INTEGER): LONGINT; (* $EntryExitCode- *)
  34.   BEGIN
  35.     sys.INLINE(0C1C1H,04E75H);   (* MULS D1,D0 *)
  36.   END Muls;
  37.  
  38.  
  39.   PROCEDURE Mulu (i{0},j{1}: INTEGER): LONGINT; (* $EntryExitCode- *)
  40.   BEGIN
  41.     sys.INLINE(0C0C1H,04E75H);   (* MULU D1,D0 *)
  42.   END Mulu;
  43.  
  44.  
  45.   PROCEDURE Sqrt (i{1}: LONGINT): INTEGER; (* $EntryExitCode- *)
  46.   BEGIN
  47.     sys.INLINE(07000H,            (*      MOVEQ  #0,D0     *)
  48.                02201H,            (*      MOVE.L D1,D1     *)
  49.                0671CH,            (*      BEQ.S  ende      *)
  50.                0303CH,08000H,     (*      MOVE.W #32768,D0 *)
  51.                03400H,            (*      MOVE.W D0,D2     *)
  52.                0E24AH,            (*      LSR.W  #1,D2     *)
  53.                03600H,            (* loop MOVE.W D0,D3     *)
  54.                0C6C3H,            (*      MULU   D3,D3     *)
  55.                0B283H,            (*      CMP.L  D3,D1     *)
  56.                06708H,            (*      BEQ.S  \2        *)
  57.                06204H,            (*      BHI.S  \1        *)
  58.                09042H,            (*      SUB.W  D2,D0     *)
  59.                06002H,            (*      BRA.S  \2        *)
  60.                0D042H,            (* \1   ADD.W  D2,D0     *)
  61.                0E24AH,            (* \2   LSR.W  #1,D2     *)
  62.                066ECH,            (*      BNE.S  loop      *)
  63.                04E75H)            (* ende RTS              *)
  64.   END Sqrt;
  65.  
  66.  
  67.   PROCEDURE Swap (VAR x,y: INTEGER);
  68.     VAR t: INTEGER;
  69.   BEGIN
  70.     t := x; x := y; y := t;
  71.   END Swap;
  72.  
  73.  
  74.   PROCEDURE BitRev (i{1},n{2}: INTEGER): INTEGER;  (* $EntryExitCode- *)
  75.   BEGIN
  76.     sys.INLINE(07000H,        (*      MOVEQ   #0,D0   *)
  77.                05342H,        (*      SUBQ.W  #1,D2   *)
  78.                06B08H,        (*      BMI.S   ende    *)
  79.                0E289H,        (* loop LSR.L   #1,D1   *)
  80.                0E390H,        (*      ROXL    #1,D0   *)
  81.                051CAH,0FFFAH, (*      DBRA    D2,loop *)
  82.                04E75H)        (* ende RTS             *)
  83.   END BitRev;
  84.  
  85.  
  86.   PROCEDURE FFT * (VAR xreal,ximag: ARRAY OF INTEGER; n: INTEGER);
  87.  
  88.     VAR
  89.       log2n,n2,i,l,k,kn2    : INTEGER;
  90.       treal,timag           : LONGINT;
  91.       c,s                   : INTEGER;
  92.       argCos,argSin,argInc  : INTEGER;
  93.       shift                 : BOOLEAN;
  94.       xrkn2,xikn2           : POINTER TO INTEGER;
  95.       xrk,xik               : POINTER TO INTEGER;
  96.       temp                  : LONGINT;
  97.  
  98.   BEGIN
  99.     IF (n <= 0) OR (n > sinTabLen) THEN sys.INLINE(TRAP0) END;
  100.     log2n := 0;
  101.     WHILE ASH(LONG(1),log2n) < n DO INC(log2n) END;
  102.     IF n # ASH(LONG(1),log2n) THEN sys.INLINE(TRAP0) END;
  103.     IF (LEN(xreal) < n) OR (LEN(ximag) < n) THEN sys.INLINE(TRAP0) END;
  104.  
  105.     k := n;
  106.     REPEAT
  107.       DEC(k);
  108.       i := BitRev(k,log2n);
  109.       IF i < k THEN
  110.         Swap(xreal[k],xreal[i]);
  111.         Swap(ximag[k],ximag[i]);
  112.       END;
  113.     UNTIL k = 0;
  114.  
  115.     n2 := 1;
  116.     argInc := sinTabLen;
  117.     shift := TRUE;
  118.  
  119.     l := 0; WHILE l < log2n DO
  120.       IF l + 1 = log2n THEN shift := FALSE END;
  121.       k := 0;
  122.       argSin := 0;
  123.       argCos := sinTabLen DIV 4;
  124.       argInc := argInc DIV 2;
  125.  
  126.       REPEAT
  127.         i := n2;
  128.         REPEAT
  129.           s := sinTab[argSin MOD sinTabLen];
  130.           c := sinTab[argCos MOD sinTabLen];
  131.  
  132.           kn2 := k + n2;
  133.  
  134.           xrkn2 := sys.ADR(xreal[kn2]);
  135.           xikn2 := sys.ADR(ximag[kn2]);
  136.  
  137.           xrk := sys.ADR(xreal[k]);
  138.           xik := sys.ADR(ximag[k]);
  139.  
  140.           treal := ASH(Muls(xrkn2^,c) + Muls(xikn2^,s),-15);
  141.           timag := ASH(Muls(xikn2^,c) - Muls(xrkn2^,s),-15);
  142.  
  143.           IF shift THEN
  144.  
  145.             xrkn2^ := SHORT(ASH(LONG(xrk^) - treal,-1));
  146.             xikn2^ := SHORT(ASH(LONG(xik^) - timag,-1));
  147.  
  148.             xrk^ := SHORT(ASH(LONG(xrk^) + treal,-1));
  149.             xik^ := SHORT(ASH(LONG(xik^) + timag,-1));
  150.  
  151.           ELSE
  152.  
  153.             temp := LONG(xrk^) - treal;
  154.             IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
  155.             IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
  156.             xrkn2^ := SHORT(temp);
  157.  
  158.             temp := LONG(xik^) - timag;
  159.             IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
  160.             IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
  161.             xikn2^ := SHORT(temp);
  162.  
  163.             temp := LONG(xrk^) + treal;
  164.             IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
  165.             IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
  166.             xrk^ := SHORT(temp);
  167.  
  168.             temp := LONG(xik^) + timag;
  169.             IF temp > MAX(INTEGER) THEN temp := MAX(INTEGER) END;
  170.             IF temp < MIN(INTEGER) THEN temp := MIN(INTEGER) END;
  171.             xik^ := SHORT(temp);
  172.  
  173.           END;
  174.  
  175.           INC(argSin,argInc);
  176.           INC(argCos,argInc);
  177.  
  178.           INC(k); DEC(i);
  179.         UNTIL i = 0;
  180.  
  181.         INC(argSin,sinTabLen DIV 2);
  182.         INC(argCos,sinTabLen DIV 2);
  183.  
  184.         INC(k,n2);
  185.       UNTIL k >= n;
  186.       INC(n2,n2);
  187.     INC(l) END;
  188.   END FFT;
  189.  
  190.  
  191.   PROCEDURE Abs * (VAR xreal,ximag: ARRAY OF INTEGER; n: INTEGER);
  192.   BEGIN
  193.     IF (n > LEN(xreal)) OR (n > LEN(ximag)) THEN
  194.       sys.INLINE(TRAP0)
  195.     END;
  196.     WHILE n > 0 DO
  197.       DEC(n);
  198.       xreal[n] := SHORT(sys.LSH(Mulu(Sqrt(Sqr(xreal[n])+Sqr(ximag[n])),181),-8));
  199.     END;
  200.   END Abs;
  201.  
  202.  
  203. END FFT.
  204.  
  205.  
  206.