1\ fast fourier transform 2 3\ Copyright (C) 2005,2007 Free Software Foundation, Inc. 4 5\ This file is part of Gforth. 6 7\ Gforth is free software; you can redistribute it and/or 8\ modify it under the terms of the GNU General Public License 9\ as published by the Free Software Foundation, either version 3 10\ of the License, or (at your option) any later version. 11 12\ This program is distributed in the hope that it will be useful, 13\ but WITHOUT ANY WARRANTY; without even the implied warranty of 14\ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15\ GNU General Public License for more details. 16 17\ You should have received a copy of the GNU General Public License 18\ along with this program. If not, see http://www.gnu.org/licenses/. 19 20\ *** Fast Fourier Transformation *** 15may93py 21 22require complex.fs 23 24: Carray ( -- ) Create 0 , DOES> @ swap complex' + ; 25Carray values 26Carray expix 27 28: r+ BEGIN 2dup xor -rot and dup WHILE 1 rshift REPEAT drop ; 29: reverse ( n -- ) 2/ dup dup 2* 1 30 DO dup I < IF dup values I values 2dup z@ z@ z! z! THEN 31 over r+ LOOP 2drop ; 32 33\ reverse carry add 23sep05py 348 Value #points 35: realloc ( n addr -- ) 36 dup @ IF dup @ free throw THEN swap allocate throw swap ! ; 37: points ( n --- ) dup to #points dup complex' dup 38 ['] values >body realloc 2/ 39 ['] expix >body realloc 40 dup 0 DO 0e 0e I values z! LOOP 41 1e 0e 0 expix z! 2/ dup 2/ dup 2/ dup 1+ 1 42 ?DO pi I I' 1- 2* 2* fm*/ fsincos fswap I expix z! LOOP 43 ?DO I' I - 1- expix z@ fswap I 1+ expix z! LOOP dup 2/ 44 ?DO I' I - expix z@ fswap fnegate fswap 45 I expix z! LOOP ; 46: .values ( -- ) precision 4 set-precision 47 #points 0 DO I values z@ z. cr LOOP set-precision ; 48: .expix ( -- ) precision 4 set-precision 49 #points 2/ 0 DO I expix z@ z. cr LOOP set-precision ; 50' .values ALIAS .rvalues 51 52\ FFT 23sep05py 53 54: z2dup+ zover zover z+ ; 55 56: butterfly ( cexpix addr1 addr2 -- cexpix ) 57 zdup over z@ z* dup z@ z2dup+ z! zr- z! ; 58: butterflies ( cexpix step off end start -- ) 59 ?DO dup I + values I values butterfly 60 over +LOOP zdrop drop ; 61: fft-step ( n flag step steps -- n flag step ) 62 0 DO I 2 pick I' */ 2/ expix z@ 63 2 pick IF fnegate THEN I' 2 pick I butterflies 64 LOOP ; 65 66\ FFT 23sep05py 67 68: (fft ( n flag -- ) swap dup reverse 1 69 BEGIN 2dup > WHILE dup 2* swap fft-step 70 REPEAT 2drop drop ; 71 72: fftscale ( r -- ) 73 #points 0 DO I values dup z@ 2 fpick zscale z! LOOP fdrop ; 74: normalize ( -- ) #points s>f 1/f fftscale ; 75 76: fft ( -- ) #points true (fft ; 77: rfft ( -- ) #points false (fft ; 78 79: hamming ( -- ) #points 0 DO 80 I values dup z@ pi I #points fm*/ fsin f**2 f2* zscale z! 81 LOOP ;