Fast Fourier transform: Difference between revisions

=={{header|FreeBASIC}}==

<lang freebasic>'Graphic fast Fourier transform demo,

'press any key for the next image.

'131072 samples: the FFT is fast indeed.

'screen resolution

const dW = 800, dH = 600

'--------------------------------------

type samples

declare constructor (byval p as integer)

'sw = 0 forward transform

'sw = 1 reverse transform

declare sub FFT (byval sw as integer)

'draw mythical birds

declare sub oiseau ()

'plot frequency and amplitude

declare sub famp ()

'plot transformed samples

declare sub bird ()

as double x(any), y(any)

as integer fl, m, n, n2

end type

constructor samples (byval p as integer)

m = p

'number of points

n = 1 shl p

n2 = n shr 1

'real and complex values

redim x(n - 1), y(n - 1)

end constructor

'--------------------------------------

'in-place complex-to-complex FFT adapted from

'[ http://paulbourke.net/miscellaneous/dft/ ]

sub samples.FFT (byval sw as integer)

dim as double c1, c2, t1, t2, u1, u2, v

dim as integer i, j = 0, k, L, l1, l2

'bit reversal sorting

for i = 0 to n - 2

if i < j then

swap x(i), x(j)

swap y(i), y(j)

end if

k = n2

while k <= j

j -= k: k shr= 1

wend

j += k

next i

'initial cosine & sine

c1 = -1.0

c2 = 0.0

'loop for each stage

l2 = 1

for L = 1 to m

l1 = l2: l2 shl= 1

'initial vertex

u1 = 1.0

u2 = 0.0

'loop for each sub DFT

for k = 1 to l1

'butterfly dance

for i = k - 1 to n - 1 step l2

j = i + l1

t1 = u1 * x(j) - u2 * y(j)

t2 = u1 * y(j) + u2 * x(j)

x(j) = x(i) - t1

y(j) = y(i) - t2

x(i) += t1

y(i) += t2

next i

'next polygon vertex

v = u1 * c1 - u2 * c2

u2 = u1 * c2 + u2 * c1

u1 = v

next k

'half-angle sine

c2 = sqr((1.0 - c1) * .5)

if sw = 0 then c2 = -c2

'half-angle cosine

c1 = sqr((1.0 + c1) * .5)

next L

'scaling for reverse transform

if sw then

for i = 0 to n - 1

x(i) /= n

y(i) /= n

next i

end if

end sub

'--------------------------------------

'Gumowski-Mira attractors "Oiseaux mythiques"

'[ http://www.atomosyd.net/spip.php?article98 ]

sub samples.oiseau

dim as double a, b, c, t, u, v, w

dim as integer dx, y0, dy, i, k

'bounded non-linearity

if fl then

a = -0.801

dx = 20: y0 =-1: dy = 12

else

a = -0.492

dx = 17: y0 =-3: dy = 14

end if

window (-dx, y0-dy)-(dx, y0+dy)

'dissipative coefficient

b = 0.967

c = 2 - 2 * a

u = 1: v = 0.517: w = 1

for i = 0 to n - 1

t = u

u = b * v + w

w = a * u + c * u * u / (1 + u * u)

v = w - t

'remove bias

t = u - 1.830

x(i) = t

y(i) = v

k = 5 + point(t, v)

pset (t, v), 1 + k mod 14

next i

sleep

end sub

'--------------------------------------

sub samples.famp

dim as double a, s, f = n / dW

dim as integer i, k

window

k = iif(fl, dW / 5, dW / 3)

for i = k to dW step k

line (i, 0)-(i, dH), 1

next i

a = 0

k = 0: s = f - 1

for i = 0 to n - 1

a += x(i) * x(i) + y(i) * y(i)

if i > s then

a = log(1 + a / f) * 0.045

if k then

line -(k, (1 - a) * dH), 15

else

pset(0, (1 - a) * dH), 15

end if

a = 0

k += 1: s += f

end if

next i

sleep

end sub

sub samples.bird

dim as integer dx, y0, dy, i, k

if fl then

dx = 20: y0 =-1: dy = 12

else

dx = 17: y0 =-3: dy = 14

end if

window (-dx, y0-dy)-(dx, y0+dy)

for i = 0 to n - 1

k = 2 + point(x(i), y(i))

pset (x(i), y(i)), 1 + k mod 14

next i

sleep

end sub

'main

'--------------------------------------

dim as integer i, p = 17

'n = 2 ^ p

dim as samples z = p

screenres dW, dH, 4, 1

for i = 0 to 1

z.fl = i

z.oiseau

'forward

z.FFT(0)

'amplitude plot with peaks at the

'± winding numbers of the orbits.

z.famp

'reverse

z.FFT(1)

z.bird

cls

next i

end</lang>

<b>(Images only)</b>