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 ;