function [ x  ] = OMP( y , D  , k )
% Author: Clément DORFFER
% Date: 26/11/2018
% @: clement.dorffer@ensta-bretagne.fr
% Goal: perform a sparse approximation of a signal y using a weighted Orthogonal Matching Pursuit (OMP) technique 
% which was used in:
% 
% C. Dorffer, M. Puigt, G. Delmaire, G. Roussel, Informed Nonnegative Matrix Factorization Methods for 
% Mobile Sensor Network Calibration, IEEE Transactions on Signal and Information Processing over Networks, 
% Volume 4, Issue 4, pp. 667-682, December 2018 ;
%
% and is inspired by the original OMP technique proposed in:
%
% Y. Pati, R. Rezaiifar, P. Krishnaprasad, Orthogonal Matching Pursuit: recursive function approximation 
% with application to wavelet decomposition, in Proc. IEEE ACSSC'93, 1993.
%

residual = y; 

Atoms_sel = zeros(k,1);
D_select = zeros(size(D,1),k);
x_OMP = zeros(k,1);

for iter = 1 : k
    
   [ ~ , idx_max ] = max(abs(residual'*D));
    
    Atoms_sel(iter) = idx_max;
    D_select(:,iter) = D(:,Atoms_sel(iter));
    
    for j = 1 : iter-1
        D_select(:,iter) = D_select(:,iter)-real(D_select(:,iter)'*D_select(:,j))*D_select(:,j);
    end
    
    D_select(:,iter) = D_select(:,iter)/norm(D_select(:,iter));
    x_OMP(iter) = real(D_select(:,iter)'*residual);
    residual = residual-D_select(:,iter)*x_OMP(iter);
    
end
x = zeros(size(D,2),1);
x(Atoms_sel) = real(D(:,Atoms_sel)\y);