Results of OMP-based signal recovery.

% The true signal z is m-sparse in the domain of PSI, having coefficients
% g, m of which are nonzero and the remaining (M-m) of which are zero.
%
%  z = PSI * g
%
% The coloumns of PSI are a basis for the domain of sparsity.
%
% The signal z is subsampled by a matrix PHI, where each column of PHI
% measures z, ultimately yielding a collection of measurements y.
%
%  y = PHI' * z
%    = PHI' * PSI * g
%
% In this example, entries of z that are measured are populated into y
% while those not measured are left as zero.
% 
% Given y, PHI, and PSI, we seek to recover g (and effectively z).
% 
%  PHI: L-by-M measurement matrix
%  PSI: L-by-L sparsity basis
%  y: M-by-1 data vector (observations/measurements)
%  m: degree of sparsity of true signal in domain of PSI
% 
%  g: L-by-1 estimate of sparsity coefficients
%  H: L-by-1 sparse channel estimate (PSI * g)

clc; clearvars;

% dimensionality of true signal
L = 100;

% let's create a sparsity matrix
A = rand(L,L);
[Q,R] = qr(A); % take orthonormal matrix
PSI = Q; % sparsity basis

% sparsity level
m = 5;

% random sampler again
M = round(2*m^2); % number of samples to take
B = zeros(L,M);
idx_samples = randperm(L,M);
r = idx_samples;
c = 1:M;
B(sub2ind(size(B),r,c)) = 1; % only sample a portion
PHI = B;

% sparsity coefficients
s = zeros(L,1);
s(randperm(L,m)) = rand(m,1); % only m are nonzero

% true signal
z = PSI * s;

% measured signal
y = PHI' * z;

% estimate sparsity coefficients and resulting signal
[g,H] = omp_based_cs(PHI,PSI,y,m);

% let's see how we did
figure(1);
subplot(3,1,1);
plot(z,'-r'); hold on;
stem(idx_samples,y,'-b'); hold off; grid on;
xlabel('Sample');
ylabel('Amplitude');
legend('Truth','Observed','Location','Best');
title('Signal');
subplot(3,1,2);
stem(s,'ro'); hold on;
stem(g,'bx'); hold off; grid on;
xlabel('Sample');
ylabel('Amplitude');
legend('Truth','Predicted','Location','Best');
title('Sparse Coefficients');
subplot(3,1,3);
stem(z,'ro'); hold on;
stem(H,'bx'); hold off; grid on;
xlabel('Sample');
ylabel('Amplitude');
legend('Truth','Reconstructed','Location','Best');
title('Signal Reconstruction');
function [g,h] = omp_based_cs(PHI,PSI,y,m)
% OMP_BASED_CS Compressed sensing estimation using OMP.
%
% We observe the true signal (PSI * g) via measurements made with PHI.
%
%  y = PHI' * PSI * g 
% 
% Args:
%  PHI: L-by-M measurement matrix
%  PSI: L-by-L sparsity basis 
%  y: M-by-1 data vector (observations)
%  m: sparsity level of truth
% 
% Returns:
%  g: L-by-1 estimate of sparsity coefficients
%  h: L-by-1 sparse channel estimate (PSI * g)
%
% Reference:
%  "Signal Recovery From Random Measurements Via Orthogonal Matching 
%  Pursuit" by Joel A. Tropp and Anna C. Gilbert.

r = y; % init residual
J = []; % init index set
B = []; % init atoms

A = PHI' * PSI; 

% for each sparsity degree
for i = 1:m
    idx = argmax(A' * r); % greedy selection
    J = [J idx]; % update index set
    B = [B A(:,idx)]; % update our atoms
    x = pinv(B) * y; % projection weights
    a = B * x; % project
    r = y - a; % calculate new residual
end

% sparse coefficients
[~,L] = size(PSI);
g = zeros(L,1);
g(J,1) = x;

% signal reconstruction
h = PSI * g; 

return