%%
function alpha = Filter(L, s, Niter)
%% 输出参数
% L：行数 * 波段数
% s：波段数 * 1
% Niter: 迭代次数
%%
mu = cell(Niter, 1);
C = cell(Niter, 1);
alpha = cell(Niter, 1);
w = cell(Niter, 1);
%% 初始化

N = length(L); % 行数
mu{1} = mean(L, 'omitnan')'; % 波段数 * 1
C{1} = cov(L, 'omitrows'); % 波段数 * 波段数
alpha{1} = zeros(N, 1);
%% Albedo Factor
tempA = mu{1} .* s; % 8 * 1
tempC = C{1} \ tempA; % 8 * 1
r = (L * mu{1}) / (mu{1}' * mu{1}); % 2056 * 1
tempB = mu{1} - L'; % 8 * 2056
alpha{1} = tempB' * tempC ./ (r * tempA' * tempC);
alpha{1}(alpha{1} < 0) = 0;

%%
IsWeight = false; % 修改该参数
for k = 2 : Niter
    if(IsWeight) % 是否需要权重平衡
        epsilon = 1e-10;
        w{k} = 1 ./ (alpha{k - 1} * 2 + epsilon); % 2056 * 1
    else
        w{k} = zeros(N, 1);
    end
    mu{k} = mean(L - r .* alpha{k - 1} * mu{k - 1}' .* s', 'omitnan')';
    tempA = mu{k} .* s; % 8 * 1
    LC = L - r .* alpha{k - 1} * tempA'; % 2056 * 8
    C{k} = cov(LC, 'omitrows'); % 8 * 8
    tempC = C{k} \ tempA; % 8 * 1
    alpha{k} = zeros(N, 1);
    tempB = mu{k} - L'; % 8 * 2056
    alpha{k} = (tempB' * tempC - w{k}) ./ (r * tempA' * tempC);
    alpha{k}(alpha{k} < 0) = 0;
end
end