FFT's in Chipmunk Basic

## FFT's in Chipmunk Basic source code

For more DSP info, see: Ron's DSP Web Page

```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  fft1r()    Wrapper for a recursive FFT algorithm
rem  rec_fft()  Recursive FFT function called by fft1r()
rem  dft()      Slow plain vanilla DFT routine
rem  gmag()     Goertzel function (calculates 1 FFT bin magnitude)
rem  goertzel_mag()  Goertzel function with different calling parameters

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  2007-Nov-27  rhn

rem *** recursive out-of-place radix-2 FFT rec_fft() wrapper ***
rem	 (n,yr,yi) are local variables

sub fft1r(flag, xr, xi, m)
local n,yr,yi  : rem  local variables
n = 2^m : dim yr(n), yi(n)
rec_fft(flag, (2^m), xr,xi,0, yr,yi,0, 1,1)
mat xr = yr : mat xi = yi
end sub

rem *** recursive out-of-place radix-2 FFT ***
rem  flag = 1 for FFT, -1 for IFFT
rem  n is FFT length
rem  xr,xi are input arrays (must be dimensioned to length >= n)
rem  yr,yi are result arrays (must be dimensiond to length >= n)
rem      (kx,ky,ks,os) are sub-vector state variables
rem	 (n2,i,c,s,k1,k2,ar,ai,br,bi) are local variables
rem *** call as rec_fft(1, n, xr,xi,0, yr,yi,0, 1,1)

sub rec_fft(flag, n,xr,xi,kx,yr,yi,ky, ks,os)
local n2,i,c,s,k1,k2,ar,ai,br,bi  : rem  local variables
if (n = 1)
rem ** this does a bit-reversed-index copy and scale **
s = 1 : if (flag = -1) then s = 1/ks
yr(ky) = xr(kx) * s
yi(ky) = xi(kx) * s
else
n2 = n/2
rec_fft(flag,n2, xr,xi,kx   , yr,yi,ky,       ks*2,os)
rec_fft(flag,n2, xr,xi,kx+ks, yr,yi,ky+os*n2, ks*2,os)
for i = 0 to n2-1
c   = cos(i        * 2*pi/n)
s   = sin(i * flag * 2*pi/n)
k1 = ky+os*(i)
k2 = ky+os*(i+n2)
ar =     yr(k1)
ai =     yi(k1)
br = c * yr(k2) - s * yi(k2)
bi = c * yi(k2) + s * yr(k2)
yr(k1) = ar + br
yi(k1) = ai + bi
yr(k2) = ar - br
yi(k2) = ai - bi
next i
endif
end sub  : rem rec_fft

rem rec_fft based on an algorithm posted to comp.dsp
rem   by Steven G. Johnson on 2007-11-27

rem  *** dft subroutine ***
rem  	***  plain generic O(n^2) correlation algorithm  ***
rem  flag = 1 for fft, flag = -1 for ifft
rem  xr,xi = 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      xr, xi and stab must be arrays of length n
rem      stab must be pre-initialized with a sine lookup table

sub dft(flag, xr,xi, m)
local i,j,n,a,tmpr,tmpi,tr,ti  : rem  local variables
n = 2^m
dim tmpr(n), tmpi(n)
for j = 0 to n-1 	: rem  ** for each frequency (cycles per aperture)
tr = 0  		: rem  correlation accumulators, real & imaginary
ti = 0
for i = 0 to n-1 	: rem  ** correlate against sinusoid **
a = 2 * pi * i * j / n
tr = tr + xr(i) * cos(a)                  : rem  even sinusoid
tr = tr - xi(i) * sin(a)			: rem    complex mul
ti = ti - xr(i) * sin(a) * flag		: rem   odd sinusoid
ti = ti - xi(i) * cos(a)
next i
tmpr(j) = tr
tmpi(j) = ti
next j
a = 1 : if (flag = -1) then a = 1/n
for i = 0 to n-1 	: rem  ** scale and copy from result vector **
xr(i) = tmpr(i) * a
xi(i) = tmpi(i) * a
next i
end sub  : rem dft

rem  *** Goertzel algorithm to calculate 1-bin of DFT/FFT magnitude
rem   v is the real data vector
rem   n is the vector length
rem   f0 = frequency, sr = sample rate
rem   DFT/FFT bin number = f0 * n / sr
rem    (j,k,y1,y2,e) are local variables

sub goertzel_mag(x,n,f0,sr)
local j,k,y1,y2,e  : rem local variables
k = 2 * sin(pi * f0 / sr)
y1 = 0
y2 = 0
for j = 0 to n - 1
y2 = y2 - k * y1 + x(j)
y1 = y1 + k * y2
next j
e = sqrt(y1*y1 + y2*y2 - k*y1*y2)  : rem scaled proportional to n/2
return(e)
end sub

rem  goertzel_mag based on a rotation matrix posted to comp.dsp
rem    by Clay Turner 11/2007

rem	*** goertzel_mag() wrapper function ***
sub gmag(flag, xr, xi, m,  n, bin)
n = 2^m
bin = flag
gmag = goertzel_mag(xr, n, bin, n)
end sub

rem  Version 0.2, Copyright 2008 Ronald H Nicholson Jr.
rem  http://www.nicholson.com/rhn
rem  Consider the above code as under a BSD-style open source license,
rem  with standard disclaimers of any usefulness.

```