For more DSP info, see:
Ron's DSP Web Page
For more info on Chipmunk Basic, see the
Chipmunk Basic Home 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 recursive (out-of-place) radix-2 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.