Skip to content

Instantly share code, notes, and snippets.

@hotpaw2
Created January 30, 2018 01:58
Show Gist options
  • Save hotpaw2/eefd94972870f60ed69aa09b7edba570 to your computer and use it in GitHub Desktop.
Save hotpaw2/eefd94972870f60ed69aa09b7edba570 to your computer and use it in GitHub Desktop.
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment