+- QB64 Phoenix Edition (https://staging.qb64phoenix.com)
+-- Forum: QB64 Rising (https://staging.qb64phoenix.com/forumdisplay.php?fid=1)
+--- Forum: Prolific Programmers (https://staging.qb64phoenix.com/forumdisplay.php?fid=26)
+---- Forum: bplus (https://staging.qb64phoenix.com/forumdisplay.php?fid=36)
+---- Thread: Vince's Corner Takeout (/showthread.php?tid=270)
RE: Vince's Corner Takeout - vince - 05-10-2022
Fast Fourier Transform library
This demo shows an FFT filtering example, including how to take the inverse FFT. And includes the following three SUBs:
sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
standard forward FFT. x_r and x_i are arrays containing the input signal real and imaginary parts. The output goes in arrays xx_r and xx_i. n is the size of the array and has to be a power of two (512, 1024, etc). To take the inverse, see example
sub rfft(xx_r(), xx_i(), x_r(), n)
same as above but optimized for real only signals -- there's just one input array x_r. It is about 2x faster than above for the same n
sub dft(xx_r(), xx_i(), x_r(), x_i(), n)
this is a direct and unoptimized discrete Fourier transform added in as an example. Much slower but n doesn't have to be a power of two
Filtering example:
Code: (Select All)
const sw = 512
const sh = 600
dim shared pi as double
'pi = 2*asin(1)
pi = 4*atn(1)
declare sub rfft(xx_r(), xx_i(), x_r(), n)
declare sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
declare sub dft(xx_r(), xx_i(), x_r(), x_i(), n)
dim x_r(sw-1), x_i(sw-1)
dim xx_r(sw-1), xx_i(sw-1)
dim t as double
for i=0 to sw-1
'x_r(i) = 100*sin(2*pi*62.27*i/sw) + 25*cos(2*pi*132.27*i/sw)
x_r(i) = 100*sin(0.08*i) + 25*cos(i)
x_i(i) = 0
next
'screenres sw, sh, 32
screen _newimage(sw*2, sh, 32)
'plot input signal
pset (0, sh/4 - x_r(0))
for i=0 to sw-1
line -(i, sh/4 - x_r(i)), _rgb(255,0,0)
next
line (0, sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
color _rgb(255,0,0)
_printstring (0,0), "input signal"
fft xx_r(), xx_i(), x_r(), x_i(), sw
'plot its fft
pset (0, 50+3*sh/4 - 0.01*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,255,0)
for i=0 to sw/2
line -(i*2, 50+3*sh/4 - 0.01*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(255,255,0)
next
line (0, 50+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
'set unwanted frequencies to zero
for i=50 to sw/2
xx_r(i) = 0
xx_i(i) = 0
xx_r(sw - i) = 0
xx_i(sw - i) = 0
next
'plot fft of filtered signal
pset (sw, 50+3*sh/4 - 0.01*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,255,0)
for i=0 to sw/2
line -(sw + i*2, 50+3*sh/4 - 0.01*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(0,155,255)
next
line (sw, 50+3*sh/4)-step(sw,0), _rgb(0,155,255),,&h5555
'take inverse fft
for i=0 to sw-1
xx_i(i) = -xx_i(i)
next
fft x_r(), x_i(), xx_r(), xx_i(), sw
for i=0 to sw-1
x_r(i) = x_r(i)/sw
x_i(i) = x_i(i)/sw
next
'plot filtered signal
pset (sw, sh/4 - x_r(0))
for i=0 to sw-1
line -(sw + i, sh/4 - x_r(i)), _rgb(0,255,0)
next
line (sw, sh/4)-step(sw,0), _rgb(0,255,0),,&h5555
color _rgb(0,255,0)
_printstring (sw,0), "filtered signal"
sleep
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
sub fft(xx_r(), xx_i(), x_r(), x_i(), 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)/log(2)
'bit rev copy
for i=0 to n - 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(rev)
xx_i(i) = x_i(rev)
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 - 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
end sub
sub dft(xx_r(), xx_i(), x_r(), x_i(), n)
for i=0 to n-1
xx_r(i) = 0
xx_i(i) = 0
for j=0 to n-1
xx_r(i) = xx_r(i) + x_r(j)*cos(2*pi*i*j/n) + x_i(j)*sin(2*pi*i*j/n)
xx_i(i) = xx_i(i) - x_r(j)*sin(2*pi*i*j/n) + x_i(j)*cos(2*pi*i*j/n)
next
next
end subRE: Vince's Corner Takeout - dcromley - 05-10-2022
Interesting!
Thanks, bplus ( Do we need a "QB64 Phoenix Edition › QB64 Rising › Prolific Programmers › Vince" ? ).
Maybe I'm off-base, but would your [edit: your FFT ] program process a "square pulse" ? I.e. :
Code: (Select All)
Function f1(x) ' for 0 < x < 3
If 1 <= x and x <= 2 then f1 = 1 else f1 = 0
End FuntionRE: Vince's Corner Takeout - SMcNeill - 05-10-2022
I offered vince a spot. He declined it. It's still there, if he ever changes his mind.
RE: Vince's Corner Takeout - vince - 05-10-2022
(05-10-2022, 02:57 PM)dcromley Wrote: Maybe I'm off-base, but would your [edit: your FFT ] program process a "square pulse" ? I.e. :
Code: (Select All)Function f1(x) ' for 0 < x < 3
If 1 <= x and x <= 2 then f1 = 1 else f1 = 0
End Funtion
What do you mean by process? you can take an FFT of the pulse like any other signal. Though an ideal pulse (which doesn't really exist) contains infinite frequencies to make that infinitely steep transition.
Here is a modified example showing a square pulse and low-pass filtering which removes that 'sharpness'
![[Image: zMxvpD0.png]](https://i.imgur.com/zMxvpD0.png)
Anyways, anything specific you're trying to do?
RE: Vince's Corner Takeout - vince - 05-11-2022
B+ mods vince fine Curves code
![[Image: 0GSaR6H.png]](https://i.imgur.com/0GSaR6H.png)
Code: (Select All)
Option _Explicit
_Title "b+ mods vince fine Curves code" ' b+ 2022-02-01
DefLng A-Z
Const sw = 1024, sh = 600 ' const shared everywhere
Screen _NewImage(sw, sh, 32)
_ScreenMove 150, 60 'center stage
'put 'em all here
Dim As Long n, r, mx, my, mb, omx, omy, i, j, vs
Dim As Double bx, by, t, bin
Dim k$
ReDim x(n) As Long, y(n) As Long
vs = _NewImage(sw, sh, 32) ' vs for virtual screen
r = 5 'gap checker?
Do
Cls
k$ = InKey$
If k$ = "c" Then
_Dest vs
Line (0, 0)-(sw, sh), &HFF000000, BF
_Dest 0
Cls
End If
_PutImage , vs, 0
While _MouseInput: Wend ' poll mouse update mouse variables
mx = _MouseX: my = _MouseY: mb = _MouseButton(1)
If mb Then
n = 1
ReDim _Preserve x(n)
ReDim _Preserve y(n)
x(0) = mx - sw / 2
y(0) = sh / 2 - my
PSet (mx, my)
Do While mb
While _MouseInput: Wend ' poll mouse update mouse variables
mx = _MouseX: my = _MouseY: mb = _MouseButton(1)
Line -(mx, my), _RGB(30, 30, 30)
If (mx - omx) ^ 2 + (my - omy) ^ 2 > r ^ 2 Then
circlef mx, my, 3, _RGB(30, 30, 30)
omx = mx
omy = my
x(n) = mx - sw / 2
y(n) = sh / 2 - my
n = n + 1
ReDim _Preserve x(n)
ReDim _Preserve y(n)
End If
_Display
'_Limit 30
Loop
'close the contour
'x(n) = x(0)
'y(n) = y(0)
'n = n + 1
'redim _preserve x(n)
'redim _preserve y(n)
'redraw spline
'pset (sw/2 + x(0), sh/2 - y(0))
'for i=0 to n
'line -(sw/2 + x(i), sh/2 - y(i)), _rgb(255,0,0)
'circlef sw/2 + x(i), sh/2 - y(i), 3, _rgb(255,0,0)
'next
_Dest vs
PSet (sw / 2 + x(0), sh / 2 - y(0))
For t = 0 To 1 Step 0.001
bx = 0
by = 0
For i = 0 To n
bin = 1
For j = 1 To i
bin = bin * (n - j) / j
Next
bx = bx + bin * ((1 - t) ^ (n - 1 - i)) * (t ^ i) * x(i)
by = by + bin * ((1 - t) ^ (n - 1 - i)) * (t ^ i) * y(i)
Next
Line -(sw / 2 + bx, sh / 2 - by), _RGB(255, 0, 0)
Next
_Dest 0
End If
_Display
_Limit 30
Loop Until _KeyHit = 27
System
Sub circlef (x As Long, y As Long, r As Long, c As Long)
Dim As Long x0, y0, e
x0 = r
y0 = 0
e = -r
Do While y0 < x0
If e <= 0 Then
y0 = y0 + 1
Line (x - x0, y + y0)-(x + x0, y + y0), c, BF
Line (x - x0, y - y0)-(x + x0, y - y0), c, BF
e = e + 2 * y0
Else
Line (x - y0, y - x0)-(x + y0, y - x0), c, BF
Line (x - y0, y + x0)-(x + y0, y + x0), c, BF
x0 = x0 - 1
e = e - 2 * x0
End If
Loop
Line (x - r, y)-(x + r, y), c, BF
End SubRE: Vince's Corner Takeout - dcromley - 05-11-2022
> vince: Anyways, anything specific you're trying to do?
That is pretty much what I expected, thanks.
RE: Vince's Corner Takeout - vince - 06-01-2022
Chaotic Scattering
demo of the Gaspard-Rice system: https://en.wikipedia.org/wiki/Chaotic_scattering
left click to reposition laser
mouse wheel to change radius of reflectors
Code: (Select All)
deflng a-z
sw = 800
sh = 600
dim shared pi as double
pi = 4*atn(1)
dim as double t, a, b, a1, a2
dim as double x, y, x0, y0, x1, y1, dx, dy
r = 150
rr0 = 110
sx = 0
sy = sh/2
screen _newimage(sw, sh, 32)
do
do while _mouseinput
mw = mw + _mousewheel
loop
mx = _mousex
my = _mousey
mb = _mousebutton(1)
rr = rr0 - mw
if mb then
do while mb
do while _mouseinput
loop
mb = _mousebutton(1)
loop
valid = -1
for b = 0 to 2*pi step 2*pi/3
x1 = r*cos(b) + sw/2
y1 = r*sin(b) + sh/2
dx = mx - x1
dy = my - y1
if dx*dx + dy*dy < rr*rr then
valid = 0
exit for
end if
next
if valid then
sx = mx
sy = my
end if
end if
if mx<>old_mx or my<>old_my or mw<>old_mw then
line (0,0)-(sw,sh), _rgb(0,0,0), bf
'locate 1,1
'? mx, my, mw, mb
for b = 0 to 2*pi step 2*pi/3
circle (r*cos(b) + sw/2, r*sin(b) + sh/2), rr
next
a = _atan2(my - sy, mx - sx)
x0 = sx
y0 = sy
for t = 0 to 1000
x = t*cos(a) + x0
y = t*sin(a) + y0
for b = 0 to 2*pi step 2*pi/3
if x >= 0 and x < sw and y >=0 and y < sh then
x1 = r*cos(b) + sw/2
y1 = r*sin(b) + sh/2
dx = x - x1
dy = y - y1
if dx*dx + dy*dy < rr*rr then
a1 = _atan2(dy, dx)
a2 = 2*a1 - a - pi
line (x0, y0)-(x, y), _rgb(233,205,89)
x0 = x
y0 = y
a = a2
t = 0
exit for
end if
end if
next
next
line (x0, y0)-(x, y), _rgb(233,205,89)
end if
old_mx = mx
old_my = my
old_mw = mw
_display
_limit 50
loop until _keyhit = 27
systemRE: Vince's Corner Takeout - bplus - 06-03-2022
I think vince asked me to post this, it was a mod of his scattering that allowed laser to be set anywhere (by click of mouse) and reflect off a random arrangement of circles. As you move mouse around, the laser points at slightly different angles causing radical changes in reflection outcomes:
Code: (Select All)
_Title "*** Chaotic Scattering *** by vince and mod by bplus 2018-02-15 click mouse to reset LASER"
DefInt A-Z
Randomize Timer
Const sw = 1200
Const sh = 700
Dim Shared qb(15) As _Integer64
qb(0) = &HFF000000
qb(1) = &HFF000088
qb(2) = &HFF008800
qb(3) = &HFF008888
qb(4) = &HFF880000
qb(5) = &HFF880088
qb(6) = &HFF888800
qb(7) = &HFFCCCCCC
qb(8) = &HFF888888
qb(9) = &HFF0000FF
qb(10) = &HFF00FF00
qb(11) = &HFF00FFFF
qb(12) = &HFFFF0000
qb(13) = &HFFFF00FF
qb(14) = &HFFFFFF00
qb(15) = &HFFFFFFFF
Const nCircs = 25
Const r = 150
Const maxr = 100
Type circles
x As Integer
y As Integer
r As Integer
c As _Integer64
End Type
Dim Shared cs(nCircs) As circles
Dim i As Integer
Dim c As Integer
Dim ck As Integer
For i = 1 To nCircs
cs(i).r = Rnd * (maxr - 20) + 20
cs(i).c = qb(Int(Rnd * 15) + 1)
If i > 1 Then
ck = 0
While ck = 0
cs(i).x = Int(Rnd * (sw - 2 * cs(i).r)) + cs(i).r
cs(i).y = Int(Rnd * (sh - 2 * cs(i).r)) + cs(i).r
ck = 1
For c = 1 To i - 1
If ((cs(i).x - cs(c).x) ^ 2 + (cs(i).y - cs(c).y) ^ 2) ^ .5 < cs(i).r + cs(c).r Then ck = 0: Exit For
Next
Wend
Else
cs(i).x = Int(Rnd * (sw - 2 * cs(i).r)) + cs(i).r
cs(i).y = Int(Rnd * (sh - 2 * cs(i).r)) + cs(i).r
End If
Next
Dim t As Double
Dim a As Double, b As Double
Dim a1 As Double, a2 As Double
Dim x As Double, y As Double
Dim x0 As Double, y0 As Double
Dim x1 As Double, y1 As Double
Screen _NewImage(sw, sh, 32)
_ScreenMove 100, 20
'find a place not inside a circle
xx = sw / 2
yy = sh / 2
While checkxy%(xx, yy) = 0
xx = Int(Rnd * (sw - 2 * maxr)) + maxr
yy = Int(Rnd * (sh - 2 * maxr)) + maxr
Wend
Do
If Len(InKey$) Then
_Delay 5 'to get dang screen shot
Else
'get mouse x, y if click
Do
mx = _MouseX
my = _MouseY
mb = _MouseButton(1)
Loop While _MouseInput
End If
'cls with Fellippes suggestion
Line (0, 0)-(sw, sh), _RGBA32(0, 0, 0, 30), BF
'draw circles
For c = 1 To nCircs
Color cs(c).c
fcirc cs(c).x, cs(c).y, cs(c).r
Next
'if click make sure click was not inside one of the circles
If mb Then
Do While mb
Do
mb = _MouseButton(1)
Loop While _MouseInput
Loop
f = checkxy%(mx, my)
If f Then
xx = mx
yy = my
f = -1
End If
End If
x0 = xx
y0 = yy
a = _Atan2(my - yy, mx - xx)
t = 0
Do
t = t + 1
x = t * Cos(a) + x0
y = t * Sin(a) + y0
If x < 0 Or x > sw Or y < 0 Or y > sh Then Exit Do
For c = 1 To nCircs
If (x - cs(c).x) ^ 2 + (y - cs(c).y) ^ 2 < cs(c).r * cs(c).r Then
a1 = _Atan2(y - cs(c).y, x - cs(c).x)
a2 = 2 * a1 - a - _Pi
Line (x0, y0)-(x, y), qb(14)
x0 = x
y0 = y
a = a2
t = 0
Exit For
End If
Next
Loop
Line (x0, y0)-(x, y), qb(14)
_Display
_Limit 50
Loop Until _KeyHit = 27
System
Function checkxy% (x, y)
Dim c As Integer
For c = 1 To nCircs
If (x - cs(c).x) ^ 2 + (y - cs(c).y) ^ 2 < cs(c).r * cs(c).r Then checkxy% = 0: Exit Function
Next
checkxy% = 1
End Function
'Steve McNeil's copied from his forum note: Radius is too common a name
Sub fcirc (CX As Long, CY As Long, R As Long)
Dim subRadius As Long, RadiusError As Long
Dim X As Long, Y As Long
subRadius = Abs(R)
RadiusError = -subRadius
X = subRadius
Y = 0
If subRadius = 0 Then PSet (CX, CY): Exit Sub
' Draw the middle span here so we don't draw it twice in the main loop,
' which would be a problem with blending turned on.
Line (CX - X, CY)-(CX + X, CY), , BF
While X > Y
RadiusError = RadiusError + Y * 2 + 1
If RadiusError >= 0 Then
If X <> Y + 1 Then
Line (CX - Y, CY - X)-(CX + Y, CY - X), , BF
Line (CX - Y, CY + X)-(CX + Y, CY + X), , BF
End If
X = X - 1
RadiusError = RadiusError - X * 2
End If
Y = Y + 1
Line (CX - X, CY - Y)-(CX + X, CY - Y), , BF
Line (CX - X, CY + Y)-(CX + X, CY + Y), , BF
Wend
End SubIt's a nice effect and might be used by Indiana Jones to unlock a treasure with a beam of light ;-))
Or maybe laser printers work like this?
RE: Vince's Corner Takeout - triggered - 06-04-2022
This is the most impressive thing I've seen in qb64 to date, the responsiveness and correctness of this is amazing. Did you do this using snell's law? Where can I learn about doing reflections like this?
P.S. Is this the same as "ray tracing"?
RE: Vince's Corner Takeout - vince - 06-07-2022
Flower of Justice
this is a JB original
Code: (Select All)
'blossoming flower of justice
const sw = 800
const sh = 600
dim shared pi as double
pi = 4*atn(1)
screen _newimage(sw, sh, 32)
r = 100
do
for a = 0 to 1 step 0.01
cls
fcirc sw/2, sh/2, a*r + (1-a)*1.5*r, a, 3
_display
_limit 5
next
_delay 3
for a = 1 to 0 step -0.01
cls
fcirc sw/2, sh/2, a*r + (1-a)*1.5*r, a, 3
_display
_limit 5
next
_delay 3
loop
system
sub fcirc (x, y, r, a, n)
if not n > 0 then exit sub
for t=0 to 2*pi step 2*pi/6
xx = x + r*cos(t)
yy = y + r*sin(t)
circle (xx, yy), r
fcirc xx, yy, a*r, a, n - 1
next
end sub