Difference between movmad and med(abs(x-med(x)))

Hello everyone,

I'm trying to implement my own version of the hampel filter and I'm having problems with the movmad (Moving Median Absolute Deviation) calculation, used in the original to define the threshold with the variable xmad.

Original hampel.m code (lines 74-84)

% compute the median absolute deviations and the corresponding medians over
% the size of the filter:  ignore samples that contain NaN and truncate
% at the borders of the input
[xmad,xmedian] = movmadmed(x,2*k+1,1,'omitnan','central');

% scale the MAD by ~1.4826 as an estimate of its standard deviation
scale = -1 /(sqrt(2)*erfcinv(3/2));
xsigma = scale*xmad;

% identify points that are either NaN or beyond the desired threshold
xi = ~(abs(x-xmedian) <= nsigma*xsigma);

I have tested that xmad is equal to:

xmad=movmad(x,2*k+1);

As far as I understand, the documentation of movmad says that, in the below code, xmad and my_xmad should be equal:

size=2*k+1;

xmad=movmad(x,size);
my_xmad=medfilt1(abs(x-medfilt1(x,size)),size);

However, as you can see, I obtain different results:

Does anyone know what's the difference between movmad and the operations to obtain my_xmad?

Thanks in advance!

Answer 1

[Cris LaPierre]

You are comparing 2 different equations. The documentation says the equation used in movmad is

You are using medfilt1, which is different than median (applies a third-order one-dimensional median filter). Test your comparison using median and see if it is closer.

[Alvaro Fernandez]

Using median I still have differences. I tested this two approaches:

    function [output] = window_median(input,dim)
    
    rows=size(input,1);
    cols=size(input,2);
    output=zeros(rows,cols);
    
    for f=dim+1:rows-dim
        for c=1:cols
            output(f,c)=median(input(f-dim:f+dim,c));
        end
    end
    
    end
    %--------------------------------------

size=2*k+1;
xmad=movmad(x,size);
%movmad with custom moving median
my_xmad=window_median(abs(x-window_median(x,k)),k);

%movmad using the median of the whole columns (ignoring 0s)
                                           x_nan=x; x_nan(x==0)=NaN;
my_xmad_nozeros=window_median(abs( x-median(x_nan,'omitnan').*(x>0) ),k); 

--------

I used medfilt1 because movmad calculates the median inside a moving window. When medfilt1 is applied to an image, you obtain, for each pixel, the median of the values inside the specified window.

When you use median, it is applied to the whole array, therefore with a matrix you obtain a vector and with a vector you get an unique value.

Is there something I'm getting wrong?

[Cris LaPierre]

You've specified a window, so your implementation of movmad would be more like this

A = [1 2 4 -1 -2 -3 -1 3 2 1];
for i = 1:length(A)
    s = max(1,i-1);
    e = min(length(A),i+1);
    myM(i) = median(abs(A(s:e)-median(A(s:e))));
end

myM = 1×10    
    0.5000    1.0000    2.0000    1.0000    1.0000    1.0000    2.0000    1.0000    1.0000    0.5000

Compare those results to the ones your get from movmad

M = movmad(A,3)

M = 1×10    
    0.5000    1.0000    2.0000    1.0000    1.0000    1.0000    2.0000    1.0000    1.0000    0.5000

[Alvaro Fernandez]

Yes! Your implementation is right. I misunderstood the movmad method.

Thank you very much!