moving median with variable window

4 次查看（过去 30 天）
Michal
Michal 2024-9-23
评论： Bruno Luong 2024-9-25
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 个评论
Image Analyst
Image Analyst 2024-9-23
Can you explain the use case? Why do you want to do this?
Michal
Michal 2024-9-24
编辑：Michal 2024-9-24
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).

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Bruno Luong
Bruno Luong 2024-9-23
编辑：Bruno Luong 2024-9-23
Ran in:
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 个评论
显示 5更早的评论
Matt J
Matt J 2024-9-24
编辑：Matt J 2024-9-24
Ran in:
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

请先登录，再进行评论。

Matt J
Matt J 2024-9-23
编辑：Matt J 2024-9-24
Ran in:
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 2024-9-24
编辑：Matt J 2024-9-24
Ran in:
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 个评论
显示 3更早的评论
Michal
Michal 2024-9-25
@Bruno Luong Could you add the reference on the Weiss paper?
Bruno Luong
Bruno Luong 2024-9-25
Done; somehow this Answers forum and firefox have issue when I edit it, must use another browser.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Digital Filtering 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by