Fourier Basis Vectors¶

In [2]:
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=4, suppress=True)

N-th roots of unity¶

Consider taking the $N^\mathrm{th}$ root of 1. It can be written $$ W_N = e^{2\pi i\frac{1}{N}} $$ In fact, there are $N$ such roots, $$ W_N = e^{2\pi i\frac{n}{N}} \quad \text{for} \ n=0, \ldots , N-1 $$

In [3]:
N = 11
W = np.exp(2.j*np.pi/N)
print(f'One of the {N}-th roots of unity is {W}')

One of the 11-th roots of unity is (0.8412535328311812+0.5406408174555976j)
In [4]:
# Raise it to the power N to get 1
W**11
Out[4]:
(1-3.0531133177191805e-16j)
In [5]:
# This is also an N-th root of unity
(W**3)
Out[5]:
(-0.14231483827328506+0.9898214418809327j)
In [6]:
# Raise it to the power N to get 1
(W**3)**11
Out[6]:
(0.9999999999999994-8.881784197001252e-16j)
In [7]:
# This function produces the Fourier basis vector W_N(k)
def Wvec(N,k):
    W = np.exp(2.j*np.pi/N)
    Wk = []
    for n in range(N):
        Wk.append(W**(n*k))
    return np.array(Wk)
In [9]:
Wk = Wvec(8,3)
print(Wk)

[ 1.    +0.j     -0.7071+0.7071j -0.    -1.j      0.7071+0.7071j
 -1.    +0.j      0.7071-0.7071j  0.    +1.j     -0.7071-0.7071j]
In [13]:
Wk = Wvec(8,3)
pts = []
for idx in range(len(Wk)):
    pts.append([np.real(Wk[idx]), np.imag(Wk[idx])])
    plt.plot(np.real(Wk[idx]), np.imag(Wk[idx]), 'o', color=idx/len(Wk)*np.ones(3));
pts = np.array(pts)
plt.plot(pts[:,0], pts[:,1],'-', color=(0.8,0.8,0.8,0.5));
circle = plt.Circle((0,0), 1, edgecolor='grey', linestyle='--', fill=False)
plt.gca().add_patch(circle)
plt.axis('square');
In [14]:
# If we plot the real and imaginary components of W_N(k), they are cosines and sines
# of frequency k.
k = 3
N = 128
Wk = Wvec(N, k)
plt.figure(figsize=(12,4));
plt.plot(np.real(Wk), 'bo');
plt.plot(np.imag(Wk), 'o', color='lightblue');
plt.title(r'$W_{'+str(N)+'}('+str(k)+')$');
In [16]:
k = 0
Wk = Wvec(N, k)
plt.figure(figsize=(12,4));
plt.plot(np.real(Wk), 'bo');
plt.plot(np.imag(Wk), 'o', color='lightblue');
plt.title(r'$W_{'+str(N)+'}('+str(k)+')$');
In [17]:
k = 124
Wk = Wvec(N, k)
plt.figure(figsize=(12,4));
plt.plot(np.real(Wk), 'bo');
plt.plot(np.imag(Wk), 'o', color='lightblue');
plt.title(r'$W_{'+str(N)+'}('+str(k)+')$');

Let's try out the orthogonality¶

In [18]:
# Fourier basis vectors of DIFFERENT frequencies
Wvec(32, 2)@np.conj(Wvec(32,4))
Out[18]:
(-1.2678496299609423e-15-7.753399305115232e-15j)
In [19]:
# Fourier basis vectors of the SAME frequencie
Wvec(32, 4)@np.conj(Wvec(32,4))
Out[19]:
(31.999999999999893+0j)

DFT Matrix¶

In [20]:
def DFTMatrix(N):
    M = []
    for k in range(N):
        M.append(np.conj(Wvec(N,k)))
    return np.array(M)
In [21]:
N = 5
M = DFTMatrix(N)
In [22]:
print( np.round(M, decimals=3) )

[[ 1.   -0.j     1.   -0.j     1.   -0.j     1.   -0.j     1.   -0.j   ]
 [ 1.   -0.j     0.309-0.951j -0.809-0.588j -0.809+0.588j  0.309+0.951j]
 [ 1.   -0.j    -0.809-0.588j  0.309+0.951j  0.309-0.951j -0.809+0.588j]
 [ 1.   -0.j    -0.809+0.588j  0.309-0.951j  0.309+0.951j -0.809-0.588j]
 [ 1.   -0.j     0.309+0.951j -0.809+0.588j -0.809-0.588j  0.309-0.951j]]
In [23]:
print( np.round(M@np.conj(M), decimals=2) )

[[ 5.+0.j  0.+0.j  0.+0.j  0.+0.j  0.-0.j]
 [ 0.+0.j  5.+0.j  0.+0.j -0.+0.j  0.-0.j]
 [ 0.-0.j  0.-0.j  5.+0.j  0.+0.j  0.+0.j]
 [ 0.-0.j -0.-0.j  0.-0.j  5.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.-0.j  0.-0.j  5.+0.j]]

Use the matrix to compute the DFT¶

In [24]:
# Start with some random signal
f = np.random.rand(N)
print(f)

[0.3044 0.9141 0.9844 0.5337 0.3221]
In [25]:
# Multiply by M, thereby computing the DFT
F = M@f
print(F)

[ 3.0587+0.j     -0.5417-0.8279j -0.2266+0.0806j -0.2266-0.0806j
 -0.5417+0.8279j]
In [26]:
# Compare to the built-in DFT function
F = np.fft.fft(f)
print(F)

[ 3.0587+0.j     -0.5417-0.8279j -0.2266+0.0806j -0.2266-0.0806j
 -0.5417+0.8279j]

We can compute its inverse¶

In [27]:
Minv = np.conj(M) / N
In [28]:
# Show that M and Minv are inverses of each other
print(np.round(M@Minv, decimals=1))

[[ 1.+0.j -0.-0.j  0.+0.j  0.+0.j  0.-0.j]
 [-0.+0.j  1.+0.j  0.+0.j  0.+0.j  0.-0.j]
 [ 0.-0.j  0.-0.j  1.+0.j  0.+0.j  0.+0.j]
 [ 0.-0.j  0.-0.j  0.-0.j  1.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.-0.j  0.-0.j  1.-0.j]]
In [ ]: