Powered by ViewCVS 1.0-dev
| 1 : | anton | 1.2 | \ fast fourier transform |
| 2 : | |||
| 3 : | anton | 1.5 | \ Copyright (C) 2005,2007 Free Software Foundation, Inc. |
| 4 : | anton | 1.2 | |
| 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 : | anton | 1.6 | \ as published by the Free Software Foundation, either version 3 |
| 10 : | anton | 1.2 | \ 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 : | anton | 1.6 | \ along with this program. If not, see http://www.gnu.org/licenses/. |
| 19 : | anton | 1.2 | |
| 20 : | pazsan | 1.1 | \ *** 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 : | |||
| 72 : | pazsan | 1.4 | : 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 : | pazsan | 1.1 | |
| 76 : | : fft ( -- ) #points true (fft ; | ||
| 77 : | : rfft ( -- ) #points false (fft ; | ||
| 78 : | |||
| 79 : | pazsan | 1.3 | : hamming ( -- ) #points 0 DO |
| 80 : | I values dup z@ pi I #points fm*/ fsin f**2 f2* zscale z! | ||
| 81 : | LOOP ; |
|
CVS Admin Powered by ViewCVS 1.0-dev |
ViewCVS and CVS Help |