% Version: 1.0
% Reviser: Huayiwang
% Feature: Denoising Matched Filter
% Date: 2024-4-30

function result = DMF(swir, s)
[row, col, wav] = size(swir); % 尺寸
result = zeros(row, col);
alpha = zeros(row, 1); % 初始化
window_size = 50;
mu_img = zeros(size(swir));
for i = 1 % 迭代次数
    % 估计背景光谱值
    for ii = 1:wav
        img = swir(:,:,ii);
        img_filtered = imfilter(img, fspecial('average',window_size), 'replicate');
        mu_img(:,:,ii) = img_filtered;
    end
    % 按列反演
    for j = 1 : col % 按列进行处理
        L = swir(:, j, :); % 取出一列元素
        L = squeeze(L); % 去除列向量层级，变为行*波段
        mu0 = mu_img(:, j, :);
        mu0 = squeeze(mu0); % 行*波段
        alpha = Filter(L, s, alpha, mu0);
        result(:,j) = alpha;
    end
end

end
%%
function alpha = Filter(L, s, alpha0, mu0)
%% 输出参数
% L：行数 * 波段数
% s：波段数 * 1
wav = length(s);
row = length(L);
%% 初始化
d0 = L - repmat(alpha0, [1 wav]) .* mu0 .* repmat(s', [row, 1]);
mu0 = trimmean(d0, 10)';
tempA = mu0 .* s; % 8 * 1
C = cov(d0); % 8 * 8
tempC = pinv(C) * tempA; % 8 * 1
tempB = mu0 - L'; % 8 * 2056
alpha = (tempB' * tempC) ./ (tempA' * tempC);
alpha(alpha < 0) = 0;
end