fft1.bas

rem Some tutorial DSP subroutines for Chipmunk Basic
rem
rem  genstab() Generates a sine table for use by the fft() subroutine
rem  fft()      Implementation of a common in-place DIT FFT algorithm

rem *** genstab(m) - a subroutine to generate sine table of length 2^m ***
rem	*** must be called before calling fft() for
rem         the first time, and if the length 2^m changes

sub genstab(m)
  local i,a,n  : rem local variables
  n = 2^m
  dim stab(n)
  for i=0 to n-1
    a = 2 * pi * i/n
    stab(i) = sin(a)
  next i
end sub

rem  *** fft subroutine ***
rem     ***  in-place (presort) decimation-in-time radix-2 FFT
rem 	***  uses a pre-generated twiddle/trig factor table
rem  flag = 1 for fft, flag = -1 for ifft
rem  vr,vi = real & imaginary part of data, both 1d vector/arrays
rem  m = log2(n)  e.g. the power of 2 of the fft length n
rem    vr, vi and stab must all be arrays of length n
rem    stab() array must be pre-initialized with the sine lookup table
rem    i,j,k,n,ki,kr,istep,astep,a,wr,wi,tr,ti,qr,qi are local variables

sub fft(flag, vr,vi,m)
  local i,j,k,n,ki,kr,istep,astep,a,wr,wi,tr,ti,qr,qi 
  n = 2^m
  rem  		*** shuffle data using bit reversed addressing ***
  for k=0 to n-1
    rem  	*** generate a bit reversed address vr k ***
    ki = k
    kr = 0
    for i=1 to m
      kr = kr * 2 	: rem  **  left shift result kr by 1 bit
      if ki mod 2 = 1 then kr = kr + 1 
      ki = int(ki/2)  	: rem  **  right shift temp ki by 1 bit
    next i
    rem	 	*** swap data vr k to bit reversed address kr
    if (kr > k)
      tr = vr[kr] : vr[kr] = vr[k] : vr[k] = tr
      ti = vi[kr] : vi[kr] = vi[k] : vi[k] = ti
    endif
  next k

  rem  		*** do fft butterflys in place ***
  istep = 2
  while ( istep <= n ) : rem  *** layers 2,4,8,16, ... ,n ***
    is2 = istep/2
    astep = n/istep
    for km = 0 to (is2)-1  	: rem 	*** outer row loop ***
      a  = km * astep		: rem  twiddle angle index
      wr =  stab(a+(n/4))	: rem  get sin from table lookup
      wi =  stab(a) 	     	: rem pos for fft , neg for ifft
      rem   stab(a) == sin(2 * pi * km / istep)  by table lookup
      if (flag = -1) then wi = -wi
      for ki = 0 to (n - istep) step istep : rem  *** inner column loop ***
        i = km + ki
        j = (is2) + i
        tr = wr * vr[j] - wi * vi[j] 	: rem ** complex multiply **
        ti = wr * vi[j] + wi * vr[j]
        qr = vr[i]
        qi = vi[i]
        vr[j] = qr - tr
        vi[j] = qi - ti
        vr[i] = qr + tr
        vi[i] = qi + ti
      next ki
    next km
    istep = istep * 2
  wend
  rem  *** scale fft (or ifft) by n ***
  if (flag = -1)  : rem  ifft scaling or test flag = 1 for fft scaling
    a = 1/n
    for i=0 to n-1 : vr(i) = vr(i)*a : vi(i) = vi(i)*a : next i
  endif
end sub  : rem fft

rem -- source code originally from http://www.nicholson.com/dsp.fft1.html
rem -- BSD 2-clause license