In [8]:
import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt
from numpy.fft import fftshift
%matplotlib inline

Problem 1

(a) & (b)


In [138]:
y_in = np.array([0,0,0,0,1,0,0,0])
y_in_2 = np.array([0, 0, 1, 1, 1, 1, 0, 0])
y_out = fftshift(fft.fft(fftshift(y_in)))
y_out_2 = fftshift(fft.fft(fftshift(y_in_2)))
f0, [ax0, ax1] = plt.subplots(1,2, figsize=(12,6));
ax0.plot(y_out)
ax1.plot(y_out_2.real, label ='real' )
ax1.plot(y_out_2.imag, label ='imag')
ax1.legend()
ax0.set_title('part (a)')
ax1.set_title('part (b)')


Out[138]:
<matplotlib.text.Text at 0x10d935a50>

(d)


In [139]:
x = np.arange(65)
y = np.sin(x/64.*6*np.pi)
f1, [ax0,ax1] = plt.subplots(1,2, figsize = (12,6))
ax0.plot(x,y)
ax0.set_title('Original function')
f_y= fft.fft(y)
ax1.plot(f_y)
ax1.set_title('FFT matrix')


Out[139]:
<matplotlib.text.Text at 0x10e42bcd0>

It seems the first and last points are non zero.

Problem 2

(a)


In [141]:
N = 2
H_DFT_2 = fftshift(fft.fft(fftshift(np.eye(N))))
print 'condition number is ', np.linalg.cond(H_DFT_2)

N = 4
H_DFT_4 = fftshift(fft.fft(fftshift(np.eye(N))))
print 'condition number is ', np.linalg.cond(H_DFT_4)


condition number is  1.0
condition number is  1.0

Problem 3

(a)


In [142]:
a = np.array([1.,0.,0.,0.])
b = np.arange(1.,5.)
c_conv = np.convolve(a,b)
print 'c is', c_conv
print 'length of c is ', c_conv.shape[0]


c is [ 1.  2.  3.  4.  0.  0.  0.]
length of c is  7

(b)


In [146]:
#pad 3 zeros to make it have the same number of samples as c
zero_pad = np.zeros(3)
azp = np.concatenate((a,zero_pad))
bzp = np.concatenate((b,zero_pad))
azp_FT = fft.fft(azp)
bzp_FT = fft.fft(bzp)
c_FT = azp_FT*bzp_FT
c = fft.ifft(c_FT)
print 'c from FFT is ',c
print 'difference from the true c is ', np.sum(np.abs(c-c_conv))


 c from FFT is  [  1.00000000e+00+0.j   2.00000000e+00+0.j   3.00000000e+00+0.j
   4.00000000e+00+0.j   8.62801893e-15+0.j  -1.26882631e-15+0.j
  -1.64947421e-15+0.j]
difference from the true c is  1.93178806285e-14

Now I get the same value as in (a)

Problem 4

(a)


In [147]:
H1 = np.array([[1.0, 0.20],[ 0, 1.0]])
H2 = np.array([[1.0, 0],[ 0.9, 0.1]])
H1_cond = np.linalg.cond(H1)
print 'Condition number of H1 is ', H1_cond
H2_cond = np.linalg.cond(H2)
print 'Condition number of H2 is ', H2_cond
Nx = 100
stdev = 0.2
x = np.linspace(0,1,Nx)
f = np.array([np.cos(2*np.pi*x),np.sin(2*np.pi*x)])
noise_1 = stdev*np.random.randn(f.shape[0],f.shape[1])
noise_2 = stdev*np.random.randn(f.shape[0],f.shape[1])
y1 = H1.dot(f)+noise_1
y2 = H1.dot(f)+noise_2
f_est_1 = np.linalg.pinv(H1).dot(y1)
f_est_2 = np.linalg.pinv(H2).dot(y2)
f0, [ax0,ax1] = plt.subplots(2,1, figsize=(10,10))
ax0.plot(x, f_est_1[0,:], x, f_est_1[1,:])
ax0.legend(['${cos(x)}$', '${sin(x)}$'])
ax0.set_title('Estimate from y1')
ax1.plot(x, f_est_2[0,:], x, f_est_2[1,:])
ax1.set_title('Estimate from y2')
ax1.legend(['${cos(x)}$', '${sin(x)}$'])
f_1_noise = (f_est_1-f)
f_2_noise = (f_est_2-f)
print 'Noise gain for H1 is', f_1_noise.std()/0.2
print 'H1\'s condition number is', H1_cond
print 'Noise gain for H2 is', f_2_noise.std()/0.2
print 'H2\'s condition number is', H2_cond


Condition number of H1 is  1.22099751242
Condition number of H2 is  18.144888059
Noise gain for H1 is 1.04493965756
H1's condition number is 1.22099751242
Noise gain for H2 is 30.5706647169
H2's condition number is 18.144888059

So the condition number does predict the noise gain. The higher the condition number is, the larger the noise gain


In [150]:
H1 = np.array([[1.0, 0.20],[ 0, 1.0]])
H2 = np.array([[1.0, 0],[ 0.9, 0.1]])
H1_cond = np.linalg.cond(H1)
print 'Condition number of H1 is ', H1_cond
H2_cond = np.linalg.cond(H2)
print 'Condition number of H2 is ', H2_cond
Nx = 200 # Nx is changed to 200
stdev = 0.2
x = np.linspace(0,1,Nx)
f = np.array([np.cos(2*np.pi*x),np.sin(2*np.pi*x)])
noise_1 = stdev*np.random.randn(f.shape[0],f.shape[1])
noise_2 = stdev*np.random.randn(f.shape[0],f.shape[1])
y1 = H1.dot(f)+noise_1
y2 = H1.dot(f)+noise_2
f_est_1 = np.linalg.pinv(H1).dot(y1)
f_est_2 = np.linalg.pinv(H2).dot(y2)
f0, [ax0,ax1] = plt.subplots(2,1, figsize=(10,10))
ax0.plot(x, f_est_1[0,:], x, f_est_1[1,:])
ax0.legend(['${cos(x)}$', '${sin(x)}$'])
ax0.set_title('Estimate from y1')
ax1.plot(x, f_est_2[0,:], x, f_est_2[1,:])
ax1.set_title('Estimate from y2')
ax1.legend(['${cos(x)}$', '${sin(x)}$'])
f_1_noise = (f_est_1-f)
f_2_noise = (f_est_2-f)
print 'Noise gain for H1 is', f_1_noise.std()/0.2
print 'H1\'s condition number is', H1_cond
print 'Noise gain for H2 is', f_2_noise.std()/0.2
print 'H2\'s condition number is', H2_cond


Condition number of H1 is  1.22099751242
Condition number of H2 is  18.144888059
Noise gain for H1 is 1.0098193647
H1's condition number is 1.22099751242
Noise gain for H2 is 30.3991699786
H2's condition number is 18.144888059

No. The number Nx doesn't affect the noise gain because how we sample the signal won't change the condition number of the system