( DFT_Rob)   ( Not an FFT - floating point)

 ( fk=f{xk}, where xk=2*PI*k/npts ) 
2 FLD !

64 VARIABLE npts
0 VARIABLE k   ( index into waveform array)
0 VARIABLE n   ( nth coefficient)
2 VARIABLE w   ( used for eg Sin{wx})
( : robarray <BUILDS 4 * ALLOT DOES> SWAP 4 * + ; )
npts @ ( robarray) ARRAY waveform
npts @ ( robarray) ARRAY DFT->f(x)
: zero_waveform npts @ 0 DO 0 I waveform ! LOOP -1 DPL ! ;
: zero_DFT->f(x) npts @ 0 DO 0 I DFT->f(x) ! LOOP -1 DPL ! ;

: 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 Cn  ( complex array of Fourier coefficients) 

( Array Cn utilities)
: zero_Cn  npts @ 0 DO 0 0 I Cn   ( 0,0,addr a, addr b ...)
           ROT SWAP               ( 0,addr a,0,addr b ...)
           ! ! LOOP ;

: rd_Cn     ( i ... a,b)
  Cn        ( addr a, addr b ...)
  @ SWAP @  ( b,a...)
  SWAP ;    ( ... a,b)

: wr_Cn     ( a,b,i ....)
  Cn        ( a,b,addr a, addr b ...)
  ROT SWAP  ( a,addr a,b,addr b ...)  
  ! ! ;     ( ...)

: add_Cn    ( a,b,i ...)
  Cn        ( a,b,addr a,addr b ...)
  ROT SWAP  ( a,addr a,b,addr b ...)
  F+!        ( a,addr a ...)
  F+! ;      ( ....)

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

: theta   ( xk,n ... value)
    FINT F* ;

( Fourier Transform - f{x} to Cn)
: DFT
   npts @ 0 DO                ( 'n=0 TO N-1')
     npts @ 0 DO              ( 'k=0 to N-1')
      I npts @ xk             ( xk value ...) ( waveform value)
      I waveform @            ( xk,fk ... )  ( fk is waveform value)
      SWAP                    ( fk value,xk ...)
      J FINT F*               ( fk,theta ...)
      DUP FCOS SWAP FSIN      ( fk,cos,sin ...)
      3 PICK F*               ( fk,cos,fk*sin ... )
      ROT ROT F*              ( fk*sin,fk*cos ... )   ( b,a)
      SWAP FMNF               ( a,-b ...)
      J add_Cn                ( increment Cn)
     LOOP 
   LOOP ;

( Cn complex array back to f{x})
: Reverse_DFT
   zero_DFT->f(x)
   npts @ 0 DO                  ( 'n=0 TO N-1')
     npts @ 0 DO                ( 'k=0 to N-1')
        I rd_Cn                 ( i ... a,b)
        I npts @ xk             ( a,b,xk ...)
        J FINT F*               ( a,b,theta...)
        DUP FCOS SWAP FSIN      ( a,b,cos,sin ...)
        z1*z2                   ( c,d ... )
        DROP                    ( real ...)
        npts @ FINT F/          ( real/64 ...)
        J DFT->f(x) F+!
      LOOP
    LOOP ;


( Results display section ... )
( show waveform values)
: nCR 0 DO CR LOOP ;
: see_waveform  15 nCR ." Original Waveform" CR npts @ 0 DO I 8 + I DO I waveform @ F. SPACE LOOP CR 8 +LOOP ;

( show Real/Imaginary values - 8 values/line display)
: view_real CR npts @ 0 DO I 8 + I DO I rd_Cn DROP F. SPACE  LOOP CR 8 +LOOP ;
: view_imag CR npts @ 0 DO I 8 + I DO I rd_Cn F. SPACE DROP  LOOP CR 8 +LOOP ;
: view_Cn  22 nCR ." Real Cn (Cosines):" view_real 22 nCR ." Imaginary Cn (Sines):" view_imag ;

( show reverse DFT)
: view_DFT->f(x) 22 nCR ." DFT->f(x) - Reverse DFT" CR npts @ 0 DO I 8 + I DO I DFT->f(x) @ F. SPACE LOOP CR 8 +LOOP ;
: TFD Reverse_DFT view_DFT->f(x) ;

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






