additional findchangepts function output

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?

Answer 1

[Greg Dionne]

#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'))

[Michal]

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

[Greg Dionne]

Looks good.

Answer 2

[Greg Dionne]

#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.

[Michal]

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?

Answer 3

[Greg Dionne]

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.

[Michal]

OK ... here is the link with some example real data .

Array temperature contains 210 temperature signals (sampling period 1sec). Values "-2" in signals mean "no signal status".

Let's show results for 160th signal with nice jumps:

findchangepts(temperatures(:,160), 'Statistic','linear','MinThreshold',0.5,'MinDistance',10)

I am looking for jumps with step value greater than MinStep > 0.5 degC, for example. Statistic "linear" is required, because the whole signal is slowly modulated by global trends, so I do not need to detect artificial steps slow trends in data.

Note: for proper functionality of findchangepts use:

temperatures = double(temperatures);

Function has problems with "single" type of input arrays (already reported bug).

[Greg Dionne]

Thanks for the data. If I'm understanding you correctly, it seems you want to recover the quantized steps, removing just the noise component and leaving the global trend intact. If so, do you also have the global trend that distorts the otherwise flat quantization levels? (see temperatures(:,16) for an example of this) If not, it will be considerably harder.

[Michal]

Regarding steps detection only, you are right. But, as you can see, the steps are performed typically between small number of levels modulated by global trend. And the robust estimation of these levels need to remove global trends, which are always present. This is the main reason, why I am using the "linear" statistic, because in this case the function "fitchangepts" provides segment between change points with constant slope. The knowledge of this information, is probably sufficient information for reliable global trends removing. After that is possible to apply stnadard histogram technique to establish number of step levels by significant peaks in histogram.

So finally, do you have any idea how to remove global trends from signals? I think, that "linear" fitchangepts output will be good info for this task.

What is your opinion about that?

[Greg Dionne]

The good news is that your trend is slowly varying, so it has a good chance of not surviving operations with diff(). My first thought would be to perform diff, then somehow quantize the results to units of 0.5, then restore via cumsum(). As to the correct DC offset to use, we could use median() on a segment that is fairly long.... maybe find ways to do this taking into account adjacent samples (e.g. groups-of-three). I'll need to think about this.

[Greg Dionne]

This might get you started(?):

 findchangepts(unwrap(mod(temperatures(1:3000,16),.0125)*(2*pi/0.0125)),'statistic','linear','Minthreshold',22)

Basically compute the signal modulo .0125, back out the slow trend. use findchangepoints to fix where the trend has residual jumps, then subtract it off. Then round to nearest .0125 after all is said and done...

[Michal]

Thanks ... but I am afraid I do not understand well what exactly you suggest me. Could you show me your idea by a simple script, with all steps you mentioned?

[Greg Dionne]

I'll take a look at this over the weekend...

[Michal]

I think the problem is harder than appears before, so any viable solution is not available. Am I right?

[Greg Dionne]

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.