home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 202.lha / FFT / fft.hmsl < prev    next >
Encoding:
Text File  |  1988-12-28  |  10.9 KB  |  340 lines
  1. \                                FFT.HMSL
  2.  
  3. \ Nick Didkovsky, Robert Marsanyi                         August 21, 1988
  4.  
  5.  
  6.  
  7. \ These routines perform Fast Fourier Transform and inverse FFT
  8. \ on HMSL array or sample objects. Partials and phase angle data are made
  9. \ available in the shape editor for manipulation.
  10.  
  11. \ Usage is         sample-in fft.h
  12. \ or               sample-out ifft.h
  13.  
  14. \ This version runs end-to-end windows across the sample.
  15. \ Window size defaults to 32, changeable.
  16.  
  17. INCLUDE? FFT 4TH:FFT/FFT.JF 
  18. INCLUDE? ARCSIN 4TH:FFT/ARCSIN
  19. EXISTS? TERM.FFT.HMSL .IF TERM.FFT.HMSL .THEN
  20.  
  21. ANEW TASK-FFT_HMSL
  22.  
  23. OB.SAMPLE SAMPLE-IN
  24. OB.SAMPLE SAMPLE-OUT
  25. OB.BARRAY FFT-IN
  26. OB.BARRAY IFFT-OUT
  27. OB.SHAPE MAGNITUDES
  28. OB.SHAPE PHASE-INDICES
  29.  
  30. VARIABLE MAX-SAMPLE-SIZE ( 2048 data point dafault)
  31. VARIABLE WINDOW-SIZE     ( 32 default)
  32. VARIABLE SAMPLE-MAX      ( 128 for 8 bit sampling)
  33. VARIABLE DATA-WIDTH      ( 1 for 8 bit sampling)
  34. VARIABLE FFT-INITED?
  35. FFT-INITED? OFF
  36.  
  37. 256 128 XVECTOR WINDOW     ( initialize: 256 is alloted to WINDOW.)
  38.                            ( DON'T make window bigger than this or )
  39.                            ( you'll clobber your dictionary!!)
  40. : FREE.ALL.FFT
  41.   CLEAR: SHAPE-HOLDER
  42.   FREE: FFT-IN
  43.   FREE: IFFT-OUT
  44.   FREE: SAMPLE-OUT
  45.   FREE: MAGNITUDES
  46.   FREE: PHASE-INDICES
  47. ;
  48.  
  49. : SET.LIMITS.PHASE
  50.   WINDOW-SIZE @ 2/ 0 DO
  51.      -999 999 I PUT.DIM.LIMITS: PHASE-INDICES
  52.   LOOP
  53. ;
  54.  
  55.  
  56. : INIT.ALL.FFT.OBJECTS
  57.   FREE.ALL.FFT
  58.   MAX-SAMPLE-SIZE @ NEW: SAMPLE-OUT
  59.   MAX-SAMPLE-SIZE @ NEW: IFFT-OUT
  60.   MAX-SAMPLE-SIZE @ NEW: FFT-IN
  61.   MAX-SAMPLE-SIZE @ WINDOW-SIZE @ /  DUP DUP DUP ( #els for next 2 shapes)
  62.   WINDOW-SIZE @ 2/ 1+ NEW: MAGNITUDES     ( #partials = window-size/2+1 )   
  63.   SET.MANY: MAGNITUDES 
  64.   " HANDS OFF" WINDOW-SIZE @ 2/ PUT.DIM.NAME: MAGNITUDES
  65.   WINDOW-SIZE @ NEW: PHASE-INDICES ( 2 phase dims per partial, see below) 
  66. \  CLEAR: PHASE-INDICES
  67.   SET.MANY: PHASE-INDICES
  68.   SET.LIMITS.PHASE
  69.   MANY: SAMPLE-IN IF SAMPLE-IN ADD: SHAPE-HOLDER THEN ( anything slicker?)
  70.   SAMPLE-OUT ADD: SHAPE-HOLDER
  71.   PHASE-INDICES ADD: SHAPE-HOLDER
  72.   MAGNITUDES ADD: SHAPE-HOLDER
  73. ;
  74.  
  75. \ note that the WINDOWSIZE/2+1st dimension of MAGNITUDES contains the
  76. \ WINDOWSIZE/2+1st real value resulting from FFT
  77.  
  78. : SET.WINDOW.SIZE ( width -- )
  79.   DUP 256 > IF ABORT" Window width MUST be 256 or less!" THEN
  80.   DUP WINDOW-SIZE ! 
  81.   WINDOW !        ( -- , first word in WINDOW is number of entries)
  82.   INIT.ALL.FFT.OBJECTS
  83. ;
  84.  
  85. : SET.SAMPLE.SIZE ( size -- )
  86.   WINDOW-SIZE @ / WINDOW-SIZE @ * ( round off points to multiple of window)
  87.   MAX-SAMPLE-SIZE !
  88.   INIT.ALL.FFT.OBJECTS
  89. ;
  90.  
  91. : REMIND.ME 
  92.   cr cr
  93.   ." SET.WINDOW.SIZE ( size -- )  will change WINDOW width." cr
  94.   ." SET.SAMPLE.SIZE ( size -- )  changes maximum length of sample" cr
  95.   ." <sample> FFT.H runs the Fast Fourier Transform, " CR
  96.   ." <target> IFFT.H is the inverse." CR CR
  97.   ." REMIND.ME to see this message again" cr 
  98. ;
  99.  
  100. : INIT.FFT.HMSL
  101.   2048 MAX-SAMPLE-SIZE !
  102.   128 SAMPLE-MAX !
  103.   1 DATA-WIDTH !
  104.   cr cr
  105.   ." SAMPLE-MAX defaulting to 128 (8 bit sampling)" cr 
  106.   ." DATA-WIDTH defaulting to 1   (8 bit sampling)" cr
  107.   ." DON'T CHANGE ONE WITHOUT THE OTHER!!!" cr cr
  108.   ." WINDOW-SIZE defaulting to 32" cr 
  109.   ." MAX-SAMPLE-SIZE defaulting to 2048" cr
  110.   REMIND.ME 
  111.   32 SET.WINDOW.SIZE
  112.   MAX-SAMPLE-SIZE @ NEW: SAMPLE-OUT
  113.   MAX-SAMPLE-SIZE @ NEW: IFFT-OUT
  114.   MAX-SAMPLE-SIZE @ NEW: FFT-IN
  115.   LOADTABLE                  ( sine lookup table)
  116.   FFT-INITED? ON
  117. ;
  118.  
  119. : TERM.FFT.HMSL
  120.   FREE: SAMPLE-IN
  121.   FREE.ALL.FFT
  122. ;
  123.  
  124.  
  125.  
  126.  
  127. \ **************************** FFT code starts ****************************
  128.  
  129.  
  130. \ WINDOW contains WINDOW-SIZE spectrum entries. We only want the first half,
  131. \ since the second half is just complex conjugate of first half.
  132.  
  133.  
  134.  
  135.  
  136. : ROUND ( tens-int -- ones-int)
  137.   DUP 0< IF NEGATE  5 + 10 / NEGATE ELSE 5 + 10 / THEN
  138. ;
  139.  
  140.  
  141. : |X|.SCALE.UP 
  142.   |X|^2 10000 * SQRT ;  
  143.  
  144. : LOAD.WINDOW ( el# -- , load 32 points into WINDOW starting at el#)
  145.   DUP WINDOW-SIZE @ + SWAP DO 
  146.       I AT: FFT-IN b->s ( load up stack with signed audio data)
  147.   LOOP
  148.   SAMPLE-MAX @ WINDOW-SIZE @ WINDOW VECTOR! ( interleaved load to WINDOW)
  149. ;
  150.  
  151. \ These words load up one time slice. ie, all elements in one dimension.
  152. : LOAD.MAGNITUDES ( el# --, load magnitudes of partials into MAGNITUDES)
  153.   WINDOW-SIZE @ 2/ 0 DO
  154.      I 1+ WINDOW X@V |X| ( complex fetch, calc magnitude)
  155.      OVER I ED.TO: MAGNITUDES ( el# -- store |X| into MAGNITUDES)
  156.   LOOP
  157. \ Now for that weird window-size/2+1st value !
  158. \ Note that (for example) the 17th element of WINDOW goes to 
  159. \ the 16th dimension of MAGNITUDES
  160.   WINDOW-SIZE @ 2/ DUP 1+ ( -- el# window-index window-index+1)
  161.   WINDOW X@V NIP          ( -- el# window-index real-value)
  162.   -ROT ED.TO: MAGNITUDES
  163. ;
  164.  
  165. \ PHASE-INDICES has 2 dimensions for each partial.  Partial 0 maps to
  166. \ dimension 0 and 1, where 0 holds the sine-table index, and 1 holds
  167. \ a flag for the real component's sign.  Since we are using a sine table,
  168. \ we calculate phase angle from imaginary component and magnitude.  This
  169. \ way, the sign of the real component can get lost, unless it is flagged
  170. \ some other way.  Here's the way.
  171.  
  172. VARIABLE SIGN-FLAG
  173.  
  174. : LOAD.PHASE-INDICES ( el# -- )
  175.   0 SIGN-FLAG !
  176.   WINDOW-SIZE @ 2/ 0 DO                  ( loop through partials)
  177.      DUP I ED.AT: MAGNITUDES ( -- el# |X|, borrow magnitude from MAGNITUDES)
  178.      DUP 0= IF DROP 0 OVER I 2* ED.TO: PHASE-INDICES ( no mag? no phase!)
  179.                     0 OVER I 2* 1+ ED.TO: PHASE-INDICES ( -- el#, zero flag)    
  180.        ELSE 
  181.          I 1+ WINDOW X@V          ( -- el# |X| imag real)
  182.          0< SIGN-FLAG !     ( -- el# |X| imag, save sign of real component)
  183.                                  ( -- el# |X| imag, only want imaginary part)
  184.          10000 ROT               ( -- el# imag 10000 |X|)
  185.          */                      ( scaled divide)
  186.          ARCSIN.INDEX            ( -- el# arcsinindex)
  187.          OVER I 2* ED.TO: PHASE-INDICES  ( -- el# )
  188.          SIGN-FLAG @ OVER I 2* 1+ ED.TO: PHASE-INDICES ( -- el#)
  189.        THEN
  190.    LOOP
  191.    DROP ( don't need el#)
  192. ;
  193.  
  194. \ This word sets data-width (1 for 8 bit sampling)
  195. \ Don't change data-width without changing SAMPLE-MAX !!!
  196. \ Note that the data crammed into FFT-IN will be unsigned.
  197. \ It becomes signed in LOAD.WINDOW 
  198.  
  199. : COPY.TO.FFT.IN ( sample_addr -- )
  200.   DUP DUP        ( -- addr addr addr)
  201.   WIDTH: [] DATA-WIDTH !  ( -- addr addr) 
  202.   FREE: FFT-IN
  203.   DATA-WIDTH @ SET.WIDTH: FFT-IN 
  204.   DATA-WIDTH @ MAX-SAMPLE-SIZE @ * NEW: FFT-IN 
  205. \  CLEAR: FFT-IN
  206.   MANY: []   ( -- addr many-in-input)
  207.   MAX-SAMPLE-SIZE @ MIN ( -- addr min-number-of-data-to-move)
  208.   DATA-WIDTH @ *        ( -- addr min-number-scaled)
  209.   SWAP
  210.   data.addr: []  data.addr: fft-in rot move
  211. ;
  212.  
  213. \ ******************** THIS IS THE MAIN FFT WINDOW LOOP *******************
  214.  
  215. \ see if the imaginary components are indeed always 0 in 1st and 17th
  216.  
  217. : FUNKY  WINDOW XVECTOR. ; ( print WINDOW)
  218.  
  219. \ The window scans through sample end-to-end.
  220. \ Usage is:   sample-name fft.h
  221.  
  222. : FFT.H ( sample_addr -- )
  223.   cr
  224.   FFT-INITED? @ NOT IF INIT.FFT.HMSL THEN
  225.   DEPTH IF COPY.TO.FFT.IN ELSE ." Need a sample name!" ABORT THEN
  226.   MAX-SAMPLE-SIZE @ WINDOW-SIZE @ / 0 DO   ( # windows loop)
  227.       I WINDOW-SIZE @ *  ( starting point)
  228.       LOAD.WINDOW
  229.       WINDOW FFT   ( fft done in place! WINDOW contains spectrum now)
  230.       I LOAD.MAGNITUDES
  231.       I LOAD.PHASE-INDICES
  232.       ." ." FLUSHEMIT
  233.   LOOP
  234.   CR
  235. ;
  236. \ ************************** IFFT code starts **************************
  237.  
  238. \ real and imaginary components are retrieved using phase angle and
  239. \ magnitude.  The code grabs 1 magnitudes and 1 phase angle
  240. \ from MAGNITUDES and PHASE-INDICES respectively,
  241. \ and derives real and complex components for FIRST half of WINDOW.
  242. \ The SECOND half of WINDOW consists of complex conjugates of the 1st half.
  243. \ remember window indexes from 1..WINDOW-SIZE, NOT 0..WINDOW-SIZE-1
  244. \ This routine also takes advantage of the odd/even of the el#, inverting
  245. \ the data on odd values, leaving even values alone.
  246.  
  247. VARIABLE SIGN-BACK
  248.  
  249. : INV.WINDOW.1 ( el# --)
  250.   0 SIGN-BACK !
  251.   WINDOW-SIZE @ 2/ 0 DO ( loop through dimensions)
  252.        DUP I 2DUP  ( -- el# el# i el# i)
  253.        ED.AT: MAGNITUDES ( -- el# el# i |X|)             
  254.        DUP 0= IF NIP NIP 0  ( -- el# 0 0, don't bother multiplying by 0)
  255.        ELSE
  256.          -ROT                 ( -- el# |X| el# i)
  257.          2DUP                 ( -- el# |X| el# i el# i)
  258.          2* ED.AT: PHASE-INDICES ( -- el# |X| el# i index-to-table)
  259.          -ROT                    ( -- el# |X| index el# i)
  260.          2* 1+ ED.AT: PHASE-INDICES SIGN-BACK ! ( set flag for cosine)
  261.          2DUP                 ( -- el# |X| index |X| index)
  262.          SIN.LOOKUP           ( -- el# |X| ind |X| sin) 
  263.          1000  */  ROUND      ( -- el# |X| ind imaginary-component)
  264.          -ROT
  265.           COS.LOOKUP
  266.          1000 */  ROUND       ( -- el# imag real)
  267.          SIGN-BACK @ IF NEGATE THEN    ( -- el# imag sign-corrected-real)
  268.        THEN
  269.        I 1+ WINDOW X!V      ( -- el#, store result to window)
  270.   LOOP
  271. \ The next two lines retrieve the WINDOW-SIZE/2+1st real value.
  272. \ Same notice: there is an index difference between where MAGNITUDES
  273. \ stores the value and index where it goes in WINDOW.
  274. \ If the window size is 32, the 17th entry into WINDOW is stored in the
  275. \ 16the dimension of MAGNITUDES.
  276.   WINDOW-SIZE @ 2/ DUP 1+ >R ED.AT: MAGNITUDES 0 SWAP 
  277.   R> WINDOW X!V  
  278. ; 
  279.  
  280. \ The remaining complex numbers are conjugates of the first half
  281. \ example, for 32 wide window: 2->32, 3->31, 4->30, 5->29, etc etc 
  282. \ 17th value was retrieved from 17th dimension of MAGNITUDES. 
  283. : INV.WINDOW.2 
  284.   WINDOW-SIZE @ 2/ 1+ 2 DO
  285.       I WINDOW X@V X'   ( take conjugate)
  286.      WINDOW-SIZE @ 2+ I - WINDOW X!V ( store symmetrically)
  287.   LOOP
  288. ;
  289.  
  290. \ Take the result of one window IFFT and load up a chunk of IFFT-OUT
  291. : STUFF.FROM.WINDOW  ( ifft# -- )
  292.   WINDOW-SIZE @ * ( starting index in IFFT-OUT )
  293.   WINDOW-SIZE @ 1+ 1 DO
  294.      I WINDOW X@V NIP OVER I + 1- TO: IFFT-OUT ( store real part)
  295.   LOOP
  296.   DROP
  297. ;
  298.  
  299. \ This word moves the result of IFFT from array IFFT-OUT to target object.
  300. \ Most likely this will be a sample object, but it could be an array or ... 
  301. \ just in case it's SAMPLE-OUT you're using, I save you the trouble
  302. \ doing the put.#oneshot: and set.many: on it.
  303.  
  304. : MOVE.TO.TARGET ( target.addr -- )
  305.   DUP            ( -- addr addr )
  306.   FREE: []
  307.   DATA-WIDTH @ OVER SET.WIDTH: []
  308.   MAX-SAMPLE-SIZE @ OVER NEW: []
  309.   DATA.ADDR: [] DATA.ADDR: IFFT-OUT SWAP MAX-SAMPLE-SIZE @ DATA-WIDTH @ * MOVE
  310.   MAX-SAMPLE-SIZE @ SET.MANY: SAMPLE-OUT
  311.   MAX-SAMPLE-SIZE @ PUT.#ONESHOT: SAMPLE-OUT 
  312. ;
  313.  
  314. : IFFT.H ( target.object -- )
  315.   cr FREE: IFFT-OUT
  316.   DATA-WIDTH @ SET.WIDTH: IFFT-OUT 
  317.   MAX-SAMPLE-SIZE @ NEW: IFFT-OUT
  318. \  CLEAR: IFFT-OUT
  319.   MAX-SAMPLE-SIZE @ WINDOW-SIZE @ /  0 DO
  320.   I INV.WINDOW.1  ( el# -- )
  321.       INV.WINDOW.2
  322.       WINDOW IFFT ( ifft done in place! WINDOW contains data points now)
  323.       I STUFF.FROM.WINDOW 
  324.       ." ." FLUSHEMIT
  325.   LOOP
  326.   MOVE.TO.TARGET ( target.object -- )
  327. ;
  328.  
  329.  
  330.      
  331. CR CR  
  332.    ." ***************************** FFT ************************************"
  333. CR CR
  334.    ."   Be sure to LOAD: a sample into SAMPLE-IN before you INIT.FFT.HMSL!" 
  335. CR CR
  336.    ." **********************************************************************"
  337. CR CR
  338.  
  339.  
  340.