Complex arrays & the FFT

Packed complex arrays, and a Fourier transform that works at every length and gives identical output in every engine. Run a cell to see its output beside the expected one.

A complex array

czeros allocates; indexing gives a complex number, and storing a real promotes it.

Viv [1]

a = czeros(3)
a[0] = 1.0 + 2.0i
a[1] = 5.0
show a
show len(a)

A literal with a complex in it

One complex element makes the whole literal a packed complex array; the reals promote.

Viv [2]

z = [1.0 + 2.0i, 3.0, 0.0 - 1.5i]
show z

Real, imaginary, magnitude

re, im and abs over a complex array give a REAL array of the same shape. abs is what you plot as a spectrum.

Viv [3]

z = [3.0 + 4.0i, 0.0 - 2.0i]
show re(z)
show im(z)
show abs(z)

A delta transforms to all ones

A unit impulse at index 0 transforms to all ones, which also confirms that fft does not scale its output.

Viv [4]

d = czeros(8)
d[0] = 1.0
fft(d)
show d

A pure tone shows up in two bins

A cosine at frequency 3 on 16 points puts N/2 = 8 into bins 3 and 13. The test checks that each of the other 14 bins is below 1e-12.

Viv [5]

PI = 3.141592653589793
n = 16
t = czeros(n)
for j to n {
    t[j] = cos(2.0 * PI * 3.0 * j / n)
}
fft(t)
m = abs(t)
show m[3]
show m[13]
test "every other bin is rounding noise" {
    for k to n {
        if k != 3.0 and k != 13.0 {
            assert_lt(m[k], 1.0e-12)
        }
    }
}

ifft inverts fft

ifft carries the 1/N, so the round trip returns what you started with — the same convention as numpy and FFTW.

Viv [6]

x = [1.0 + 0.5i, 2.0, 0.0 - 1.0i, 4.0]
y = [1.0 + 0.5i, 2.0, 0.0 - 1.0i, 4.0]
fft(y)
ifft(y)
test "the round trip is the identity" {
    for j to 4 {
        assert_lt(abs(y[j] - x[j]), 1.0e-14)
    }
}

Any length, not just powers of two

Seven is prime. Lengths that are not a power of two go through Bluestein's algorithm, so they cost N log N too and nothing about the call changes.

Viv [7]

PI = 3.141592653589793
n = 7
x = czeros(n)
for j to n {
    x[j] = cos(2.0 * PI * 2.0 * j / n)
}
fft(x)
m = abs(x)
show m[2]
test "a prime-length tone still lands in its bin" {
    assert_near(m[2], 3.5, 1.0e-12)
    assert_lt(m[1], 1.0e-12)
}

A real array is not a complex one

fft refuses a real array and says what to do instead. The two are different types and never mix silently.

Viv [8]

x = zeros(8)
fft(x)

expected error: `fft` needs a complex array…

← Complex numbers Linear algebra →