05-03-2022, 01:10 AM
Sliding window FFT example
This program demonstrates some of the algorithms useful for audio or other signal processing and particularly for music visualizers. Shows the effects of a short-time Fourier transform with rectangular windowing, Gaussian windowing, as well as tone detection. Features my optimized SUB rfft, or the Fast Fourier Transform -- a fast algorithm for evaluating discrete Fourier transforms. This one is particularly optimized for purely real signals. The tone detector code is meant for detecting pure sine waves in noise to high precision with spectral interpolation -- this could be useful for something like a musical instrument tuner.
This program demonstrates some of the algorithms useful for audio or other signal processing and particularly for music visualizers. Shows the effects of a short-time Fourier transform with rectangular windowing, Gaussian windowing, as well as tone detection. Features my optimized SUB rfft, or the Fast Fourier Transform -- a fast algorithm for evaluating discrete Fourier transforms. This one is particularly optimized for purely real signals. The tone detector code is meant for detecting pure sine waves in noise to high precision with spectral interpolation -- this could be useful for something like a musical instrument tuner.
Code: (Select All)
const sw = 2048
const sh = 600
dim shared pi as double
pi = 4*atn(1)
'declare sub rfft(xx_r(), xx_i(), x_r(), n)
dim x_r (sw-1), x_i (sw-1)
dim xx_r(sw-1), xx_i(sw-1)
dim st_x_r (512-1), st_x_i (512-1)
dim st_xx_r(512-1), st_xx_i(512-1)
dim st_x_r2 (sw-1), st_x_i2 (sw-1)
dim st_xx_r2(sw-1), st_xx_i2(sw-1)
dim t as double
'create signal consisting of three sinewaves in RND noise
for i=0 to sw/3-1
x_r(i) = 100*sin(2*pi*(sw*1000/44000)*i/sw) '+ (100*rnd - 50)
next
for i=sw/3 to 2*sw/3-1
x_r(i) = 100*sin(2*pi*(sw*2000/44000)*i/sw) '+ (100*rnd - 50)
next
for i=2*sw/3 to sw-1
x_r(i) = 100*sin(2*pi*(sw*8000/44000)*i/sw) '+ (100*rnd - 50)
next
screen _newimage(sw/2, sh, 32),,1,0
'plot signal
pset (0, sh/4 - x_r(0))
for i=0 to sw/2 - 1
line -(i, sh/4 - x_r(i*2)), _rgb(70,0,0)
next
line (0, sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
color _rgb(255,0,0)
_printstring (0, 0), "2048 samples of three sine waves (1 kHz, 2 kHz, 8 kHz) in RND noise sampled at 44 kHz"
rfft xx_r(), xx_i(), x_r(), sw
'plot its fft
'pset (0, 70+3*sh/4 - 0.005*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) )
for i=0 to sw/2
pset (i*2, 70 + 3*sh/4), _rgb(70,70,0)
line -(i*2, 70+3*sh/4 - 0.005*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(70,70,0)
next
line (0, 70+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
color _rgb(70,70,0)
_printstring (0, sh/2), "its entire FFT first half"
color _rgb(70,0,0)
_printstring (0, sh/2 + 16), "rectangular short time FFT"
color _rgb(0,70,0)
_printstring (0, sh/2 + 32), "gaussian short time FFT"
screen ,,0,0
pcopy 1,0
mx = 0
do
do
mx = _mousex
my = _mousey
mbl = _mousebutton(1)
mbr = _mousebutton(2)
mw = mw + _mousewheel
loop while _mouseinput
pcopy 1,0
'draw windows
if mx > sw/2-256 then mx = sw/2 - 256 - 1
if mx < 0 then mx = 0
'''rectangular window
line (mx,1)-step(256,sh/4 - 1),_rgb(255,0,0),b
'''gaussian window
z = (0 - mx - 128)/(128/2)
pset (mx, sh/4 - (sh/4)*exp(-z*z/2))
for i=0 to sw/2-1
z = (i - mx - 128)/(128/2)
line -(i, sh/4 - (sh/4)*exp(-z*z/2)),_rgb(0,255,0)
next
'take it's windowed short time FFT
for i=0 to 512-1
'rectangular window -- do nothing
st_x_r(i) = x_r(mx*2 + i)
next
for i=0 to sw - 1
'gaussian window -- smooth out the edges
z = (i - mx*2 - 256)/(128/2)
st_x_r2(i) = x_r(i)*exp(-z*z/2)
next
'''plot signal rectangular
pset (mx, sh/4 - st_x_r(0))
for i=0 to 256 -1
line -(mx + i, sh/4 - st_x_r(i*2)), _rgb(255,0,0)
next
line (0, sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
'''plot signal gaussian
pset (0, sh/4 - st_x_r2(0))
for i=0 to sw/2 - 1
line -(i, sh/4 - st_x_r2(i*2)), _rgb(0,255,0)
next
line (0, sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
rfft st_xx_r(), st_xx_i(), st_x_r(), 512
rfft st_xx_r2(), st_xx_i2(), st_x_r2(), sw
'plot its short time fft rectangular
pset (0, 70+3*sh/4 - 0.015*sqr(st_xx_r(0)*st_xx_r(0) + st_xx_i(0)*st_xx_i(0)) )
for i=0 to 128
'pset (i*8, 70 + 3*sh/4), _rgb(256,256,0)
line -(i*8, 70+3*sh/4 - 0.015*sqr(st_xx_r(i)*st_xx_r(i) + st_xx_i(i)*st_xx_i(i)) ), _rgb(256,0,0)
next
'''parabolic tone finder
dim max as double, d as double
max = 0
m = 0
for i=0 to 256
d = sqr(st_xx_r(i)*st_xx_r(i) + st_xx_i(i)*st_xx_i(i))
if d > max then
max = d
m = i
end if
next
dim c as double
dim u_r as double, u_i as double
dim v_r as double, v_i as double
u_r = st_xx_r(m - 1) - st_xx_r(m + 1)
u_i = st_xx_i(m - 1) - st_xx_i(m + 1)
v_r = 2*st_xx_r(m) - st_xx_r(m - 1) - st_xx_r(m + 1)
v_i = 2*st_xx_i(m) - st_xx_i(m - 1) - st_xx_i(m + 1)
c = (u_r*v_r + u_i*v_i)/(v_r*v_r + v_i*v_i)
color _rgb(70,70,0)
_printstring (sw/4, sh/2), "spectral parabolic interpolation tone detector"
color _rgb(255,0,0)
_printstring (sw/4, sh/2 + 16), "f_peak = "+str$((m + c)*44000/512)+" Hz"
i = m
pset ((i + c)*8, 70 + 3*sh/4), _rgb(256,256,0)
line -((i + c)*8, sh ), _rgb(256,0,0)
'plot its short time fft gaussian
pset (0, 70+3*sh/4 - 0.03*sqr(st_xx_r2(0)*st_xx_r2(0) + st_xx_i2(0)*st_xx_i2(0)) )
for i=0 to sw/2
'pset (i*8, 70 + 3*sh/4), _rgb(256,256,0)
line -(i*2, 70+3*sh/4 - 0.03*sqr(st_xx_r2(i)*st_xx_r2(i) + st_xx_i2(i)*st_xx_i2(i)) ), _rgb(0,256,0)
next
'''parabolic tone finder
max = 0
m = 0
for i=0 to sw/2
d =sqr(st_xx_r2(i)*st_xx_r2(i) + st_xx_i2(i)*st_xx_i2(i))
if d > max then
max = d
m = i
end if
next
u_r = st_xx_r2(m - 1) - st_xx_r2(m + 1)
u_i = st_xx_i2(m - 1) - st_xx_i2(m + 1)
v_r = 2*st_xx_r2(m) - st_xx_r2(m - 1) - st_xx_r2(m + 1)
v_i = 2*st_xx_i2(m) - st_xx_i2(m - 1) - st_xx_i2(m + 1)
c = (u_r*v_r + u_i*v_i)/(v_r*v_r + v_i*v_i)
color _rgb(0,256,0)
_printstring (sw/4, sh/2 + 32), "f_peak = "+str$((m + c)*44000/sw)+" Hz"
i = m
pset ((i + c)*2, 70 + 3*sh/4), _rgb(0,256,0)
line -((i + c)*2, sh ), _rgb(0,256,0)
_display
_limit 30
loop until _keyhit=27
system
sub rfft(xx_r(), xx_i(), x_r(), n)
dim w_r as double, w_i as double, wm_r as double, wm_i as double
dim u_r as double, u_i as double, v_r as double, v_i as double
log2n = log(n/2)/log(2)
for i=0 to n/2 - 1
rev = 0
for j=0 to log2n - 1
if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
next
xx_r(i) = x_r(2*rev)
xx_i(i) = x_r(2*rev + 1)
next
for i=1 to log2n
m = 2^i
wm_r = cos(-2*pi/m)
wm_i = sin(-2*pi/m)
for j=0 to n/2 - 1 step m
w_r = 1
w_i = 0
for k=0 to m/2 - 1
p = j + k
q = p + (m \ 2)
u_r = w_r*xx_r(q) - w_i*xx_i(q)
u_i = w_r*xx_i(q) + w_i*xx_r(q)
v_r = xx_r(p)
v_i = xx_i(p)
xx_r(p) = v_r + u_r
xx_i(p) = v_i + u_i
xx_r(q) = v_r - u_r
xx_i(q) = v_i - u_i
u_r = w_r
u_i = w_i
w_r = u_r*wm_r - u_i*wm_i
w_i = u_r*wm_i + u_i*wm_r
next
next
next
xx_r(n/2) = xx_r(0)
xx_i(n/2) = xx_i(0)
for i=1 to n/2 - 1
xx_r(n/2 + i) = xx_r(n/2 - i)
xx_i(n/2 + i) = xx_i(n/2 - i)
next
dim xpr as double, xpi as double
dim xmr as double, xmi as double
for i=0 to n/2 - 1
xpr = (xx_r(i) + xx_r(n/2 + i)) / 2
xpi = (xx_i(i) + xx_i(n/2 + i)) / 2
xmr = (xx_r(i) - xx_r(n/2 + i)) / 2
xmi = (xx_i(i) - xx_i(n/2 + i)) / 2
xx_r(i) = xpr + xpi*cos(2*pi*i/n) - xmr*sin(2*pi*i/n)
xx_i(i) = xmi - xpi*sin(2*pi*i/n) - xmr*cos(2*pi*i/n)
next
'symmetry, complex conj
'for i=0 to n/2 - 1
' xx_r(n/2 + i) = xx_r(n/2 - 1 - i)
' xx_i(n/2 + i) =-xx_i(n/2 - 1 - i)
'next
end sub