Clear Filters
Clear Filters

moving median with variable window

5 views (last 30 days)
Michal
Michal on 23 Sep 2024
Commented: Bruno Luong on 25 Sep 2024
Is there any way how to effectively generalize movmedian function to work with variable window length or local variable k-point median values, where k is vector with the same length as length of input vector (lenght(x) = lenght(k))?
Example:
x = 1:6
k = 2,3,3,5,3,2
M = movmedian_vk(x,k)
M = 1, 2, 3, 4, 5, 5.5
My naive solution looks like:
function M = movmedian_vk(x,k)
if length(k) ~= length(x)
error('Incomaptible input data')
end
M = zeros(size(x));
[uk,~,ck] = unique(k);
for i = 1:length(uk)
M_i = movmedian(x,uk(i));
I_i = (ck == i);
M(I_i) = M_i(I_i);
end
end
  2 Comments
Image Analyst
Image Analyst on 23 Sep 2024
Can you explain the use case? Why do you want to do this?
Michal
Michal on 24 Sep 2024
Edited: Michal on 24 Sep 2024
Robust and effective trend extraction in a case of a priori known 1-D signal parts with high slope changes (typicaly by active control of dynamic system). Median filter then used short window in a case of active control and long windows in opposite case.
But after some additional test I learned that this naive approch is not suitable for reliable trend estimation.
Anyway, I will be very happy for any hint how to apply robust median filter on my use case, where separate parts of signal shoud be filtered with different filter windows (something like weighting).

Sign in to comment.

Answers (3)

Bruno Luong
Bruno Luong on 23 Sep 2024
Edited: Bruno Luong on 23 Sep 2024
One way (for k not very large)
x = 1:6
x = 1×6
1 2 3 4 5 6
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
k = [2,3,4,5,3,2]; % Note: I change k(3) to 4
winmedian(x,k)
ans = 1×6
1.0000 2.0000 2.5000 4.0000 5.0000 5.5000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function mx = winmedian(x,k)
x = reshape(x, 1, []);
k = reshape(k, 1, []);
K = max(k);
p = floor(K/2);
q = K-p;
qm1 = q-1;
r = [x(q:end), nan(1,qm1)];
c = [nan(1,p), x(1:q)];
X = hankel(c,r);
i = (-p:qm1).';
kb = floor(k/2);
kf = k-1-kb;
mask = i < -kb | i > kf;
X(mask) = NaN;
mx = median(X,1,'omitnan');
end
  7 Comments
Matt J
Matt J on 24 Sep 2024
Edited: Matt J on 24 Sep 2024
When there are a small number of unique k(i), yes, yours is best. However, more generally, Bruno's is faster:
k = randi(30,1,1e5);
x = rand(1,1e5);
timeit(@() winmedianMichal(x,k))
ans = 0.0579
timeit(@() winmedianBruno(x,k))
ans = 0.0371
function M = winmedianMichal(x,k)
if length(k) ~= length(x)
error('Incomaptible input data')
end
M = zeros(size(x));
[uk,~,ck] = unique(k);
for i = 1:length(uk)
M_i = movmedian(x,uk(i));
I_i = (ck == i);
M(I_i) = M_i(I_i);
end
end
function mx = winmedianBruno(x,k)
x = reshape(x, 1, []);
k = reshape(k, 1, []);
K = max(k);
p = floor(K/2);
q = K-p;
qm1 = q-1;
r = [x(q:end), nan(1,qm1)];
c = [nan(1,p), x(1:q)];
X = hankel(c,r);
i = (-p:qm1).';
kb = floor(k/2);
kf = k-1-kb;
mask = i < -kb | i > kf;
X(mask) = NaN;
mx = median(X,1,'omitnan');
end

Sign in to comment.


Matt J
Matt J on 23 Sep 2024
Edited: Matt J on 24 Sep 2024
x = rand(1,6)
x = 1×6
0.1034 0.9884 0.6244 0.0233 0.2999 0.6556
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
k = [2,3,3,5,3,2];
n=numel(x);
J=repelem(1:n,k);
I0=1:numel(J);
splitMean=@(vals,G) (accumarray(G(:),vals(:))./accumarray(G(:),ones(numel(vals),1)))';
cc=repelem( round(splitMean( I0,J )) ,k);
zz=min(max(I0-cc+J+1,1),n+2);
vals=[nan,x,nan];
vals=vals(zz);
I=I0-repelem( find(diff([0,J]))-1 ,k);
X=accumarray([I(:),J(:)], vals(:), [max(k),n],[],nan);
M = median(X,1,'omitnan')
M = 1×6
0.1034 0.6244 0.6244 0.6244 0.2999 0.4777
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Matt J
Matt J on 24 Sep 2024
Edited: Matt J on 24 Sep 2024
Anyway, I will be very happy for any hint how to apply robust median filter on my use case, where separate parts of signal shoud be filtered with different filter windows (something like weighting).
If your movmedian windows are simply varying over a small sequence of consecutive intervals, then the code below shows a little bit of speed-up. It won't give the exact same output near the break points between intervals, but it should be fairly close.
x = rand(1,1e5);
k = 8000*ones(1,1e5);
k(20000:30000) =50;
k(18000:20000) =200;
k(30000:32000) =200;
timeit(@() winmedianMichal(x,k))
ans = 0.0134
timeit(@() winmedianMatt(x,k))
ans = 0.0097
function M = winmedianMichal(x,k)
if length(k) ~= length(x)
error('Incomaptible input data')
end
M = zeros(size(x));
[uk,~,ck] = unique(k);
for i = 1:length(uk)
M_i = movmedian(x,uk(i));
I_i = (ck == i);
M(I_i) = M_i(I_i);
end
end
%Requires download of groupConsec
%https://www.mathworks.com/matlabcentral/fileexchange/78008-tools-for-processing-consecutive-repetitions-in-vectors
function M = winmedianMatt(x,k)
M=splitapply(@(a,b){movmedian(a,b(1))}, x,k, groupConsec(k));
M=[M{:}];
end
  5 Comments
Michal
Michal on 25 Sep 2024
@Bruno Luong Could you add the reference on the Weiss paper?
Bruno Luong
Bruno Luong on 25 Sep 2024
Done; somehow this Answers forum and firefox have issue when I edit it, must use another browser.

Sign in to comment.

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!