function [RR_refined, RR_Normal_Idx, RR_Large_Idx, RR_Small_Idx, RR_Left_max, RR_Right_max] = hampel(RR)

X=(1:length(RR))';
DX=3;
T=3;

S0 = NaN(size(RR));
RR0 = S0;

if numel(X) > 1
    for i = 1:numel(RR)
        [RR0(i), S0(i)]  = localwindow(X, RR, DX, i);
    end
end

RR_Normal_Idx = find(abs(RR - RR0) <= T*S0);
RR_Large_Idx = find(RR - RR0 > T*S0);
RR_Small_Idx = find(RR - RR0 < -T*S0);

RR_Left_max = RR0(1) + T*S0(1);
RR_Right_max = RR0(end) + T*S0(end);

RR_refined = RR;
RR_refined(RR_Large_Idx) = RR0(RR_Large_Idx);
RR_refined(RR_Small_Idx) = RR0(RR_Small_Idx);
end

function [RR0, S0] = localwindow(X, RR, DX, i)

    Idx_Local_Win = X(i) - DX <= X & X <= X(i) + DX;
    RR0 = round(median(RR(Idx_Local_Win)));
    S0 = 1.4826*median(abs(RR(Idx_Local_Win) - RR0));
end

