library(bcs)

# size of the basis function expansion
N <- 64
# generate sparse coefficient vector
w <- rep(0,N)
w[sample(1:N,10)] <- runif(10,-1,1)
# create wavelet basis trasform matrix
wavelet.basis <- WaveletBasis(N)
# generate actual signal
signal <- wavelet.basis%*%w
# now we try and recover 'w' and 'signal' from samples
num_samps <- 25
# create random measurement matrix
measure.mat <- matrix(runif(num_samps*N),num_samps,N)
measure.mat <- measure.mat/matrix(rep(sqrt(apply(measure.mat^2,2,sum)),
                                      num_samps),num_samps,N,byrow=TRUE);
PHI <- measure.mat%*%wavelet.basis
# actual samples we see
y <- measure.mat%*%signal
# use fast Laplace algorithm
w_est <- FindSparse(PHI, y)

# estimate signal
signal_est <- wavelet.basis%*%w_est
# Root mean squared error of estimate
error <- sqrt(mean((signal - signal_est)^2))