3
$\begingroup$

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?

$\endgroup$
4
  • $\begingroup$ While I have written a sliding $\min()$ or sliding $max()$ algorithm that has complexity proportional to $\log_2(B)$ where $B$ is the buffer length, and I understand how to split that into two simultaneous sliding min and max for the top and bottom halves of the data (and the sliding median would be the average of the max of the bottom half buffer and the min of the top half buffer), i have never written the program that is an efficient sliding median. $\endgroup$ – robert bristow-johnson Sep 16 '20 at 0:59
  • 1
    $\begingroup$ @robertbristow-johnson most popular implementation of the sliding-median is the median-heap. I think that there are many implementations that you can find on the web... $\endgroup$ – Gideon Genadi Kogan Sep 17 '20 at 6:50
  • $\begingroup$ You could also review the source code for the median filter implemented in python scipy.ndimage.median_filter $\endgroup$ – Dan Boschen Oct 26 '20 at 14:09
  • $\begingroup$ @Dan Boschen , I was looking for an efficient implementation of the median filter - which I have. I was looking for an efficient implementation of the rolling MAD, which is closely related to the Hampel filter. Unlike the median filter, the Hampel filter modifies the signal in only at the locations which are marked as spike $\endgroup$ – Gideon Genadi Kogan Oct 26 '20 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.