home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 326.lha / KFFT_v1.1 / fft1 < prev    next >
Text File  |  1989-12-23  |  6KB  |  205 lines
  1. \ FFT1 - Fast Fourier Transform JForth* source code.
  2. \
  3. \ KFFT V1.1 (C)Copyright 1989, Jerry Kallaus.  All rights reserved. 
  4. \ May be freely redistributed for non-commercial use (FREEWARE).
  5. \
  6. \ Conditionally compiles as various floating point and fixed point
  7. \ versions of a Cooley-Tukey radix-2 FFT algorithm.
  8. \ Interface/entry words FFT and IFFT are at the end of this file.
  9. \ File FFTinc contains INCLUDEs for all needed files.
  10. \ File FFTControls contains compilation controls and parameters.
  11. \ File Cmplx contains complex arithmetic and data definition words.
  12. \ File FFT.asm contains support assembly source code and
  13. \ file FFT.cod contains the corresponding 'object' code.
  14. \ Words containing curly braces after the word name use JForth's
  15. \ local variables capability.
  16. \ See accompaning file FFT.doc for details.
  17. \
  18. \ * JForth is a product of Delta Research.
  19.  
  20.  
  21. ANEW task-fft1
  22.  
  23. VARIABLE reversal-fft         VARIABLE inv-fft
  24. reversal-fft  OFF
  25.  
  26. \ Variables specific to fixed point versions.
  27.  
  28. float_fft? NOT .IF
  29.         VARIABLE inshift-fft   VARIABLE shifts-fft   VARIABLE outbits-fft
  30.         inshift-fft   OFF      shifts-fft OFF        16 outbits-fft !
  31.         VARIABLE evens-fft     VARIABLE odds-fft     VARIABLE blk-exp-fft
  32.         evens-fft     OFF      odds-fft   ON
  33. .THEN
  34.  
  35.  
  36. \ Perform sanity check on fft arguments.
  37.  
  38. : CHECK.INPUTS.FFT { a m  --- }  ( a is data addr, m is log2-fft-size )
  39.   a 1 AND
  40.          IF  CR ." FFT odd array address ????" ABORT THEN
  41.   m 2 max_log2_fft  WITHIN? NOT
  42.   IF  CR ." FFT size is 2**" m . ."  = "  m 2**  .
  43.        ." ????"
  44.       CR ." Set max_log2_fft for sizes greater than 2**" max_log2_fft .
  45.       ABORT
  46.   THEN
  47. ;
  48.  
  49. \ Either perform bit reversal of data; or make reversal map if the
  50. \ variable reversal-fft contains address to hold reversal map.
  51.  
  52. : BIT.REVERSAL  { a m | jj n/2 ai aj rend --- }
  53.                                   ( a is data addr, m is log2-fft-size)
  54.    a m  CHECK.INPUTS.FFT
  55.    reversal-fft @  -> rend
  56.    0 -> jj
  57.    m 2**  CELLS DUP -> n/2
  58.    CELL- 2* 0
  59.    DO
  60.      jj I > IF
  61.               rend IF
  62.                      jj I rend Z!                ( make reversal map )
  63.                      2CELL rend +  -> rend
  64.                      0 0 rend Z!
  65.                    ELSE                          ( else scramble data )
  66.                      jj I a DUP  D+
  67.                      -> ai  DUP  -> aj
  68.                         Z@  ai Z@
  69.                      aj Z!  ai Z!
  70.                    THEN
  71.             THEN
  72.       n/2
  73.       BEGIN DUP DUP jj <=             ( k  = jj   )
  74.       WHILE
  75.            jj - NEGATE  -> jj         ( jj = jj-k )
  76.            U2/                        ( k  = k/2  )
  77.       REPEAT
  78.       DROP jj +  -> jj                ( jj = j+k  )
  79.    2CELL +LOOP
  80. ;
  81.  
  82. asm_fft? NOT .IF
  83.  
  84. : QUICK.REVERSAL { a r | ai aj --- }  ( data-addr  reversal-map-addr -- )
  85.     r 2CELL-    -> r
  86.     BEGIN r 2CELL+ DUP  -> r  Z@ DUP
  87.     WHILE
  88.          a DUP  D+
  89.          -> ai  DUP  -> aj
  90.             Z@  ai Z@
  91.          aj Z!  ai Z!
  92.     REPEAT  2DROP
  93. ;
  94. .THEN
  95.  
  96. \ ************************  Begin  FFT1  ****************************
  97. \
  98. \ Cooley-Tukey FFT algorithm
  99.  
  100. : FFT1 { a m | n le le1 ss hi 2 u 2 w --- }
  101. \  ( cmplx-array-address  m -- )  ( m is log base 2 of fft size )
  102. \  Byte address indexing used for efficiency.
  103.  
  104.   a m  CHECK.INPUTS.FFT
  105.   reversal-fft @
  106.     IF    a reversal-fft @ QUICK.REVERSAL
  107.     ELSE  a m  BIT.REVERSAL
  108.     THEN
  109.     m 2**  ZCELLS  -> n
  110.     0  -> hi
  111. [ auto_scale_fft? ]  .IF  shifts-fft OFF  .THEN
  112.  
  113. \ ***********************  Begin Outer Loop  *************************
  114.   m 0
  115.   DO
  116.  
  117. [ auto_scale_fft? ] .IF
  118.     I  IF                                ( auto scale-down control:       )
  119.       [ 14 2** ] LITERAL                 ( Shift input to new stage right )
  120.       hi <   IF 1 ELSE 0 THEN            ( 0,1 or 2 places depending on   )
  121.       [ 15 2** ] LITERAL                 ( highest bit set on output of   )
  122.       hi <   IF 1+       THEN            ( last stage.  if > 2**14 use 1  )
  123.       DUP  -> ss   shifts-fft +!         (              if > 2**15 use 2  )
  124.     ELSE   inshift-fft @       -> ss     ( if stage 0, use inshift-fft    )
  125.     THEN
  126.     0  -> hi
  127. .THEN
  128.  
  129. [ auto_scale_fft? float_fft? OR NOT ]  .IF
  130.          I      1 AND odds-fft  @ AND               ( manual even/odd )
  131.          I COMP 1 AND evens-fft @ AND OR  -> ss     ( stage control )
  132.          I IF-NOT inshift-fft @ ss +      -> ss
  133.                   inshift-fft OFF
  134.            THEN   
  135. .THEN
  136.  
  137.     I 2**      DUP  -> le1
  138.     2* ZCELLS  -> le
  139. [ float_fft? w_table_fft? NOT AND ]
  140. .IF
  141.     le1  FLOAT  pi_fft SWAP F/
  142.     ZEXP                           ( w=[cos,sin] )
  143. .ELSE
  144.     I ZCELLS  w-table-fft + Z@     ( w=[cos,sin] )
  145. .THEN
  146.     inv-fft @ IF-NOT  CONJG  THEN  ( reverse oscilator for forward fft )
  147.         -> w
  148.     Z1  -> u
  149. \ *************************  Begin Middle Loop  ***************
  150.  
  151.     le1 ZCELLS DUP -> le1  0
  152.     DO
  153. \ *************************  Begin Inner Loop  ****************
  154.  
  155. [ inner_asm_fft? ]
  156. .IF
  157.       u le ss n a I le1   INNER.LOOP
  158. .THEN
  159.  
  160. [ auto_scale_fft? ]                         .IF  hi OR  -> hi  .THEN
  161. [ auto_scale_fft? NOT inner_asm_fft? AND ]  .IF  drop          .THEN
  162.  
  163. [ inner_asm_fft? NOT ]  .IF
  164.     n I
  165.       DO
  166.         a I + DUP>R
  167.         Z@ ss ZSCALE.DOWN
  168.         R@ le1 + DUP>R
  169.         Z@ ss ZSCALE.DOWN u Z*
  170.         Z2DUP
  171.         Z- R> Z!
  172.         Z+ R> Z!
  173.         le +LOOP
  174. \ ************************  End of Inner Loop  **********************
  175. .THEN
  176.     w u Z*  -> u
  177.     ZCELL +LOOP
  178.   LOOP
  179.  
  180. \ *************************  End of Middle Loop  **********************
  181. \ *************************  End of Outer  Loop  **********************
  182.  
  183. [ auto_scale_fft? ]  .IF     ( Normalize data to outbits-fft )
  184.   m 2**  2*  -> n            ( Set n to number of cells )
  185.   hi                         ( hi=0 => data all zero => nothing to do )
  186.   IF
  187.     a  n
  188.     outbits-fft @ hi NSBITS  -
  189.     DUP shifts-fft +!  ASHIFT.ARRAY
  190.   THEN
  191.   shifts-fft @  blk-exp-fft +!           ( update block-flt-pt exponent )
  192. .THEN
  193. ;
  194.  
  195. \ *************************  End of FFT1  *****************************
  196.  
  197. \ m is log base 2 of fft-size
  198.  
  199. : INIT.MAP.FFT  { map-addr m --- }
  200.      map-addr  reversal-fft !  0 m BIT.REVERSAL ;
  201.  
  202. : FFT  ( cmplx-array-addr m -- )  inv-fft OFF FFT1 ;  ( forward fft )
  203.  
  204. : IFFT ( cmplx-array-addr m -- )  inv-fft ON  FFT1 ;  ( inverse fft )
  205.