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