additional findchangepts function output

10 views (last 30 days)
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?

Accepted Answer

Greg Dionne
Greg Dionne on 23 Nov 2016
#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 Comments
Michal
Michal on 24 Nov 2016
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

Sign in to comment.

More Answers (2)

Greg Dionne
Greg Dionne on 23 Nov 2016
#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 Comment
Michal
Michal on 24 Nov 2016
Edited: Michal on 24 Nov 2016
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?

Sign in to comment.


Greg Dionne
Greg Dionne on 28 Nov 2016
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 Comments
Michal
Michal on 19 Dec 2016
I think the problem is harder than appears before, so any viable solution is not available. Am I right?
Greg Dionne
Greg Dionne on 18 Jan 2017
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!