additional findchangepts function output

12 次查看（过去 30 天）
Michal
Michal 2016-11-23
评论： Greg Dionne 2017-1-18
I have two questions regarding function "findchangepts" (DSP Toolbox):
1. How to effectively create additional output vector "x_hat" from standard output "ipt" of the "findchangepts" function at the same sample grid as input data "x". The "x_hat" vector corresponds to function which is piecewise constant or linear approximation of the input "x" signal with jumps at detected change points "ipt" and is produced only as graphics output by "findchangepts" in case of no output variables (see internal function "cpplot" at "findchangepts.m" source code file).
2. Any idea how to choose the input parameters of the "findchangepts" function to restrict output change points only for jump steps values less than some threshold value?

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

采纳的回答

Greg Dionne
Greg Dionne 2016-11-23
#1 Maybe something like:
function y = fitchangepts(x, icp, statistic)
y = nan(size(x));
K = length(icp);
nseg = K+1;
istart = [1; icp(:)];
istop = [icp(:)-1; length(x)];
if strcmp(statistic,'mean') || strcmp(statistic,'std')
for s=1:nseg
ix = (istart(s):istop(s))';
y(ix) = mean(x(ix));
end
elseif strcmp(statistic,'rms')
for s=1:nseg
ix = (istart(s):istop(s))';
y(ix) = rms(x(ix));
end
else % linear
for s=1:nseg
ix = (istart(s):istop(s))';
y(ix) = polyval(polyfit(ix,x(ix),1),ix);
end
end
Test it:
load engineRPM.mat
plot(fitchangepts(x,findchangepts(x,'Statistic','linear','MinThreshold',var(x)/2),'linear'))
  2 个评论
Michal
Michal 2016-11-24
Thanks Greg, I just made some small additional modifications:
function y = fitchangepts(x, icp, statistic)
% size of input data vector
xs = size(x);
% conditional transposition to column vector (polyfit)
if isrow(x)
x = x';
end
% auxilliary vars
y = zeros(xs);
K = length(icp);
nseg = K+1;
istart = [1; icp(:)];
istop = [icp(:)-1; length(x)];
% find trend for each segment
if strcmp(statistic,'mean') || strcmp(statistic,'std')
for s=1:nseg
ix = (istart(s):istop(s))';
y(ix) = mean(x(ix));
end
elseif strcmp(statistic,'rms')
for s=1:nseg
ix = (istart(s):istop(s))';
y(ix) = rms(x(ix));
end
else % linear
for s=1:nseg
ix = (istart(s):istop(s))';
y(ix) = polyval(polyfit(ix,x(ix),1),ix);
end
end
end

请先登录，再进行评论。

更多回答（2 个）

Greg Dionne
Greg Dionne 2016-11-23
#2 You can get close to this by running FINDCHANGEPTS once with a given threshold, finding all segments that are too long, and re-running on each of these with lower thresholds.
You'll probably want to explain the motivation behind the request though, so the solution works for you.
  1 个评论
Michal
Michal 2016-11-24
编辑：Michal 2016-11-24
I am looking for method which is able to detect jumps (steps) in noised signal. My question is: How to transform the minimum step size threshold to parameter MinThreshold (one of findchangepts function parameters), which represents minimum improvement in total residual error for each changepoint?

请先登录，再进行评论。

Greg Dionne
Greg Dionne 2016-11-28
I don't think I have a good answer to this. The 'mean' option works by performing a sum residual square error, introducing a constant penalty for each break. So if we have a signal with a small shift in mean over a large number of samples, the sum of the residuals would eventually swamp the computation and force a break (no matter how small the shift). This doesn't seem like what you want. The only other option which takes mean into account is the 'std' option; if your noise is distributed uniformly over all segments, maybe that could work(?). No promises of course, but if you share your data I can try to come up with something practical.
  9 个评论
显示 7更早的评论
Michal
Michal 2016-12-19
I think the problem is harder than appears before, so any viable solution is not available. Am I right?
Greg Dionne
Greg Dionne 2017-1-18
The main problem occurs when the slope of the trend is steep and the quantized levels are moving in the opposing sense. Then it becomes difficult to extract.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Measurements and Statistics 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by