Return to fft.fs CVS log

Up to [gforth] / gforth

Added hamming filter to FFT

    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 ;