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>