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:
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