Mathematical Explorations with Julia Language

We are going to start exploring the Discrete Fourier Transform (DFT) as Change of Basis.

Exploring the Fourier Basis

DFT is a powerful tool for analyzing signal. At its core, the DFT is a simple change of Basis. In this blog entry, we will explore and generate some of the Fourier Basis. We will also show that these bases are, in fact, orthogonal.

In “signal” notation, the Fourier Basis can be defined as:

Here the index k will denote the vector in the family of vectors in C^(N)

The index n will denote the element within the vector.

n, k will vary according to the rule:

Julia Language implementation

Let’s start by loading some of the packages we will need later on as we get into the coding to start exploring.

using MySignalProcessing
using Plots
using LinearAlgebra
using Symbolics

First Vector of Basis, k=0

Now we will start to generate Basis Vectors. Using the formula given above for the Basis, we plug in k=0. In Julia, math notation is very much close to the theory, so it is a pleasure to investigate and play with these concepts.

k = 0; # create a variable k and set it to zero
w₀ = [exp(((2pi/64)*im)*n*k) for n in 1:64-1] # generate 64-element vector with k=0
n = [x for x in range(0, 64 - 1)]; #create vector n

Now we can plot both Real and the Imaginary part of the vector:

p1 = plot(n , real(w₀), line = :stem, marker = :o, color=:red)
p2 = plot(n , imag(w₀), line = :stem, marker = :o, color=:black)
plot(p1, p2, layout = (2,1))

For the first vector of the Basis, we can see Real part all one, Imaginary part all zero, as we set the k to zero in the formula.

Second Vector of Basis, k=1

Here we do the same thing, but plug in k=1 in the formula and construct the second vector:

k = 1;
w₁ = [exp(((2pi/64)*im)*n*k) for n in range(0, 64-1)];
n = [x for x in range(0, 64 - 1)];

To plot, we write:

plotlyjs()
l = layout = (2,1)
p1 = plot(n , real(w₁), line = :stem, marker = :o, color=:red)
p2 = plot(n , imag(w₁), line = :stem, marker = :o, color=:black)
plot(p1, p2, layout = l)

Things are getting a bit more interesting. Now we have a full period on both the real and imaginary parts. As we go up in the basis index, the speed of this vector will increase as it goes around the unit circle.

Third Vector of the Basis, k=2

k = 2;
w₃ = [exp(((2pi/64)*im)*n*k) for n in range(0, 64-1)];
n = [x for x in range(0, 64 - 1)];
### plot
l = layout = (2,1)
p1 = plot(n , real(w₃), line = :stem, marker = :o, color=:red)
p2 = plot(n , imag(w₃), line = :stem, marker = :o, color=:black)
plot(p1, p2, layout = l)

4th Vector of the Basis, k=3

k = 3;
w₄ = [exp(((2pi/64)*im)*n*k) for n in range(0, 64-1)];
n = [x for x in range(0, 64 - 1)];

l = layout = (2,1)
p1 = plot(n , real(w₄), line = :stem, marker = :o, color=:red)
p2 = plot(n , imag(w₄), line = :stem, marker = :o, color=:black)
plot(p1, p2, layout = l)

I believe you are getting the idea: as we go up in the basis index, the vector will move faster and faster. Some of the vectors are pretty interesting. Let us jump a bit and go straight to vector with index k=16.

Vector K=16 of the Basis

k = 16;
w16  = [exp(((2pi/64)*im)*n*k) for n in range(0, 64-1)];
n = [x for x in range(0, 64 - 1)];

l = layout = (2,1)
p1 = plot(n , real(w16), line = :stem, marker = :o, color=:red)
p2 = plot(n , imag(w16), line = :stem, marker = :o, color=:black)
plot(p1, p2, layout = l)

This vector is going so fast that we can’t even discern this as sinusoidal signal. The vector is positive and then almost immediately jumps to fill negative.

Highest Frequency Signal: k=32

As a matter of fact, the signal will reach its highest possible frequency when k=32:

k = 32;;
w32 = [exp(((2pi/64)*im)*n*k) for n in range(0, 64-1)];
n = [x for x in range(0, 64 - 1)];

l = layout = (2,1)
p1 = plot(n , real(w32), line = :stem, marker = :o, color=:red)
p2 = plot(n , round.(imag(w32)), line = :stem, marker = :o, color=:black)
plot(p1, p2, layout = l)

After this the vectors are actually going to appear slower than the previous ones, going backwards. I will leave you to try some of the remaining ones.

Now these vectors form a Basis and can Span the space, so let us focus on the Orthogonality of these vectors in the basis. Fourier Basis are Orthogonal, so we will need to see how to check that.

When two vectors are Orthogonal, their dot (inner) product is zero; therefore, if we can take all the vectors in the basis and computer their dot product with each other, we should get zero.

Orthogonality of Basis Vectors

I am going to take the dot product of the basis vectors we have generated above. Let us define some symbolics to study all the different combinations we can have for the inner product operation.

using Symbolics
#defining the symbolics vectors w₀, w₁, w₃, w₄, w16, w32
@variables w₀ w₁ w₃ w₄ w16 w32 

With Julia, we can create a Set with the basis vectors we have defined:

V₁ = Set([w₀, w₁, w₃, w₄, w16, w32])

Now we can take the product of all the elements of this Set, to get all the different combinations and groupings to compare:

M = Iterators.product(V₁, V₁) |> collect

Some of these groupings for the dot product have vectors grouped with itself. A dot product of a vector with itself is going to give the magnitude of the vector. We are not really interested in this, but we can keep that in our minds as we check the result of all these groupings.

### generate all the vectors we have seen so far again:
w₀ = [exp(((2pi/64)*im)*n*0) for n in range(0, 64-1)];
w₁ = [exp(((2pi/64)*im)*n*1) for n in range(0, 64-1)];
w₃ = [exp(((2pi/64)*im)*n*3) for n in range(0, 64-1)];
w₄ = [exp(((2pi/64)*im)*n*4) for n in range(0, 64-1)];
w16 = [exp(((2pi/64)*im)*n*16) for n in range(0, 64-1)];
w32 = [exp(((2pi/64)*im)*n*32) for n in range(0, 64-1)];
# make a set of all these vectors:
V₁ = Set([w₀, w₁, w₃, w₄, w16, w32])

Now we are ready to take the product. This time, we are not using symbolics, but using the actual vectors defined above. We will collect the result in a matrix M.

M = Iterators.product(V₁, V₁) |> collect

6×6 Matrix{Tuple{Vector{ComplexF64}, Vector{ComplexF64}}}: ([1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im … 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, e-m, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im, 1.0+0.0im]) -1.0+6.12323e-16im, 1.0–7.34788e-16im, -1.0+8.57253e-16im, …

Now we can iterate over the collection M and take the dot product of each grouping of vectors.

r = [round(M[i][1] ⋅ M[i][2] , digits = 2) for i in range(1, length(M))]

The dot product of a basis vector with itself converges to 64 since the elements in the sum will be equal to 1. All the other vectors have a dot product of 0, therefore they form an Orthogonal Basis of vectors in C^(N).

Julia Language is a powerful tool for these kinds of analysis. It makes it easy to explore complex mathematical ideas since the notation is so close to the theory.

I hope you enjoyed it!