Week 9 - Transforms¶

Problem 12.1:¶

Prove that the DFT is unitary.¶

$$ M_{fn} = \frac {1}{\sqrt {N}} e^{2\pi i\frac{fn}{N}} $$$$ M^\dagger_{f'n'} = M^*_{n'f'} = \frac {1}{\sqrt {N}} e^{-2\pi i\frac{f'n'}{N}} $$$$ \begin{split} (M^\dagger M)_{kl} &= \sum_{j=0}^{N-1} M^\dagger_{kj} M_{jl} = \sum_{j=0}^{N-1} \frac {1}{N} e^{2\pi i(l-k)\frac{j}{N}} \\ &= \frac {1}{N} \sum_{j=0}^{N-1} [\cos (2\pi (l-k)\frac{ j}{N}) + i\sin (2\pi (l-k)\frac{j}{N})] \\ &= \frac {1}{N} \sum_{j=1}^{N} [\cos (2\pi (l-k)\frac{j}{N}) + i\sin (2\pi (l-k)\frac{j}{N})] \end{split}$$

If $l = k$ we get: $$ (M^\dagger M)_{kl} = 1 $$

Otherwise, using Lagrange's trigonometric identities, we get:

$$ \begin{split} (M^\dagger M)_{kl} &= \frac{1}{N} \left [ -\frac{1}{2} + \frac{\sin ((N+1/2)2\pi \frac{(l-k)}{N})}{2 \sin \pi \frac{(l-k)}{N}} + i \frac {cos \pi \frac{(l-k)}{N} - \cos ((N+1/2)2\pi \frac{(l-k)}{N})} {2 \sin \pi \frac{(l-k)}{N}} \right ]\\ \\ &= \frac{1}{N} \left [ \frac{\sin (2\pi (l-k) + \pi \frac{(l-k)}{N}) - \sin \pi \frac{(l-k)}{N}}{2 \sin \pi \frac{(l-k)}{N}} + i \frac{\cos \pi \frac{(l-k)}{N} - \cos (2\pi (l-k) + \pi \frac{(l-k)}{N})}{2 \sin \pi \frac{(l-k)}{N}} \right ] \\ \\ &= 0 \end{split}$$

Therefore: $$ M^\dagger M = I $$

Problem 12.2:¶

Calculate the inverse wavelet transform, using Daubechies fourth-order coefficients, of a vector of length $2^{12}$, with a $1$ in the 5th and 30th places and zeros elsewhere.¶

Let $M^T$ be the inverting matrix and $\tilde{x}$ be a vector of length $2^{12}$, with a $1$ in the 5th and 30th places and zeros elsewhere. Then: $$ x = M^T \tilde{x} $$ $$ x_i = M^T_{i,5} + M^T_{i,30} $$ $$ x_5 = c_0 = \frac{1+\sqrt{3}}{4\sqrt{2}} $$ $$ x_6 = c_1 = \frac{3+\sqrt{3}}{4\sqrt{2}} $$ $$ x_7 = c_2 = \frac{3-\sqrt{3}}{4\sqrt{2}} $$ $$ x_8 = c_3 = \frac{1-\sqrt{3}}{4\sqrt{2}} $$ $$ x_{29} = c_3 = \frac{1-\sqrt{3}}{4\sqrt{2}} $$ $$ x_{30} = -c_2 = -\frac{3-\sqrt{3}}{4\sqrt{2}} $$ $$ x_{31} = c_1 = \frac{3+\sqrt{3}}{4\sqrt{2}} $$ $$ x_{32} = -c_0 = -\frac{1+\sqrt{3}}{4\sqrt{2}} $$

and zeros elsewhere.

x=np.array([0]*50, dtype=np.float64)
s = 3**.5
x[5] = 1+s
x[6] = 3+s
x[7] = 3-s
x[8] = 1-s
x[29] = 1-s
x[30] = -3+s
x[31] = 3+s
x[32] = -1-s

x *= 1/(4*2**.5)

plt.plot(range(50), x)
plt.show()

Problem 12.3:¶

Consider a measurement of a three-component vector $\vec x$, with $x_1$ and $x_2$ being drawn independently from a Gaussian distribution with zero mean and unit variance, and $x_3 = x_1 + x_2$.¶

(a) Analytically calculate the covariance matrix of $\vec x$.¶

$$ E[x_3] = E[x_1 + x_2] = E[x_1] + E[x_2] = 0 $$$$ \Sigma_{ij} = cov(x_i, x_j) = E[x_i x_j] - E[x_i] E[x_j] $$$$ \begin{split} \Sigma &= \left [ \begin{matrix} E[x_1^2] - E[x_1]^2 & E[x_1 x_2] - E[x_1] E[x_2] & E[x_1 x_3] - E[x_1] E[x_3] \\ E[x_2 x_1] - E[x_2] E[x_1] & E[x_2^2] - E[x_2]^2 & E[x_2 x_3] - E[x_2] E[x_3] \\ E[x_3 x_1] - E[x_3] E[x_1] & E[x_3 x_2] - E[x_3] E[x_2] & E[x_3^2] - E[x_3^2] \end{matrix} \right ] \\ &= \left [ \begin{matrix} Var[x_1] & E[x_1]E[x_2] & E[x_1 (x_1+x_2)] \\ E[x_2]E[x_1] & Var[x_2] & E[x_2 (x_1+x_2)] \\ E[(x_1+x_2) x_1] & E[(x_1+x_2) x_2] & E[(x_1+x_2)^2] - E[(x_1+x_2)]^2 \end{matrix} \right ] \\ &= \left [ \begin{matrix} 1 & 0 & Var[x_1]] \\ 0 & 1 & Var[x_2] \\ Var[x_1] & Var[x_2] & Var[x_1]+Var[x_2] \end{matrix} \right ] \\ &= \left [ \begin{matrix} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 2 \end{matrix} \right ] \\ \end{split} $$

(b) What are the eigenvalues?¶

$$ \begin{split} | \Sigma - \lambda I | &= \left [ \begin{matrix} 1-\lambda & 0 & 1 \\ 0 & 1-\lambda & 1 \\ 1 & 1 & 2-\lambda \end{matrix} \right ] = (1-\lambda)[(1-\lambda)(2-\lambda) - 1] - (1-\lambda)\\ &= (1-\lambda)[2-3\lambda+\lambda^2 - 2] = \lambda (1-\lambda)(\lambda-3) = 0 \end{split} $$

And the eigenvalues are $3$, $1$, and $0$.

(c) Numerically verify these results by drawing a data set from the distribution and computing the covariance matrix and eigenvalues.¶

import numpy as np

N = 10**5

X = np.zeros((3, N))
X[0] = np.random.normal(size=N)
X[1] = np.random.normal(size=N)
X[2] = X[0] + X[1]

cov = np.cov(X)
print('covariance matrix:')
print(np.round(cov, 3))

covariance matrix:
[[ 1.002 -0.002  1.   ]
 [-0.002  0.999  0.997]
 [ 1.     0.997  1.997]]

eigenvalues, eigenvectors = np.linalg.eig(cov)
print('eigenvalues:')
print(np.round(eigenvalues, 3))

eigenvalues:
[ 2.996  1.002  0.   ]

(d) Numerically find the eigenvectors of the covariance matrix, and use them to construct a transformation to a new set of variables $\vec y$ that have a diagonal covariance matrix with no zero eigenvalues. Verify this on the data set.¶

print('eigenvvectors:')
print(eigenvectors)

eigenvvectors:
[[-0.40728934  0.70765957  0.57735027]
 [-0.40920649 -0.7065527   0.57735027]
 [-0.81649583  0.00110687 -0.57735027]]

Y = np.matmul(np.transpose(eigenvectors), X)

covY = np.cov(Y)

print('new covariance matrix:')
print(np.round(covY, 5))

new covariance matrix:
[[ 2.99561 -0.      -0.     ]
 [-0.       1.00203  0.     ]
 [-0.       0.       0.     ]]

Yeigenvalues, Yeigenvectors = np.linalg.eig(covY)
print('new eigenvalues:')
print(np.round(Yeigenvalues, 5))

new eigenvalues:
[ 2.99118  0.       1.00026]

Problem 12.4:¶

Generate pairs of uniform random variables $\{s_1, s_2\}$ with each component contained in $[0, 1]$.¶

N = 1000

S = np.random.rand(2, N)

(a) Plot these data.¶

import matplotlib.pyplot as plt

plt.scatter(S[0], S[1])
plt.show()

(b) Mix them $(\vec x = A · \vec s)$ with a square matrix $A = \left [ \begin{matrix} 1 & 2 \\ 3 & 1 \end{matrix} \right ]$ and plot.

A = np.array([[1, 2], [3, 1]])

X = np.matmul(A, S)

plt.scatter(X[0], X[1])
plt.show()

(c) Make $\vec x$ zero mean, diagonalize with unit variance, and plot.¶

X_ = (X - np.mean(X, 1, keepdims=True)) 

covX = np.cov(X_)
Xeigenvalues, Xeigenvectors = np.linalg.eig(covX)

Y = np.matmul(np.transpose(Xeigenvectors), X_)
Y_ = Y / np.sqrt(Y.var(axis=1, keepdims=True))

plt.scatter(Y_[0], Y_[1])
plt.show()

print('mean: \n', np.round(np.mean(Y_, 1),10))
print('covariance: \n', np.round(np.cov(Y_), 10))

mean: 
 [ 0. -0.]
covariance: 
 [[ 1.001001  0.      ]
 [ 0.        1.001001]]