( FFT_Rob)   ( FFT - floating point)

 ( fk=f{xk}, where xk=2*PI*k/npts ) 
2 FLD ! 
1 VARIABLE yscale  ( to fit graph to vertical space)
1 VARIABLE ymax      ( ABS max Y range of plot)

1000 CONSTANT X0
addrmode 8 + @ DUP 2 5 */ CONSTANT Y1
DUP 4 5 */ CONSTANT Y2
DUP 6 5 */ CONSTANT Y3 
DUP 8 5 */ CONSTANT Y4       ( screen coords for 4 graphs)
5 / CONSTANT Y_window ( half screen height of graph windows)

( stuff from FFT by Noble: ********* with Rob mods)
0 VARIABLE DIRECTION?
: FORWARD 0 DIRECTION? ! ;
: REVERSE 1 DIRECTION? ! ;  ( was -1, for TRUE)
1 VARIABLE MMAX
0 VARIABLE ISTEP
0 VARIABLE n.bits
128 VARIABLE npts  ( was 1024)
0 VARIABLE k   ( index into waveform array)
0 VARIABLE IP
1 FINT 0 FINT zVARIABLE W
0 FINT DUP zVARIABLE WP     ( WP used to hold Cos,Sin) 

: zARRAY ( complex nos - adresses of pairs of a,b)
   <BUILDS  8 * ALLOT DOES>    ( i, addr ....)
   SWAP 8 *                    ( addr, ith offset ...)
   +                           ( addr a ...)
   DUP 4+ ;                    ( addr a, addr b...)

npts @ zARRAY waveform ( complex array of input data - imag=0)

: xk     ( k,npts .... xk)
   FINT SWAP FINT SWAP F/        ( Fk/Fnpts .. )
   2*PI F* ;                       ( .... 2*PI*k/npts)

( Array waveform utilies)
: zero_waveform npts @ 0 DO 0 0 I waveform   ( 0,0,addr a, addr b ...)
            ROT SWAP               ( 0,addr a,0,addr b ...)
            ! ! LOOP ;

: rd_waveform   ( i ... a,b)
  waveform      ( addr a, addr b ...)
  @ SWAP @  ( b,a...)
  SWAP ;    ( ... a,b) 

: wr_waveform     ( a,b,i ....)
  waveform        ( a,b,addr a, addr b ...)
  ROT SWAP    ( a,addr a,b,addr b ...)  
  ! ! ;       ( ...)

: add_waveform    ( a,b,i ...)
  waveform        ( a,b,addr a,addr b ...)
  ROT SWAP    ( a,addr a,b,addr b ...)
  F+!         ( a,addr a ...)
  F+! ;       ( ....)

( : see_waveform  CR ." Real " CR 
   npts @ 0 DO I 8 + I DO I rd_waveform DROP F. SPACE LOOP CR 8 +LOOP 
   CR ." Imaginary" CR 
   npts @ 0 DO I 8 + I DO I rd_waveform F. DROP SPACE LOOP CR 8 +LOOP ;)

: calc_nbits         ( n.bits=LOG{npts}/LOG{2})
  npts @ FINT FLOG
  2 FINT FLOG F/
  [ 0.5 FLOAT ] LITERAL F+ FIX n.bits ! ;

(    Bit Reversal Fn   )
: B.R   ( n ... n')
  0 SWAP      ( set up stack)
  n.bits @ 0 DO DUP 1 AND   ( pick out 1's bit)
   ROT 2* +                ( left shift 1 and add 1's bit)
   SWAP 2/                 ( right shift n)
   LOOP DROP ;

: BIT.REVERSE 0 IP !
   npts @ 0 DO
             I B.R IP !
             IP @ I < 0=   ( 0= equiv to NOT)
      IF 
       IP @ rd_waveform  I rd_waveform
       IP @ wr_waveform  I wr_waveform 
      ENDIF
  LOOP ;   ( end of reversal - npts times)

: z-+   ( a,b,c,d ... a-c,b-d,a+c,b+d)
   zOVER zOVER ( a,b,c,d,a,b,c,d....)
   z+          ( a,b,c,d,a+c,b+d....)
   >R >R       ( a,b,c,d............)
   z-          ( a-c,b-d............)
   R> R> ;     ( a-c,b-d,a+c,b+d....) 

: THETA Pi MMAX @ FINT F/ ;  ( THETA =Pi/MMAX)

: INIT.TRIG  THETA DUP FSIN FMNF SWAP FCOS SWAP WP z!  1 FINT 0 FINT ;    

: NEW.W ( w...w') WP z@ z1*z2 ;

: INNERLOOP   ( 1.0,0.0...w)
   ( npts I) DO MMAX @ I + IP !
           zDUP                       ( 1.0,0.0,1.0,0.0...)
           IP @ rd_waveform           ( 1.0,0.0,1.0,0.0,a,b...) ( w,w,waveform{IP}...)
           z1*z2                      ( 1.0,0.0,r1,i1...)       ( w,temp...)
           I rd_waveform              ( 1.0,0.0,r1,i1,c,d...)   ( w,temp,zi...)
           zSWAP                      ( w,zi,temp...)
           z-+                        ( w,zi-temp,zi+temp...)
           I wr_waveform              ( w,zi-temp...)
           IP @ wr_waveform           ( w...)
          ISTEP @ +LOOP ;


: FFT
  1 MMAX !         ( 1 IS MMAX )
  calc_nbits       ( npts LG2 IS N.BITS)
  BIT.REVERSE
  BEGIN
   npts @ MMAX @ >
  WHILE
   INIT.TRIG                          ( 1.0,0.0...)
   MMAX @ 2* ISTEP !
   MMAX @ 0 DO
             npts @ I INNERLOOP
             NEW.W
            LOOP zDROP
   ISTEP @ MMAX !
  REPEAT ;

( End of Noble stuff ********** )
( Inverse FFT now follows ---------------)
: IFFT
   npts @ 0 DO
          I DUP rd_waveform zconj   ( I,a,-b...)
          ROT wr_waveform
         LOOP
   FFT
   npts @ DUP FINT SWAP 0 DO DUP          ( npts,npts...)
    I DUP rd_waveform zconj  ( npts,npts,I,a,-b...)
    zSWAP                    ( npts,a,-b,npts,I...)
    >R z/f R>                ( npts,a/npts,-b/npts,I...)
    wr_waveform              ( npts...) 
   LOOP DROP ;

( Results display section ... )
( show waveform values)
: nCR 0 DO CR LOOP ;

( show Real/Imaginary values - 8 values/line display)
: see_waveform  5 VDU1 0 Y4 100 + MOVE ." Original Waveform" CR
     npts @ 0 DO I 8 + I DO I rd_waveform DROP F. SPACE LOOP CR 8 +LOOP  4 VDU1 ;

: view_real CR npts @ 0 DO I 8 + I DO I rd_waveform DROP F. SPACE  LOOP CR 8 +LOOP ;
: view_imag CR npts @ 0 DO I 8 + I DO I rd_waveform F. SPACE DROP  LOOP CR 8 +LOOP ;

: view_fft 5 VDU1 0 Y3 100 + MOVE ." Real Cn (Cosines):" view_real 
           0 Y2 100 + MOVE ." Imaginary Cn (Sines):" view_imag  4 VDU1 ;

( show reverse DFT)
: view_FFT->f(x) 5 VDU1 0 Y1 100 + MOVE ." FFT->f(x) - Reverse FFT" 
     CR npts @ 0 DO I 8 + I DO I rd_waveform DROP F. SPACE LOOP CR 8 +LOOP  4 VDU1 ;

( words for power results)    
: power CR ." power" CR npts @ 0 DO I 8 + I DO I rd_waveform 
   DUP F* SWAP DUP F* F+ F. SPACE LOOP CR 8 +LOOP ;


( Test data)
( 16 npts ! calc_nbits )
: see_wave  CR ." Complex Input" CR npts @ 0 DO I 8 + I DO I rd_waveform SWAP F. SPACE F. CR LOOP CR 8 +LOOP ;
: set_wave  CR npts @ 0 DO I FINT 0 FINT I wr_waveform LOOP ;  ( put 0,1,2... in waveform - to check reordering)
