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?

Improve this question
$\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, 2020 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, 2020 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, 2020 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, 2020 at 17:37

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.