Sliding Window

[_vince]

On the topic of music visualizers, you can see the effect of windowing

(this wasn't optimized for visualization -- everything is calculated on the fly. it ran okay on my modern laptop)

Code: [Select]

'defdbl a-z

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 (512-1), st_x_i2 (512-1)
dim st_xx_r2(512-1), st_xx_i2(512-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 - 256/2)/(128/2)
pset (mx, sh/4 - (sh/4)*exp(-z*z/2))
for i=0 to 256
z = (i - 256/2)/(128/2)
line -(mx + 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)

'gaussian window -- smooth out the edges
z = (i - 512/2)/(256/2)
st_x_r2(i) = x_r(mx*2 + 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 (mx, sh/4 - st_x_r2(0))
for i=0 to 256 -1
line -(mx + 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(), 512


'plot its short time fft rectangular
pset (0, 70+3*sh/4 - 0.005*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.005*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.005*sqr(st_xx_r2(0)*st_xx_r2(0) + st_xx_i2(0)*st_xx_i2(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.005*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 256
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/512)+" Hz"

i = m
pset ((i + c)*8, 70 + 3*sh/4), _rgb(0,256,0)
line -((i + c)*8, 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


 

[_vince]

Here is one with a non (less) truncated gaussian and no noise, note how the peaks become perfect gaussians.  The red peaks start to show a train of peaks due to the sharp transitions between sines

Code: [Select]

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


 

[bplus]

Not sure what to do with it, but as always, nice interesting graphic @_vince

[OldsCool]

Why did he make it?
Becasue he could !!!
Very Nice graphic example as well as educaational !
It would actually make a good display as part of a futuristic cockpit and would be cool to alter it so as to render a continuous random readout across the screen.

PS, I finally found someone else who uses the variable 'i' as a temporary for/next loop as much as I always did! I guess that's why I'm OldsCool...
thanks,
DB

[johnno56]

Reminded me of those old audio visualisations that WinAmp used...

Nice job, Vince!