Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Vince's Corner Takeout
#6
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.

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
Reply


Messages In This Thread
Vince's Corner Takeout - by bplus - 04-29-2022, 02:12 PM
RE: Vince's Corner Takeout - by vince - 04-29-2022, 09:34 PM
RE: Vince's Corner Takeout - by vince - 05-02-2022, 03:10 AM
RE: Vince's Corner Takeout - by bplus - 05-02-2022, 04:25 AM
RE: Vince's Corner Takeout - by vince - 05-02-2022, 11:16 PM
RE: Vince's Corner Takeout - by vince - 05-03-2022, 01:10 AM
RE: Vince's Corner Takeout - by bplus - 05-03-2022, 01:15 AM
RE: Vince's Corner Takeout - by vince - 05-03-2022, 04:26 AM
RE: Vince's Corner Takeout - by bplus - 05-03-2022, 03:32 PM
RE: Vince's Corner Takeout - by vince - 05-10-2022, 03:41 AM
RE: Vince's Corner Takeout - by vince - 05-10-2022, 03:57 AM
RE: Vince's Corner Takeout - by dcromley - 05-10-2022, 02:57 PM
RE: Vince's Corner Takeout - by vince - 05-10-2022, 08:14 PM
RE: Vince's Corner Takeout - by SMcNeill - 05-10-2022, 02:59 PM
RE: Vince's Corner Takeout - by vince - 05-11-2022, 01:13 AM
RE: Vince's Corner Takeout - by dcromley - 05-11-2022, 01:58 AM
RE: Vince's Corner Takeout - by vince - 06-01-2022, 09:05 AM
RE: Vince's Corner Takeout - by vince - 08-11-2022, 02:51 AM
RE: Vince's Corner Takeout - by bplus - 06-03-2022, 02:47 PM
RE: Vince's Corner Takeout - by triggered - 06-04-2022, 02:00 AM
RE: Vince's Corner Takeout - by vince - 06-07-2022, 02:02 AM
RE: Vince's Corner Takeout - by bplus - 06-07-2022, 02:15 AM
RE: Vince's Corner Takeout - by vince - 07-13-2022, 05:23 AM
RE: Vince's Corner Takeout - by BSpinoza - 07-14-2022, 04:54 AM
RE: Vince's Corner Takeout - by bplus - 07-14-2022, 04:35 PM
RE: Vince's Corner Takeout - by aurel - 08-11-2022, 01:02 PM
RE: Vince's Corner Takeout - by bplus - 08-11-2022, 04:22 PM
RE: Vince's Corner Takeout - by aurel - 08-11-2022, 05:33 PM
RE: Vince's Corner Takeout - by BSpinoza - 08-12-2022, 03:44 AM
RE: Vince's Corner Takeout - by vince - 08-11-2022, 08:42 PM
RE: Vince's Corner Takeout - by vince - 08-19-2022, 05:00 AM
RE: Vince's Corner Takeout - by bplus - 08-19-2022, 06:33 PM
RE: Vince's Corner Takeout - by vince - 08-23-2022, 10:04 PM
RE: Vince's Corner Takeout - by vince - 11-04-2022, 01:48 AM
RE: Vince's Corner Takeout - by vince - 03-31-2023, 11:07 PM
RE: Vince's Corner Takeout - by vince - 09-18-2023, 11:45 PM
RE: Vince's Corner Takeout - by Dav - 09-19-2023, 12:54 AM
RE: Vince's Corner Takeout - by bplus - 09-19-2023, 01:37 AM
RE: Vince's Corner Takeout - by GareBear - 09-19-2023, 03:56 PM
RE: Vince's Corner Takeout - by bplus - 09-19-2023, 04:47 PM
RE: Vince's Corner Takeout - by vince - 09-19-2023, 06:54 PM
RE: Vince's Corner Takeout - by bplus - 09-19-2023, 09:02 PM
RE: Vince's Corner Takeout - by vince - 01-13-2024, 07:15 PM
RE: Vince's Corner Takeout - by bplus - 01-13-2024, 07:59 PM
RE: Vince's Corner Takeout - by GareBear - 01-13-2024, 10:54 PM
RE: Vince's Corner Takeout - by vince - 02-16-2024, 04:01 AM
RE: Vince's Corner Takeout - by bplus - 02-16-2024, 02:27 PM
RE: Vince's Corner Takeout - by vince - 02-16-2024, 07:16 PM
RE: Vince's Corner Takeout - by Sprezzo - 02-17-2024, 03:04 AM
RE: Vince's Corner Takeout - by bplus - 02-17-2024, 02:44 PM
RE: Vince's Corner Takeout - by vince - 10-07-2024, 08:11 AM
RE: Vince's Corner Takeout - by bplus - 10-07-2024, 01:32 PM
RE: Vince's Corner Takeout - by vince - 10-11-2024, 12:05 AM
RE: Vince's Corner Takeout - by vince - 10-11-2024, 12:16 AM
RE: Vince's Corner Takeout - by PhilOfPerth - 10-11-2024, 03:51 AM
RE: Vince's Corner Takeout - by vince - 10-11-2024, 07:55 PM
RE: Vince's Corner Takeout - by PhilOfPerth - 10-11-2024, 11:30 PM
RE: Vince's Corner Takeout - by vince - 10-15-2024, 10:54 AM
RE: Vince's Corner Takeout - by bplus - 10-15-2024, 01:40 PM



Users browsing this thread: 14 Guest(s)