function x_hat = OMP(y, A, sparsityLevel)
    
    [~, dimSparsity] = size(A); 

    support = [];
    r = y;

    for idx_L = 1:sparsityLevel
        
        gamma = A'*r;

        [~, idx_max] = max(abs(gamma)); 
        support = [support, idx_max];

        As = A(:, support);
        xs = As \ y;

        y_hat = As*xs; 
        r = y - y_hat; 
    end

    x_hat = zeros(dimSparsity, 1); 
    x_hat(support) = xs; 
end
