Vince's Corner Takeout

[vince]

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 sub

[dcromley]

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 Funtion

[SMcNeill]

I offered vince a spot. He declined it. It's still there, if he ever changes his mind.

[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 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'



Anyways, anything specific you're trying to do?

[vince]

B+ mods vince fine Curves code

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 Sub

[dcromley]

> vince: Anyways, anything specific you're trying to do?

That is pretty much what I expected, thanks.

[vince]

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
system

[bplus]

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 Sub

It'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?

[triggered]

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"?

[vince]

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