Annotation of gforth/fft.fs, revision 1.6

1.2       anton       1: \ fast fourier transform
                      2: 
1.5       anton       3: \ Copyright (C) 2005,2007 Free Software Foundation, Inc.
1.2       anton       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
1.6     ! anton       9: \ as published by the Free Software Foundation, either version 3
1.2       anton      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
1.6     ! anton      18: \ along with this program. If not, see http://www.gnu.org/licenses/.
1.2       anton      19: 
1.1       pazsan     20: \           *** Fast Fourier Transformation ***        15may93py
                     21: 
                     22: require complex.fs
                     23: 
                     24: : Carray ( -- )  Create 0 ,  DOES> @ swap complex' + ;
                     25: Carray values
                     26: Carray 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
                     34: 8 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: 
1.4       pazsan     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 ;
1.1       pazsan     75: 
                     76: : fft  ( -- )  #points  true (fft ;
                     77: : rfft ( -- )  #points false (fft ;
                     78: 
1.3       pazsan     79: : hamming ( -- )  #points 0 DO
                     80:        I values dup z@ pi I #points fm*/ fsin f**2 f2* zscale z!
                     81:     LOOP ;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>