function [ x  ] = OMP( y , D  , k )

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);