Return to fft.fs CVS log

Up to [gforth] / gforth

Annotation of gforth/fft.fs, revision 1.1

1.1     ! pazsan      1: \           *** Fast Fourier Transformation ***        15may93py
        !             2: 
        !             3: require complex.fs
        !             4: 
        !             5: : Carray ( -- )  Create 0 ,  DOES> @ swap complex' + ;
        !             6: Carray values
        !             7: Carray expix
        !             8: 
        !             9: : r+ BEGIN 2dup xor -rot and dup WHILE 1 rshift REPEAT  drop ;
        !            10: : reverse  ( n -- )  2/ dup dup 2* 1
        !            11:   DO  dup I < IF  dup values I values 2dup z@ z@ z! z! THEN
        !            12:       over r+  LOOP  2drop ;
        !            13: 
        !            14: \ reverse carry add                                    23sep05py
        !            15: 8 Value #points
        !            16: : realloc ( n addr -- )
        !            17:     dup @  IF  dup @ free throw  THEN  swap allocate throw swap ! ;
        !            18: : points  ( n --- )  dup to #points dup complex' dup
        !            19:   ['] values >body realloc  2/
        !            20:   ['] expix  >body realloc
        !            21:   dup 0 DO  0e 0e I values z!  LOOP
        !            22:   1e 0e 0 expix z! 2/ dup 2/ dup 2/ dup 1+ 1
        !            23:   ?DO  pi I I' 1- 2* 2* fm*/ fsincos fswap   I expix z!  LOOP
        !            24:   ?DO  I' I - 1- expix z@ fswap    I 1+ expix z!  LOOP  dup 2/
        !            25:   ?DO  I' I -    expix z@ fswap fnegate fswap
        !            26:                                     I    expix z!  LOOP ;
        !            27: : .values  ( -- )  precision  4 set-precision
        !            28:   #points 0 DO  I values z@ z. cr  LOOP  set-precision ;
        !            29: : .expix  ( -- )   precision  4 set-precision
        !            30:   #points 2/ 0 DO  I expix z@ z. cr  LOOP  set-precision ;
        !            31: ' .values ALIAS .rvalues
        !            32: 
        !            33: \ FFT                                                  23sep05py
        !            34: 
        !            35: : z2dup+ zover zover z+ ;
        !            36: 
        !            37: : butterfly ( cexpix addr1 addr2 -- cexpix )
        !            38:   zdup over z@ z* dup z@ z2dup+ z! zr- z! ;
        !            39: : butterflies ( cexpix step off end start -- )
        !            40:   ?DO  dup I + values I values butterfly
        !            41:   over +LOOP  zdrop drop ;
        !            42: : fft-step ( n flag step steps -- n flag step )
        !            43:   0 DO  I 2 pick I' */ 2/ expix z@
        !            44:         2 pick IF  fnegate  THEN  I' 2 pick I butterflies
        !            45:   LOOP ;
        !            46: 
        !            47: \ FFT                                                  23sep05py
        !            48: 
        !            49: : (fft ( n flag -- )  swap dup reverse 1
        !            50:   BEGIN  2dup >  WHILE  dup 2* swap fft-step
        !            51:   REPEAT  2drop drop ;
        !            52: 
        !            53: : normalize ( -- )  #points dup s>f 1/f
        !            54:   0 DO  I values dup z@ 2 fpick zscale z!  LOOP  fdrop ;
        !            55: 
        !            56: : fft  ( -- )  #points  true (fft ;
        !            57: : rfft ( -- )  #points false (fft ;
        !            58: