2011-07-25 18:57:42 +08:00
|
|
|
"""
|
|
|
|
|
===========================
|
|
|
|
|
Orthogonal Matching Pursuit
|
|
|
|
|
===========================
|
|
|
|
|
|
|
|
|
|
Using orthogonal matching pursuit for recovering a sparse signal from a noisy
|
|
|
|
|
measurement encoded with a dictionary
|
|
|
|
|
"""
|
|
|
|
|
print __doc__
|
|
|
|
|
|
|
|
|
|
import pylab as pl
|
|
|
|
|
import numpy as np
|
2011-09-02 17:00:02 +08:00
|
|
|
from sklearn.linear_model import orthogonal_mp
|
|
|
|
|
from sklearn.datasets import make_sparse_coded_signal
|
2011-07-25 18:57:42 +08:00
|
|
|
|
2011-07-29 06:28:00 +08:00
|
|
|
n_components, n_features = 512, 100
|
2011-07-25 18:57:42 +08:00
|
|
|
n_atoms = 17
|
|
|
|
|
|
2011-07-29 06:28:00 +08:00
|
|
|
# generate the data
|
|
|
|
|
###################
|
2011-07-29 18:01:11 +08:00
|
|
|
|
2011-07-31 23:28:57 +08:00
|
|
|
# y = Dx
|
2011-07-29 18:01:11 +08:00
|
|
|
# |x|_0 = n_atoms
|
|
|
|
|
|
2011-08-04 21:50:47 +08:00
|
|
|
y, D, x = make_sparse_coded_signal(n_samples=1,
|
|
|
|
|
n_components=n_components,
|
|
|
|
|
n_features=n_features,
|
|
|
|
|
n_nonzero_coefs=n_atoms,
|
|
|
|
|
random_state=0)
|
2011-07-25 18:57:42 +08:00
|
|
|
|
2011-07-29 06:28:00 +08:00
|
|
|
idx, = x.nonzero()
|
2011-07-25 18:57:42 +08:00
|
|
|
|
2011-07-29 06:28:00 +08:00
|
|
|
# distort the clean signal
|
|
|
|
|
##########################
|
2011-07-31 23:28:57 +08:00
|
|
|
y_noisy = y + 0.05 * np.random.randn(len(y))
|
2011-07-25 18:57:42 +08:00
|
|
|
|
|
|
|
|
# plot the sparse signal
|
|
|
|
|
########################
|
|
|
|
|
pl.subplot(3, 1, 1)
|
2011-07-31 23:28:57 +08:00
|
|
|
pl.xlim(0, 512)
|
2011-07-25 18:57:42 +08:00
|
|
|
pl.title("Sparse signal")
|
2011-07-29 06:28:00 +08:00
|
|
|
pl.stem(idx, x[idx])
|
2011-07-25 18:57:42 +08:00
|
|
|
|
|
|
|
|
# plot the noise-free reconstruction
|
|
|
|
|
####################################
|
2011-08-03 22:06:35 +08:00
|
|
|
x_r = orthogonal_mp(D, y, n_atoms)
|
2011-07-29 06:28:00 +08:00
|
|
|
idx_r, = x_r.nonzero()
|
2011-07-25 18:57:42 +08:00
|
|
|
pl.subplot(3, 1, 2)
|
2011-07-31 23:28:57 +08:00
|
|
|
pl.xlim(0, 512)
|
2011-07-25 18:57:42 +08:00
|
|
|
pl.title("Recovered signal from noise-free measurements")
|
|
|
|
|
pl.stem(idx_r, x_r[idx_r])
|
|
|
|
|
|
|
|
|
|
# plot the noisy reconstruction
|
|
|
|
|
###############################
|
2011-08-03 22:06:35 +08:00
|
|
|
x_r = orthogonal_mp(D, y_noisy, n_atoms)
|
2011-07-29 06:28:00 +08:00
|
|
|
idx_r, = x_r.nonzero()
|
2011-07-25 18:57:42 +08:00
|
|
|
pl.subplot(3, 1, 3)
|
2011-07-31 23:28:57 +08:00
|
|
|
pl.xlim(0, 512)
|
2011-07-25 18:57:42 +08:00
|
|
|
pl.title("Recovered signal from noisy measurements")
|
|
|
|
|
pl.stem(idx_r, x_r[idx_r])
|
|
|
|
|
|
|
|
|
|
pl.subplots_adjust(0.06, 0.04, 0.94, 0.90, 0.20, 0.38)
|
|
|
|
|
pl.suptitle('Sparse signal recovery with Orthogonal Matching Pursuit',
|
|
|
|
|
fontsize=16)
|
2011-07-29 06:34:56 +08:00
|
|
|
pl.show()
|