looking for strictly recurrent and fast moving median implementation

I am looking for any suitable trick how to effectively compute moving window (length w) median. I know of course the movmedian function, but I need strictly recurrent native MATLAB function working sample by sample.

My naive solution, which is equivalent to the

output_median = movmedian(input_x,[w-1,0]) 

is as follows:

rng('default')
% number of samples
N = 25;
% moving median windows length
w = 5;
% init history buffer and median
x_hist = rand;
med_new = x_hist;
% init input x vector
input_x = zeros(1,N);
input_x(1) = x_hist;
% init output median vector of length N
output_median = zeros(1,N);
output_median(1) = med_new;
for i = 2:N
    x_new = rand;
    [med_new,x_hist] = moving_median(x_hist,x_new,w);
    input_x(i) = x_new;
    output_median(i) = med_new;
end

where function moving_median is here:

function [med_new,x_hist] = moving_median(x_hist,x_new,w)
% Old length of history
wo = length(x_hist);
% Update history
x_hist = [x_hist(max(1,wo-w+2):wo),x_new];     % Grow history until size w, then append new x and remove oldest x
med_new = median(x_hist);
end

Any idea how to make this algorithm more effective (faster) a still strictly recurrent?

Target use case should works with window length:

w ~ 1e3 - 1e4 (!!!)

Additional notes:

1. fast moving median computing (including MATLAB movmedian function) is always based on advanced data structures use like Heap or Queues, etc.
2. some sort-structure information could be stored in these structures and used at the next sample step to significant speed-up of median computing
3. similar approach is used in movmedian function, but this function is not directly applicable on running-data stream!

[Matt J]

I need strictly recurrent native MATLAB function working sample by sample.

I for one don't understand what that means. What is "strictly recurrent"?

Also, what does the "native" requirement mean? You just want to avoid MEX files?

[Michal]

Strictly recurrent ... means possibility to apply moving_median function on every new sample of the "measured" signal separately (sample by sample). Not to apply on complete input vector of samples. See, for-loop in my toy problem.

Native ... means primarily pure MATLAB codes, but proper MEX file is acceptable and maybe a better... :)

It is strange, that for this (very common) type of on-line processing task is not available any proper MATLAB function.

[Image Analyst]

I still don't understand. Sounds like recurrent is just a regular mdeian filter using a sliding window. Just pass each individual sample vector into movmedian. Not sure why you think that does not do the job for you.

How long does your processing take and how fast do you need it to be? There are other denoising methods other than the median filter, which is non-linear and slower than linear filters.

Answer 1

[Steven Lord]

If I understand what you're trying to do, you have a process that generates data (like streaming a video, where you don't necessarily have to download the whole video before you start playing) and as each new data point is entered, you want to take the moving median including that new data point. Is that correct? If so, the simplest solution is probably to call median on the subset of the vector containing your stream of data that is in the window. As a quick and dirty example:

x = (1:5).^2;
m = median(x); % moving median: 2 elements before and 2 elements after
for k = 6:10
    x(k) = k.^2; % If I were making this a full example I'd try to preallocate x and m
    m(k-4) = median(x(k-4:k));
end
m

m = 1×6

9 16 25 36 49 64

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Compare with taking the moving median of the entire x vector, discarding any windows that are not contained entirely inside the vector. If you want other endpoint handling, you should be able to adjust the for loop or the initial creation of the m vector.

m2 = movmedian(x, [2 2], Endpoints = "discard")

m2 = 1×6

9 16 25 36 49 64

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Michal]

The main problem is, that my moving window for median computing is pretty long (w ~ 1e3 - 1e4 samples). So, the direct median computing at each for-loop step is ineffective or simply too slow.

[Michal]

"If I understand what you're trying to do, you have a process that generates data (like streaming a video, where you don't necessarily have to download the whole video before you start playing) and as each new data point is entered, you want to take the moving median including that new data point. Is that correct?"

Yes, I have an on-line stream of 1D-signal data, with sample frequency 1Hz, and I need to compute moving-median with window lenghth w (typically ~ 1e3 samples) for each new sample I get by this stream.

Each new window contains one new (current) sample and one (oldest) sample is removed. The entire median sort should therefore be very similar to the previous one. And this similarity can be effectively exploited to achieve significant speedup of each step of the median calculation.

[Les Beckham]

If your sample rate is only 1 Hz you should have no issue with direct computation of the median of the most recent 1000 or 10000 samples being "simply too slow".

For example:

x = rand([1 1e4]);
tic;median(x);toc

Elapsed time is 0.001857 seconds.

tic;median(x);toc % repeat since the first time will be slower due to JIT processing

Elapsed time is 0.000275 seconds.

Less than half a millisecond doesn't look terribly slow to me.

[Michal]

OK :) you are right, because I tried to make my question simple enough ... but in reality, I have about 500-1000 1D streams and I need to a apply above-mentioned moving-median filtering on each stream separately!

But anyway ... so far is direct median computation on each step viable way to implement in MATLAB some functional demostration prototype ...

[Steven Lord]

Could you combine the 500-1000 individual 1-D streams into a 2-D matrix and call median with a dimension input?

A = randi([-10 10], 4, 5)

A = 4×5

-7 -1 -8 0 -4 6 7 9 -5 0 -2 1 -8 5 1 2 0 -10 -7 -6

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

median(A)

ans = 1×5

0 0.5000 -8.0000 -2.5000 -2.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

median(A, 1)

ans = 1×5

0 0.5000 -8.0000 -2.5000 -2.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

median(A, 2)

ans = 4×1

-4 6 1 -6

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 2

[Matt J]

It sounds like you should be using Simulink, in particular a median filter block,

https://www.mathworks.com/help/dsp/ref/medianfilter.html

[Michal]

This is definitely what I do not looking for... :(

[Matt J]

But you expressed surprise that Matlab doesn't have tools for handling instreaming data. Simulink was made for precisely that.

[Michal]

I do not using Simulink... That is all.

Are you sure that this Simulink block using something better then direct median computation?

[Matt J]

I've never used it, so you'll just have to try, but it is the kind of thing Simulink was buit for.