home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 202.lha / FFT / fft.jf < prev    next >
Encoding:
Text File  |  1988-12-28  |  13.9 KB  |  481 lines
  1. include? sqrt ju:sqrt  ( fast and speedy integer square root)
  2. include? X* 4th:fft/X*.ASM  ( assembler complex multiply)
  3.  
  4. anew task-FFT
  5.  
  6. \                                 FFT.JF
  7.  
  8. \ This text file contains a fixed-point Fast Fourier Transform
  9. \ program written in JForth.  To make the FFT words,
  10. \ two libraries were created containing complex math and vector
  11. \ words.
  12.  
  13. \ This code was published in the Dr. Dobb's Journal special
  14. \ FORTH issue in 1984.  
  15. \ Copyright (c) 1984 Joe Barnhart
  16.  
  17.  
  18.  
  19.  
  20. FORTH DEFINITIONS DECIMAL
  21.  
  22. : 4DROP 2DROP 2DROP ;
  23.  
  24. ( du< )
  25. ( for systems that don't have du< normally... )
  26.  
  27. FORTH DEFINITIONS DECIMAL
  28.  
  29. : FFT-DU< 
  30.   ROT SWAP 2DUP U<
  31.   IF   4DROP -1
  32.   ELSE = >R U< R> AND
  33.   THEN ;
  34.  
  35.  
  36. : SQR
  37.   DUP M* ;
  38.  
  39. ( complex data type and operations )
  40.  
  41. : COMPLEX@        2@ ;                ( addr  ---  x )
  42. : COMPLEX!        2! ;                ( x  addr --- )
  43.  
  44. : XVARIABLE             ( --- )
  45.   CREATE    CELL 2* HERE OVER  ( two elements per complex no.)
  46.             ERASE ALLOT         ( initialized to zero )
  47.   DOES> ;
  48.  
  49. : XCONSTANT             ( x --- )
  50.   CREATE    HERE CELL 2* 2DUP      ( two elements per cmplx )
  51.             ERASE ALLOT  COMPLEX!         ( save TOS into addr )
  52.   DOES>     COMPLEX@ ;                    ( return complex constant)
  53.  
  54. ( complex stack words )
  55.  
  56. : COMPLEXDUP      2DUP ;              ( x  ---  x  x )
  57. : XSWAP     2SWAP ;             ( x1  x2  ---  x2  x1 )
  58. : COMPLEXDROP     2DROP ;             ( x  --- )
  59. : XOVER     2OVER ;             ( x1 x2 --- x1 x2 x1 )
  60. : X2DUP     2OVER 2OVER ;       ( x1  x2  ---  x1  x2  x1  x2 )
  61.  
  62.  
  63.  
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. ( complex add, subtract, magnitude )
  71.  
  72. : X+        ROT + >R + R> ;             ( x1  x2  ---  x1+x2 )
  73. : X-        ROT SWAP - >R - R> ;        ( x1  x2  ---  x1-x2 )
  74. : X2/       2/ SWAP 2/ SWAP ;           ( x --- x/2 )
  75. : |X|^2     DUP * SWAP DUP * SWAP + ;   ( x --- s )
  76. : |X|       |X|^2  SQRT ;               ( x --- mag_x )
  77. : X'        SWAP NEGATE SWAP ;          ( x --- x_conjugate )
  78. : X0=       0= SWAP 0= AND ;            ( x --- t/f )
  79.  
  80.  
  81.  
  82.  
  83.  
  84.  
  85.  
  86. ( complex multiply )
  87.  
  88. \ This works.  Use it if the included assembler version gives you trouble.
  89.  
  90. \ : X*        >R  X2DUP           ( x1  x2  n ---  x1*x2/n )
  91. \             ROT R@ */                   ( real part = )
  92. \            ROT ROT R@ */  -            ( re1*re2 - im1*im2 )
  93. \            R> SWAP >R >R               ( save partial result)
  94. \            ROT ROT R@ */               ( re2*im1 )
  95. \            ROT ROT R> */  +            ( + re1*im2 )
  96. \            R> ;                        ( imag real --- )
  97.  
  98.  
  99.  
  100.  
  101.  
  102. \ Stack indexing is used to keep the original data available as
  103. \ the complex product is built.  When the product is complete, the
  104. \ original data is popped off the stack and the product pushed
  105. \ onto it.
  106. ( more complex operations )
  107.  
  108. : X*!       OVER >R >R          ( ^x1 ^x2 n --- x2 <-- x2*x1/n)
  109.             COMPLEX@  ROT COMPLEX@  R>              ( form x1*x2/n )
  110.             X*   R> COMPLEX! ;                ( save it into x2^ )
  111.  
  112. : X/        >R XSWAP XOVER      ( x1 x2  n  ---  n*x1/x2 )
  113.             X' R@ X*                    ( x1*x2' )
  114.             XSWAP |X| DUP R@ */         ( sqr|x2| )
  115.             SWAP OVER  R@ SWAP */       ( real part )
  116.             ROT ROT R> SWAP */ SWAP ;   ( imag part )
  117.  
  118.  
  119.  
  120.  
  121.  
  122. ( Vector definition )
  123.  
  124. : VECTOR        ( order scale word_size --- )
  125.   CREATE    ROT 2DUP , ,        ( store order, then wsize )
  126.             ROT ,               ( then store scale )
  127.             * DUP ,             ( bytes of entry storage )
  128.             HERE OVER ERASE     ( erase vector )
  129.             ALLOT               ( allot dictionary for entries)
  130.   DOES> ;                       ( return address of base )
  131.  
  132.  
  133.  
  134. \ This data structure holds a vector that contains "order"
  135. \ entries.  Legal indicies are in the range of 1 to "order".
  136. \ "Scale" is used in numeric operations on entries, and the
  137. \ "word_size" field describes the number of bytes per entry.
  138. ( vector operations )
  139.  
  140. : &ORDER        ( vector_descriptor --- order_of_vector )
  141.   @ ;
  142.  
  143. : &ESIZE        ( vec_desc --- size_of_entries )
  144.   CELL+ @ ;
  145.  
  146. : &SCALE        ( vec_desc --- scale_of_vec )
  147.  CELL 2* + @ ;
  148.  
  149. : &SIZE         ( vec_desc --- size_of_data_area )
  150.   CELL DUP 2* + + @ ;
  151.  
  152.  
  153.  
  154. ( matrix operations, cont. )
  155.  
  156. : [I]           ( index  vec_desc --- addr_of_entry )
  157.   DUP >R &ORDER                 ( save vec_desc, get entrysize)
  158.   OVER < OVER 1 < OR ABORT" Vector bounds exceeded... [I]"
  159.   1-  R@ &ESIZE *               ( index * entrysize )
  160.   CELL 4 * + R> + ;            ( +offset to entries )
  161.  
  162.  
  163.  
  164. \ This word calculates a storage address given a column number
  165. \ and a matrix descriptor.  If the desired address is outside
  166. \ of the alloted vector storage area, an exception is gener-
  167. \ ated.  This word is usually followed by "@" or "!".  The follow-
  168. \ ing screen contains an equivalent definition for 68000 Forth.
  169.  
  170.  
  171.  
  172. ( complex vector defs )
  173.  
  174. : XVECTOR       ( order scale --- )
  175.   CELL 2* VECTOR ;     ( declare a vector with two words/ent.)
  176.  
  177. : X@V           ( index vec_desc --- x )
  178.   [I] COMPLEX@ ;              ( fetch vector entry pointed by index)
  179.  
  180. : X!V           ( x index vec_desc --- )
  181.   [I] COMPLEX! ;              ( store into vector at position index)
  182.  
  183. \ print an xvector in table form
  184.  
  185. : XVECTOR.
  186.   CR DUP @ ." NUMBER OF SAMPLES:" SPACE . CR
  187.   CR ." REAL" CR
  188.   DUP>R DUP @ SWAP 16 + SWAP 0 DO
  189.       DUP I CELL 2* * + @ CELL 2* .R  
  190.       I 1+ CELL 2* MOD  NOT IF CR THEN
  191.   LOOP ." IMAGINARY" CR
  192.   DROP R> DUP @ SWAP 20 + SWAP 0 DO
  193.       DUP I CELL 2* * + @ CELL 2* .R  
  194.       I 1+ CELL 2* MOD  NOT IF CR THEN
  195.   LOOP
  196. DROP
  197. ;
  198.  
  199.  
  200. ( FFT constants and variables )
  201.  
  202. 10000   CONSTANT        10K             ( scaling constant )
  203.         VARIABLE        DIRECTION       ( direction of fft )
  204.         VARIABLE        VEC             ( current vector base)
  205.         VARIABLE        N               ( points in fft )
  206.         XVARIABLE       U               ( cmplx angle counter)
  207.         XVARIABLE       W               ( cmplx angle incr)
  208.  
  209.  
  210.  
  211.  
  212.  
  213.  
  214.  
  215.  
  216. ( bit reversal of vector )
  217.  
  218. : V[I][J]               ( i  j  --- [i]  [j] )
  219.   VEC @ [I] SWAP                ( calculate first addr )
  220.   VEC @ [I] SWAP ;              ( then second addr )
  221.  
  222. : XSWAPV                ( i  j  ---  )
  223.   V[I][J]                       ( calculate addresses )
  224.   2DUP SWAP >R >R               ( save one copy on ret stack)
  225.   COMPLEX@ ROT COMPLEX@  R> COMPLEX! R> COMPLEX! ;      ( swap cmplx numbers )
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232. ( bit reversal, cont. )
  233.  
  234. : BIT-REVERSE           ( --- )
  235.   N @ 1 SWAP 1 DO               ( for i=1 to size of vector )
  236.     DUP I > IF                  ( compare indicies )
  237.        DUP I XSWAPV             ( swap vector entries )
  238.     THEN
  239.     N @ 2/ SWAP                 ( let incr=points/2 )
  240.     BEGIN   2DUP <              ( while incr < indx )
  241.     WHILE   OVER - SWAP 2/ SWAP ( indx=indx-incr, incr=incr/2)
  242.     REPEAT  +                   ( indx=indx+incr )
  243.   LOOP DROP ;                   ( drop indx )
  244.  
  245.  
  246.  
  247.  
  248. ( butterfly calculation )
  249. : BFLY          ( j  i  ---  )
  250.   V[I][J] 2DUP >R >R            ( calc actual addresses )
  251.   COMPLEX@  U COMPLEX@  10K  X*      ( form temp product )
  252.   ROT COMPLEX@  XSWAP X2DUP           ( prepare stack for + and - )
  253.   X- >R >R  X+ R> R>            ( Fi-temp, Fi+temp )
  254.   DIRECTION @                         ( test for fwd/inv transform )
  255.   IF   XSWAP                    ( IFFT, no divide )
  256.   ELSE  X2/ XSWAP X2/           ( FFT, divide by POINTS???)
  257.   THEN R> COMPLEX!  R> COMPLEX!  ;          ( store products )
  258.  
  259.  
  260.  
  261.  
  262. ( 2** and log2 )
  263.  
  264. : 2**           ( n  ---  2**n )
  265.   0 MAX CELL 8 * 1- MIN        ( limit input data )
  266.   1 SWAP 0 DO 2* LOOP ;        ( DON'T loop if n=0 )
  267.  
  268. : LOG2          ( n  ---  log2[n] )
  269.   DUP 0= ABORT" Can't take log of zero. "
  270.   0 SWAP                        ( limit one word integer )
  271.   BEGIN   DUP 0>                ( test upper bit )
  272.   WHILE   2* SWAP 1+ SWAP       ( if zero, shift and incr count)
  273.   REPEAT  DROP                  ( done, drop remainder )
  274.   CELL 8 * 1- SWAP - ;         ( return log base 2 of input)
  275.  
  276.  
  277.  
  278. ( table of cosines and sines for 2**0 to 2**8 )
  279. CREATE WTAB
  280.          9 , CELL 2* , 10000 , 32 ,   ( entries, wsize, scale, data)
  281.         -10000 ,  00000 ,  00000 ,  -10000 ,  07071 ,  -07071 ,
  282.          09239 ,  -03827 ,  09808 , -01951 ,  09952 ,  -00980 ,
  283.          09988 ,  -00491 ,  09997 ,  -00245 , 09999 ,  -00123 ,
  284.  
  285.  
  286. \ This is a table of cosines and sines for angles of pi/n where n
  287. \ is an integral power of two.  The first entry is pi/1, cos is -1
  288. \ and sin is 0.  Next is pi/2, cos is 0 and sin is 1.  Entries up
  289. \ to pi/256 currently, but it can be extended for larger FFT's.
  290. \ Conjugates of the sine and cosine values are used in this table.
  291.  
  292.   
  293. ( inner-loop )
  294.  
  295. : INNER-LOOP    ( startpt increment --- )
  296.   DUP 2/                ( --s in in/2)   ( increment/2 is bfly offset)
  297.   ROT N @ 1+ SWAP       ( -- inc inc/2)    ( limits are startpt...N )
  298.   DO 
  299.    I 2DUP +        ( --inc inc/2 i inc/2+i) ( do a butterfly computation )
  300.         BFLY OVER       ( -- in in/2 in)  ( between I and I+offset )
  301.   +LOOP 2DROP ;         ( increment loop )
  302.  
  303. ( fft-kernel )
  304.  
  305. : FFT-KERNEL            ( --- ) ( N, VEC, must be set )
  306.   N @ LOG2 1+ 1                 ( do for i=1 to log2[N]-1 )
  307.   DO   
  308.         0 10K U COMPLEX!          ( init U to 1+j0 )
  309.         I WTAB X@V              ( get 1/-arg[pi/I])
  310.         DIRECTION @ IF X' THEN W COMPLEX!   ( !!!!take conjugate if IFFT )
  311.         I 2** DUP 2/ 1+ 1       ( increment and offset counter)
  312.         DO   I OVER INNER-LOOP  ( figure innermost loop )
  313.              W U 10K X*!        ( set new value for angle )
  314.         LOOP DROP               ( drop increment value )
  315.   LOOP ;
  316.  
  317.  
  318. ( fft and ifft )
  319.  
  320. : INIT          ( vector --- )
  321.   DUP VEC !
  322.   &ORDER N ! ;                  ( initialize size variable)
  323.  
  324. : DO-VECTOR     ( vector --- )
  325.   INIT                          ( initialize size field )
  326.   BIT-REVERSE                   ( bit swap input )
  327.   FFT-KERNEL ;                  ( perform butterflies )
  328.  
  329. : FFT         0 DIRECTION ! DO-VECTOR ;
  330. : IFFT       -1 DIRECTION ! DO-VECTOR ;
  331.  
  332.  
  333.  
  334.  
  335.  
  336. ( plotting words for complex vector )
  337.  
  338. VARIABLE MAG     256 MAG !     ( diff between max and min )
  339. VARIABLE MINI   -128 MINI !    ( min value in vector )
  340. VARIABLE PLOTBUF   68 ALLOT     ( buffer for building image )
  341. VARIABLE INTERVAL  1 INTERVAL ! ( allows skipping points )
  342. VARIABLE FUNC    ' NOOP FUNC !  ( controls RE, IM, or MAG plot)
  343.  
  344. : CHAR-MAP                   ( n --- n )
  345.   MINI @ -
  346.   70 MAG @ */ ;                 ( map into character coord. )
  347.  
  348.  
  349.  
  350.  
  351.  
  352. ( put char into buffer before printing line )
  353.  
  354. : BUF!                  ( char n --- )
  355.   69 MIN 0 MAX                  ( limit access to buffer )
  356.   PLOTBUF SWAP ROT FILL ;       ( store character into buffer)
  357.  
  358. : ?BUF!                 ( char n --- )
  359.   DUP 69 <= OVER 0 >= AND       ( check for out-of-bounds )
  360.   IF   PLOTBUF + C!             ( if in bounds, print char )
  361.   ELSE 2DROP
  362.   THEN ;
  363.  
  364.  
  365.  
  366.  
  367.  
  368. ( plot one line )
  369.  
  370. : PLOTDAT               ( n --- )
  371.   PLOTBUF 70 32 FILL            ( blank out plotter buffer )
  372.   ASCII : 0 CHAR-MAP ?BUF!           ( "zero" line if in range )
  373.   ASCII * SWAP CHAR-MAP BUF!         ( plot character for data )
  374.   PLOTBUF 70 -TRAILING TYPE ;   ( print the plot buffer )
  375.  
  376. : PLOTINDX              ( n --- )
  377.   5 .R                  ( print index number )
  378.   SPACE ASCII | EMIT ;       ( and right-hand bar )
  379.  
  380.  
  381.  
  382.  
  383.  
  384. ( plot entire vector )
  385.  
  386. : PLOT                  ( vec_desc --- )
  387.   CR  DUP &ORDER 1+ 1           ( do for each entry in vector)
  388.   DO  I PLOTINDX                ( get entry )
  389.       I OVER X@V                ( get entry )
  390.       FUNC @ EXECUTE            ( leave RE, IM, or MAG )
  391.       PLOTDAT CR                ( plot value )
  392.   INTERVAL @ +LOOP              ( loop by arbitrary increment)
  393.   CR DROP ;                     ( drop vec. descriptor )
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400. ( plot real, imag, or magnitude of vector )
  401.  
  402. : RE            SWAP DROP ;             ( leave real part )
  403. : IM            DROP ;                  ( leave imag part )
  404.  
  405. : PLOTMAG       ['] |X| FUNC ! PLOT ;   ( plot magnitude value)
  406. : PLOTRE        ['] RE  FUNC ! PLOT ;   ( plot real part )
  407. : PLOTIM        ['] IM  FUNC ! PLOT ;   ( plot imag part )
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416. ( load data from stack )
  417. ( TOS is vector descriptor, next is number of points, then )
  418. ( scale factor, data points follow )
  419. ( data-on-stack order #points name -- )
  420. : VECTOR!
  421.   OVER >R                    ( -- order #points name, save #points )
  422.   SWAP OVER !                ( -- order name,  #points assigned to vector)
  423.   SWAP OVER                  ( -- name order name )
  424.   CELL 2* + !                ( -- name, stores order at 8+ vector address)
  425.   1 R>  -DO                    ( for each vector entry, )
  426.      SWAP OVER               ( -- name data-point-from-stack name)
  427.      0 ROT                   ( -- name name 0 data-pt, set imag part to zero )
  428.      I                       ( -- name name 0 data-pt ith )
  429.      4 ROLL                  ( -- name 0 data-pt ith name )
  430.      X!V                     ( -- name 0, store into vector at ith position)
  431.   1 -LOOP DROP ;             ( work backwards thru vector )
  432.  
  433.  
  434. ( one cosine term, 32 data points )
  435.  
  436. 32 128 XVECTOR COSWAVE
  437. 0       128    0       0       0       0       0       0
  438. 0       0       0       0       0       0       0       0
  439. 0       0       0       0       0       0       0       0
  440. 0       0       0       0       0       0       0       0
  441. 128 32 COSWAVE VECTOR!
  442.  
  443.  
  444. ( After IFFT, should leave a cosine wave in the real array and)
  445. ( a sine wave in the imag array. )
  446.  
  447. ( five cosine terms, 32 data points )
  448.  
  449. 32 128 XVECTOR SQUARE
  450. 0       128    0       40     0       14     0       3
  451. 0       0       0       0       0       0       0       0
  452. 0       0       0       0       0       0       0       0
  453. 0       0       0       0       0       0       0       0
  454. 128 32 SQUARE VECTOR!
  455.  
  456.  
  457. ( After IFFT, should leave something like a square wave in the)
  458. ( imaginary array. )
  459.  
  460. \ add a constant across the board to a vector, used for testing and fun. 
  461. : ADDCONSTANT ( C ADR -- )
  462. DUP @ 1+ ( count)
  463. 1 DO
  464.   OVER DUP >R >R
  465.   I OVER X@V R> R> X+
  466.   I 3 PICK X!V
  467. LOOP
  468. 2DROP
  469. ;
  470.  
  471. \ 32 BENCH.FFT takes
  472. \ 14.84 seconds with X* coded in forth,
  473. \ 6.28 seconds with the assembler version of X* !!!
  474.  
  475. : BENCH.FFT
  476.  0 DO
  477.    SQUARE IFFT
  478.    SQUARE FFT
  479. LOOP
  480. ;
  481.