import numpy as np
import matplotlib.pyplot as plt
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)
# Raise it to the power N to get 1
W**11
(1-3.0531133177191805e-16j)
# This is also an N-th root of unity
(W**3)
(-0.14231483827328506+0.9898214418809327j)
# Raise it to the power N to get 1
(W**3)**11
(0.9999999999999994-8.881784197001252e-16j)
# 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)
Wk = Wvec(8,2)
print(Wk)
[ 1.00000000e+00+0.00000000e+00j 2.22044605e-16+1.00000000e+00j -1.00000000e+00+4.44089210e-16j -6.66133815e-16-1.00000000e+00j 1.00000000e+00-8.88178420e-16j 1.11022302e-15+1.00000000e+00j -1.00000000e+00+1.33226763e-15j -1.55431223e-15-1.00000000e+00j]
Wk = Wvec(16,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));
plt.axis('square');
# 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)+')$');
k = 4
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)+')$');
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)+')$');
# Fourier basis vectors of DIFFERENT frequencies
Wvec(32, 2)@np.conj(Wvec(32,4))
(-1.3292605156969836e-15-8.642008586929625e-15j)
# Fourier basis vectors of the SAME frequencie
Wvec(32, 4)@np.conj(Wvec(32,4))
(31.999999999999897-4.440892098500624e-16j)
def DFTMatrix(N):
M = []
for k in range(N):
M.append(np.conj(Wvec(N,k)))
return np.array(M)
N = 5
M = DFTMatrix(N)
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]]
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]]
# Start with some random signal
f = np.random.rand(N)
print(f)
[0.83750527 0.47707399 0.33843408 0.63070474 0.62927861]
# Multiply by M, thereby computing the DFT
F = M@f
print(F)
[2.9129967 +0.j 0.39533725+0.31654758j 0.24192758-0.18850228j 0.24192758+0.18850228j 0.39533725-0.31654758j]
# Compare to the built-in DFT function
F = np.fft.fft(f)
print(F)
[2.9129967 +0.j 0.39533725+0.31654758j 0.24192758-0.18850228j 0.24192758+0.18850228j 0.39533725-0.31654758j]
Minv = np.conj(M) / N
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]]