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:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>