2. Using iterative algo to solve NUS problem

In [1]:

    

# load environment
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import dft, hadamard, pinv
from numpy.fft import fft, ifft
from math import sqrt

from sklearn.linear_model import OrthogonalMatchingPursuit


    

In [2]:

    

# parameters
n_components = 1024    # N: components of the spectrum,
n_features = 200       # M: length of observation
n_nonzero_coefs = 10   # K: nb of non-zero components in the direct domain. 

# create peaks
randindex =  np.random.permutation(n_components)[:n_nonzero_coefs]

# spectrum
x = np.zeros(n_components)
x[randindex] = np.arange(n_nonzero_coefs) #
plt.plot(x)


    

    
Out[2]:





[<matplotlib.lines.Line2D at 0x107c2da90>]






    

In [11]:

    

# Choose the transform

# Hadamard
D = hadamard(n_components)

# Fourier
#D = dft(n_components)

# Random
#D = np.random.randn(n_components,n_components)/np.sqrt(n_components)

# combining Hadamard and random
#D1 = hadamard(n_components)
#D2 = np.random.randn(n_components,n_components)/np.sqrt(n_components)
#D = np.dot(D1,D2)

# For some transform, there is no difference between truncation and random sampling
# create peaks

sampling =  np.random.permutation(n_components)[:n_features]
Dt = D[sampling,:]

# create acquired data
spec = np.dot(Dt,x)
# eventually add some noise
noise = 0.0
spec += noise*max(spec)*np.random.randn(n_features)

print Dt.shape, x.shape, spec.shape

plt.plot(x)
plt.title("original data")
plt.figure()
plt.plot(spec)
plt.title("acquired data")


    

    




(200, 1024) (1024,) (200,)






    
Out[11]:





<matplotlib.text.Text at 0x10a8b4510>






    













    

In [12]:

    

# Least Square Solution
Dm1 = pinv(Dt)
lsq = np.dot(Dm1, spec)
plt.plot(lsq)


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0x10ce896d0>]






    

In [13]:

    

######## Reconstruction

omp = OrthogonalMatchingPursuit()
omp.fit(Dt, spec)
coef = omp.coef_
idx_r, = coef.nonzero()

plt.title("signal recovered")
plt.plot(coef)
plt.xlim(xmin = 0, xmax=n_components)

plt.figure()
plt.title("original")
plt.plot(x)
h=plt.xlim(xmin = 0, xmax=n_components)


    

In [ ]: