# Efficient implementation of spike removal

I have implemented the moving median absolute deviation (moving MAD) and it seems like bit-exact to Matlab's implementation. Nevertheless, I am sure that it is not efficient.

The usual median filter should be implemented with 2 heaps. The moving MAD can not be implemented this way since the absolute deviation vector for each element is completely different from one sample to another. This forces us to use the quick-select for every sample - which make the calculation very long...

function M = myMovmad(x, xmedian, windSize)
M = zeros(size(xmedian));
for iX=1:length(M)
ind1 = max([iX - (windSize - 1) / 2, 1]);
ind2 = min([length(x), iX + (windSize - 1) / 2]);
M(iX) = median(abs(x(ind1:ind2) - xmedian(iX)));
end
end


Somehow, Matlab managed to implement moving MAD with computation time only twice the median filter. This suggests that they somehow managed to use a double median filter. Any ideas on how to implement it?

Seems like the Hampel filter is part of a more general group of filters called Recursive Median Filters (also related to Robust Scale Estimates). Several filters from this group have a C implementation in the GNU Scientific Library based on the linked article. Yet, their implementation seems to be very similar to the presented above - which is not satisfactory. Is there any better implementation or more efficient spike removal out there?