1: \ fast fourier transform
2:
3: \ Copyright (C) 2005 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 2
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, write to the Free Software
19: \ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA.
20:
21: \ *** Fast Fourier Transformation *** 15may93py
22:
23: require complex.fs
24:
25: : Carray ( -- ) Create 0 , DOES> @ swap complex' + ;
26: Carray values
27: Carray expix
28:
29: : r+ BEGIN 2dup xor -rot and dup WHILE 1 rshift REPEAT drop ;
30: : reverse ( n -- ) 2/ dup dup 2* 1
31: DO dup I < IF dup values I values 2dup z@ z@ z! z! THEN
32: over r+ LOOP 2drop ;
33:
34: \ reverse carry add 23sep05py
35: 8 Value #points
36: : realloc ( n addr -- )
37: dup @ IF dup @ free throw THEN swap allocate throw swap ! ;
38: : points ( n --- ) dup to #points dup complex' dup
39: ['] values >body realloc 2/
40: ['] expix >body realloc
41: dup 0 DO 0e 0e I values z! LOOP
42: 1e 0e 0 expix z! 2/ dup 2/ dup 2/ dup 1+ 1
43: ?DO pi I I' 1- 2* 2* fm*/ fsincos fswap I expix z! LOOP
44: ?DO I' I - 1- expix z@ fswap I 1+ expix z! LOOP dup 2/
45: ?DO I' I - expix z@ fswap fnegate fswap
46: I expix z! LOOP ;
47: : .values ( -- ) precision 4 set-precision
48: #points 0 DO I values z@ z. cr LOOP set-precision ;
49: : .expix ( -- ) precision 4 set-precision
50: #points 2/ 0 DO I expix z@ z. cr LOOP set-precision ;
51: ' .values ALIAS .rvalues
52:
53: \ FFT 23sep05py
54:
55: : z2dup+ zover zover z+ ;
56:
57: : butterfly ( cexpix addr1 addr2 -- cexpix )
58: zdup over z@ z* dup z@ z2dup+ z! zr- z! ;
59: : butterflies ( cexpix step off end start -- )
60: ?DO dup I + values I values butterfly
61: over +LOOP zdrop drop ;
62: : fft-step ( n flag step steps -- n flag step )
63: 0 DO I 2 pick I' */ 2/ expix z@
64: 2 pick IF fnegate THEN I' 2 pick I butterflies
65: LOOP ;
66:
67: \ FFT 23sep05py
68:
69: : (fft ( n flag -- ) swap dup reverse 1
70: BEGIN 2dup > WHILE dup 2* swap fft-step
71: REPEAT 2drop drop ;
72:
73: : normalize ( -- ) #points dup s>f 1/f
74: 0 DO I values dup z@ 2 fpick zscale z! LOOP fdrop ;
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 ;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>