function [esti_ang_deg, P_OMP, test_ang_deg] = df_OMP(snapshot, fov_deg, unit_ang_deg, array_struct, sparsity_val, plot_flag)

test_ang_deg = -fov_deg : unit_ang_deg : fov_deg - unit_ang_deg;

Nant = length(array_struct.eff_ch_loc);

% Spectrum
dic_mat = zeros(Nant, length(test_ang_deg));
for ang_idx = 1 : length(test_ang_deg)
    sv = exp(-1i * 2 * pi / array_struct.lambda_c * array_struct.eff_ch_loc.' * array_struct.unit * sind(test_ang_deg(ang_idx)));    
    dic_mat(:, ang_idx) = sv;
end

mc_val = cal_mutual_coherence(dic_mat, 0)

[N, K] = size(dic_mat); % N:dim of signal, K: # atoms in dictionary
if (N ~= size(snapshot))
    error('Dimension not matched');
end

%% Initializing
x = zeros(K,1);      % coefficient (output)
r = snapshot;        % residual of b
omega = zeros(sparsity_val,1);  % selected support
A_omega = [];        % corresponding columns of A
cnt = 0;

%% Iteration
while (cnt < sparsity_val)  % choose sparsity_val atoms
    cnt = cnt+1;
    x_tmp = zeros(K,1);
    inds = setdiff([1:K],omega); % iterate all columns except for the chosen ones
    for idx = inds
        x_tmp(idx) = dic_mat(:,idx)' * r / norm(dic_mat(:,idx)); % sol of min ||a'x-b||
    end
    [~, ichosen] = max(abs(x_tmp)); % choose the maximum
    omega(cnt) = ichosen;
    A_omega = [A_omega dic_mat(:,ichosen)];
    x_ls = A_omega \ snapshot;  % Aomega * x_ls = b
    r = snapshot - A_omega * x_ls; % update r
end

for idx = 1 : sparsity_val
    x(omega(idx)) = x_ls(idx); %x_sparse(i).value;
end

P_OMP = zeros(1, K);
P_OMP(omega) = abs(x_ls).^2;

% Plot
if plot_flag == 1
    figure()
    plot(test_ang_deg, (P_OMP/max(P_OMP)));
    grid on
    xlabel('Angle [deg]')
    ylabel('Normalized Magnitude [linear]')
end

esti_ang_deg = test_ang_deg(omega);

end