Recreation of Hampel Method using "medfilt1" function

Please see the attached filles and run hampelTest.m 

I am trying to reproduce the results of running hampel(yv,3) using medfilt1 and the MAD calculation. The value of the median from hampel [shown as xm] is identical to medfilt1(yv,7). This is confirmed in medTest. The value of the std deviation (MAD), shown as xs, is different than medfilt1(abs(yv-medfilt1(yv,7),7). This is confirmed as stdTest. 

Please inform me as to how to reproduce xs from hampel using medfilt1. 

Answer 1

[MathWorks Support Team]

Here I must inform you that the results in "stdTest" is not as you have expected because you did not correctly implement the algorithm to compute the median absolute deviations (MAD) on the data set.  For the i_th MAD value, it is computed as the median of the abs value of 7 data points each subtracted by a1(i).  Instead, you have calculated a2 = abs(yv - a1), which results in the i_th MAD value calculated from the single data point at i, i.e. "yv(i) - a1(i)".  This is not the same as the MAD definition found for the "hampel" calculation

https://www.mathworks.com/help/signal/ref/hampel.html#Hampel Identifier

Since each MAD value depends on a set of "yv" values, consider using a _for-loop_ to compute the standard deviation at each point.  Below, I am providing an example of calculating the standard deviation via direct definition using the "median" function:

 

a3 = NaN(length(yv), 1);
for k = 1:length(yv)
   lowIndx = max(1, k - 3);
   highIndx = min(length(yv), k + 3);
   a3(k) = median(abs(yv(lowIndx:highIndx) - a1(k)));
end

a3 = a3*1.4826; %scale factor as per formula of MAD

stdTest=(xs-a3);

https://www.mathworks.com/help/matlab/ref/median.html

Using the "medfilt1" function is slightly more involved as you must take care of the edge cases manually (i.e. the order of the filter changes near the boundaries of "yv").  However, the idea is similar to that of using the "median" function as provided in the example.