\ Simplified demonstration, by Phil Burk, of the use of KFFT.
\
\ The KFFT package is very powerful, and has many options and
\ variations for its use.  For this reason, it can be a little
\ difficult to figure out.  I have written this simplified demo
\ for people, like me, who are not already experts in Digital
\ Signal Processing.  Once you understand this demo, please
\ refer to the excellent documents that Jerry has provided.
\
\ The FFT is a Fast way to do a Fourier Transform.  A Fourier
\ Transform is an algorithm that converts a signal to its
\ component frequencies.  Imagine a graphic spectrum display
\ on a music stereo. When the bass plays, the display indicates
\ activity at low frequencies. When the flute plays, or a cymbal
\ crashes, we see things light up at the high end of the spectrum.
\ An FFT could be used to generate a spectrum from an input
\ musical signal.  Unfortunately, the FFT requires a lot of
\ computation which means it is difficult to make it work
\ for real time music applications.  It is typically used for
\ analysing signals that are slowly changing or for "after
\ the fact" analysis of sound. Jerry has done a great
\ job, however, at making this FFT particularly fast so I
\ am excited about seeing what it can do.
\
\ FFTs are also used for analysing medical signals, ECGs and EEGs,
\ analysing seismic records, stock market data, or almost
\ any data that has periodic components. It is an extremely
\ useful tool and I am grateful to Jerry for making this available.
\
\ This demo will only use the Fixed Point, integer, version
\ of the FFT.  It will not use Floating Point.  Fixed Point
\ is faster and is often more useful for audio applications
\ which is my main interest.
\
\ The main word for doing the FFT is FFTRC which has the following
\ stack diagram:
\
\     FFTRC  ( data-address num_points_log_base2 -- )
\
\     The data-address is the address of the data array.
\         Each point is a 4 byte, 32 bit value.
\
\     The num_points_log_base2 is the log base 2 of the
\         size of the array. The size must be some power
\         of 2, eg. 8=2**3,   or 1024=2**10.
\
\ Thus for a 1024 point array at address MY-DATA, use:
\
\         CREATE MY-DATA 1024 4* ALLOT
\         MY-DATA 10 FFTRC
\
\ The output of the FFT will be left in the data array.
\
\ Author: Phil Burk
\ Copyright 1989 Phil Burk - (FreeWare)

decimal

\ UnAssign sp: so it points to current directory.
assign sp: " "

\ load required FFT code
INCLUDE?  fft    fftinc
INCLUDE?  fftrc  fftrc
INCLUDE?  ifftcr ifftcr

\ include random number generator just for this demo.
include? choose ju:random

\ The choice between floating point and fixed point is set
\ in the file "fftcontrols".  The default is fixed point.
float_fft? 
.IF ." Woops! Floating Point version! Edit fftcontrols file." abort
.THEN

ANEW TASK-DEMO_FFT1
decimal
\ This is the log base 2 of the size of the data array.
\ In other words, 2**4=16, or 2**5=32, which is the actual size.
4 constant LOG2_FFT_SIZE

\ Allocate array for the data. Each point is 4 bytes (32 bits)
log2_fft_size 2** constant DATA_SIZE
data_size array SIGNAL-DATA
\ Size is data_size longwords, or data_size*4 bytes.

\ This integer value will correspond to 1.0 or unity.
$ 1000 constant VIRTUAL_ONE
\ Thus $ 800 would be 0.5

: CLEAR.SIGNAL  ( -- , clear signal array )
    data_size 0
    DO 0 i signal-data !
    LOOP
;

: ADD.SAWTOOTHS  ( -- , add 2 sawtooths to signal)
    0 ( starting index )
    2 0
    DO  data_size 2/ 0
        DO virtual_one data_size 2/ 2/ / i * negate
           virtual_one +  ( -- index value )
           over
           signal-data +!
           1+ ( incr index )
        LOOP
    LOOP drop
;

: ADD.PULSE ( width -- , add pulse to data)
    dup 0
    DO virtual_one i signal-data +!
    LOOP
    data_size swap
    DO 0 i signal-data +!
    LOOP
;

: ADD.RANDOM ( -- , add random noise )
    data_size 0
    DO virtual_one 2* choose virtual_one -
       i signal-data +!
    LOOP
;

: PLOT.VALUE ( value min max -- , plot a value using '*' characters)
    over - >r ( -- value min, -r- range )
    - 0 max
    40 r> */  ( normalize to 0-40 chars )
    40 min 0
    DO ascii * emit
    LOOP
;

: RANGE.DATA ( address #cells -- min max , calculate range )
    >r >r
    -1 -1 shift ( maximum integer )
    dup negate
    r> r> 0
    DO  ( min max addr )
        dup @ >r  ( get value )
        cell+ -rot  ( addr min max )
        r@ max -rot
        r> min -rot
    LOOP drop
;

: PLOT.DATA ( -- , display character plot of data )
    0 signal-data data_size range.data
    data_size 0
    DO  cr i 3 .r
        i signal-data @  10 .r 2 spaces
\
\ print appropriate number of '*'s to plot data
        i signal-data @  ( min max val )
        2 pick 2 pick plot.value
    LOOP 2drop cr
;

: RUN.FFT  ( -- , calculate FFT for input signal )
\ We are not using a precalculated bit reversal map so...
    reversal-fft off
\
    ." Input signal:" cr
    plot.data ." Hit RETURN to continue:" key drop cr
\
    ." Calculating FFT" cr
    0 signal-data  ( address of first point )
    log2_fft_size  ( 2's exponent )
    FFTRC
\
    ." Fourier Transform of Signal:" cr
    plot.data ." Hit RETURN to continue:" key drop cr
\
    ." Calculating Inverse FFT" cr
    0 signal-data  ( address of first point )
    log2_fft_size
    IFFTCR
\
    ." Reconstructed version of Signal:" cr
    plot.data
;

: SAW.FFT ( -- , run FFT on sawtooth )
    clear.signal
    add.sawtooths
    run.fft
;


: PULSE.FFT ( -- , run FFT on sawtooth )
    clear.signal
    7 add.pulse
    run.fft
;

: RANDOM.FFT ( -- , run FFT on random signal )
    clear.signal
    add.random
    run.fft
;

: NOISY.FFT ( -- , run FFT on sawtooth plus noise )
    clear.signal
    add.sawtooths
    add.random
    run.fft
;


\ By modifying the spectrum, then doing an inverse
\ transform, one can simulate the effect of a complex filter.

: RUN.FFT.FILTER  ( -- , use the FFT like a filter )
\ We are not using a precalculated bit reversal map so...
    reversal-fft off
\
    ." Input signal:" cr
    plot.data ." Hit RETURN to continue:" key drop cr
\
    ." Calculating FFT" cr
    0 signal-data  ( address of first point )
    log2_fft_size
    fftrc
\
\ Apply low pass filter to spectrum.
    data_size 4 / 0
    DO data_size i - 1- signal-data off
    LOOP
\
    ." Fourier Transform of Signal with high frequencies removed:" cr
    plot.data ." Hit RETURN to continue:" key drop cr
\
    ." Calculating Inverse FFT" cr
    0 signal-data  ( address of first point )
    log2_fft_size
    ifftcr
\
    ." Reconstructed version of Signal after Filtering:" cr
    plot.data
;

: SAW.FILTER ( -- , run FFT filter on sawtooth )
    clear.signal
    add.sawtooths
    run.fft.filter
;

: PULSE.FILTER ( -- , run FFT filter on pulse )
    clear.signal
    5 add.pulse
    run.fft.filter
;

: RANDOM.FILTER ( -- , run FFT on random signal )
    clear.signal
    add.random
    run.fft.filter
;

: NOISY.FILTER ( -- , run FFT on sawtooth plus noise )
    clear.signal
    add.sawtooths
    add.random
    run.fft.filter
;

cr cr
." Enter:  SAW.FFT  to see a sawtooth transformed, then" cr
." reconstructed.  You can also enter:" cr
."     PULSE.FFT" cr
."     RANDOM.FFT" cr
."     NOISY.FFT" cr
cr
." To see how an FFT can be used as a filter, enter:  SAW.FILTER" cr
." Read this file for more information." cr