home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 326.lha / KFFT_v1.1 / Demo_FFT1 < prev    next >
Text File  |  1989-12-23  |  7KB  |  278 lines
  1. \ Simplified demonstration, by Phil Burk, of the use of KFFT.
  2. \
  3. \ The KFFT package is very powerful, and has many options and
  4. \ variations for its use.  For this reason, it can be a little
  5. \ difficult to figure out.  I have written this simplified demo
  6. \ for people, like me, who are not already experts in Digital
  7. \ Signal Processing.  Once you understand this demo, please
  8. \ refer to the excellent documents that Jerry has provided.
  9. \
  10. \ The FFT is a Fast way to do a Fourier Transform.  A Fourier
  11. \ Transform is an algorithm that converts a signal to its
  12. \ component frequencies.  Imagine a graphic spectrum display
  13. \ on a music stereo. When the bass plays, the display indicates
  14. \ activity at low frequencies. When the flute plays, or a cymbal
  15. \ crashes, we see things light up at the high end of the spectrum.
  16. \ An FFT could be used to generate a spectrum from an input
  17. \ musical signal.  Unfortunately, the FFT requires a lot of
  18. \ computation which means it is difficult to make it work
  19. \ for real time music applications.  It is typically used for
  20. \ analysing signals that are slowly changing or for "after
  21. \ the fact" analysis of sound. Jerry has done a great
  22. \ job, however, at making this FFT particularly fast so I
  23. \ am excited about seeing what it can do.
  24. \
  25. \ FFTs are also used for analysing medical signals, ECGs and EEGs,
  26. \ analysing seismic records, stock market data, or almost
  27. \ any data that has periodic components. It is an extremely
  28. \ useful tool and I am grateful to Jerry for making this available.
  29. \
  30. \ This demo will only use the Fixed Point, integer, version
  31. \ of the FFT.  It will not use Floating Point.  Fixed Point
  32. \ is faster and is often more useful for audio applications
  33. \ which is my main interest.
  34. \
  35. \ The main word for doing the FFT is FFTRC which has the following
  36. \ stack diagram:
  37. \
  38. \     FFTRC  ( data-address num_points_log_base2 -- )
  39. \
  40. \     The data-address is the address of the data array.
  41. \         Each point is a 4 byte, 32 bit value.
  42. \
  43. \     The num_points_log_base2 is the log base 2 of the
  44. \         size of the array. The size must be some power
  45. \         of 2, eg. 8=2**3,   or 1024=2**10.
  46. \
  47. \ Thus for a 1024 point array at address MY-DATA, use:
  48. \
  49. \         CREATE MY-DATA 1024 4* ALLOT
  50. \         MY-DATA 10 FFTRC
  51. \
  52. \ The output of the FFT will be left in the data array.
  53. \
  54. \ Author: Phil Burk
  55. \ Copyright 1989 Phil Burk - (FreeWare)
  56.  
  57. decimal
  58.  
  59. \ UnAssign sp: so it points to current directory.
  60. assign sp: " "
  61.  
  62. \ load required FFT code
  63. INCLUDE?  fft    fftinc
  64. INCLUDE?  fftrc  fftrc
  65. INCLUDE?  ifftcr ifftcr
  66.  
  67. \ include random number generator just for this demo.
  68. include? choose ju:random
  69.  
  70. \ The choice between floating point and fixed point is set
  71. \ in the file "fftcontrols".  The default is fixed point.
  72. float_fft? 
  73. .IF ." Woops! Floating Point version! Edit fftcontrols file." abort
  74. .THEN
  75.  
  76. ANEW TASK-DEMO_FFT1
  77. decimal
  78. \ This is the log base 2 of the size of the data array.
  79. \ In other words, 2**4=16, or 2**5=32, which is the actual size.
  80. 4 constant LOG2_FFT_SIZE
  81.  
  82. \ Allocate array for the data. Each point is 4 bytes (32 bits)
  83. log2_fft_size 2** constant DATA_SIZE
  84. data_size array SIGNAL-DATA
  85. \ Size is data_size longwords, or data_size*4 bytes.
  86.  
  87. \ This integer value will correspond to 1.0 or unity.
  88. $ 1000 constant VIRTUAL_ONE
  89. \ Thus $ 800 would be 0.5
  90.  
  91. : CLEAR.SIGNAL  ( -- , clear signal array )
  92.     data_size 0
  93.     DO 0 i signal-data !
  94.     LOOP
  95. ;
  96.  
  97. : ADD.SAWTOOTHS  ( -- , add 2 sawtooths to signal)
  98.     0 ( starting index )
  99.     2 0
  100.     DO  data_size 2/ 0
  101.         DO virtual_one data_size 2/ 2/ / i * negate
  102.            virtual_one +  ( -- index value )
  103.            over
  104.            signal-data +!
  105.            1+ ( incr index )
  106.         LOOP
  107.     LOOP drop
  108. ;
  109.  
  110. : ADD.PULSE ( width -- , add pulse to data)
  111.     dup 0
  112.     DO virtual_one i signal-data +!
  113.     LOOP
  114.     data_size swap
  115.     DO 0 i signal-data +!
  116.     LOOP
  117. ;
  118.  
  119. : ADD.RANDOM ( -- , add random noise )
  120.     data_size 0
  121.     DO virtual_one 2* choose virtual_one -
  122.        i signal-data +!
  123.     LOOP
  124. ;
  125.  
  126. : PLOT.VALUE ( value min max -- , plot a value using '*' characters)
  127.     over - >r ( -- value min, -r- range )
  128.     - 0 max
  129.     40 r> */  ( normalize to 0-40 chars )
  130.     40 min 0
  131.     DO ascii * emit
  132.     LOOP
  133. ;
  134.  
  135. : RANGE.DATA ( address #cells -- min max , calculate range )
  136.     >r >r
  137.     -1 -1 shift ( maximum integer )
  138.     dup negate
  139.     r> r> 0
  140.     DO  ( min max addr )
  141.         dup @ >r  ( get value )
  142.         cell+ -rot  ( addr min max )
  143.         r@ max -rot
  144.         r> min -rot
  145.     LOOP drop
  146. ;
  147.  
  148. : PLOT.DATA ( -- , display character plot of data )
  149.     0 signal-data data_size range.data
  150.     data_size 0
  151.     DO  cr i 3 .r
  152.         i signal-data @  10 .r 2 spaces
  153. \
  154. \ print appropriate number of '*'s to plot data
  155.         i signal-data @  ( min max val )
  156.         2 pick 2 pick plot.value
  157.     LOOP 2drop cr
  158. ;
  159.  
  160. : RUN.FFT  ( -- , calculate FFT for input signal )
  161. \ We are not using a precalculated bit reversal map so...
  162.     reversal-fft off
  163. \
  164.     ." Input signal:" cr
  165.     plot.data ." Hit RETURN to continue:" key drop cr
  166. \
  167.     ." Calculating FFT" cr
  168.     0 signal-data  ( address of first point )
  169.     log2_fft_size  ( 2's exponent )
  170.     FFTRC
  171. \
  172.     ." Fourier Transform of Signal:" cr
  173.     plot.data ." Hit RETURN to continue:" key drop cr
  174. \
  175.     ." Calculating Inverse FFT" cr
  176.     0 signal-data  ( address of first point )
  177.     log2_fft_size
  178.     IFFTCR
  179. \
  180.     ." Reconstructed version of Signal:" cr
  181.     plot.data
  182. ;
  183.  
  184. : SAW.FFT ( -- , run FFT on sawtooth )
  185.     clear.signal
  186.     add.sawtooths
  187.     run.fft
  188. ;
  189.  
  190.  
  191. : PULSE.FFT ( -- , run FFT on sawtooth )
  192.     clear.signal
  193.     7 add.pulse
  194.     run.fft
  195. ;
  196.  
  197. : RANDOM.FFT ( -- , run FFT on random signal )
  198.     clear.signal
  199.     add.random
  200.     run.fft
  201. ;
  202.  
  203. : NOISY.FFT ( -- , run FFT on sawtooth plus noise )
  204.     clear.signal
  205.     add.sawtooths
  206.     add.random
  207.     run.fft
  208. ;
  209.  
  210.  
  211. \ By modifying the spectrum, then doing an inverse
  212. \ transform, one can simulate the effect of a complex filter.
  213.  
  214. : RUN.FFT.FILTER  ( -- , use the FFT like a filter )
  215. \ We are not using a precalculated bit reversal map so...
  216.     reversal-fft off
  217. \
  218.     ." Input signal:" cr
  219.     plot.data ." Hit RETURN to continue:" key drop cr
  220. \
  221.     ." Calculating FFT" cr
  222.     0 signal-data  ( address of first point )
  223.     log2_fft_size
  224.     fftrc
  225. \
  226. \ Apply low pass filter to spectrum.
  227.     data_size 4 / 0
  228.     DO data_size i - 1- signal-data off
  229.     LOOP
  230. \
  231.     ." Fourier Transform of Signal with high frequencies removed:" cr
  232.     plot.data ." Hit RETURN to continue:" key drop cr
  233. \
  234.     ." Calculating Inverse FFT" cr
  235.     0 signal-data  ( address of first point )
  236.     log2_fft_size
  237.     ifftcr
  238. \
  239.     ." Reconstructed version of Signal after Filtering:" cr
  240.     plot.data
  241. ;
  242.  
  243. : SAW.FILTER ( -- , run FFT filter on sawtooth )
  244.     clear.signal
  245.     add.sawtooths
  246.     run.fft.filter
  247. ;
  248.  
  249. : PULSE.FILTER ( -- , run FFT filter on pulse )
  250.     clear.signal
  251.     5 add.pulse
  252.     run.fft.filter
  253. ;
  254.  
  255. : RANDOM.FILTER ( -- , run FFT on random signal )
  256.     clear.signal
  257.     add.random
  258.     run.fft.filter
  259. ;
  260.  
  261. : NOISY.FILTER ( -- , run FFT on sawtooth plus noise )
  262.     clear.signal
  263.     add.sawtooths
  264.     add.random
  265.     run.fft.filter
  266. ;
  267.  
  268. cr cr
  269. ." Enter:  SAW.FFT  to see a sawtooth transformed, then" cr
  270. ." reconstructed.  You can also enter:" cr
  271. ."     PULSE.FFT" cr
  272. ."     RANDOM.FFT" cr
  273. ."     NOISY.FFT" cr
  274. cr
  275. ." To see how an FFT can be used as a filter, enter:  SAW.FILTER" cr
  276. ." Read this file for more information." cr
  277.  
  278.