Abrupt changes in data; 'ischange', means, and 'stairs'

23 views (last 30 days)
I have multiple, large data sets that I need to acertain the means along certain linear regions. The regions are dileneated by relatively abrupt changes, as shown in the figure below that I manually marked up.
Ultimately I want a 7 element array of the means of each of these regions (aka the y-intercepts as shown). Thoughts? Other types of curve fittings that could help?
A sample data set is attached.
The below code not working for me, and I've tried a variety of max number of changes and threshold levels.
[TF,S1] = ischange(a, 'linear','MaxNumChanges',12);
plot(a, '*')
hold on
stairs(S1)

Accepted Answer

Matt J
Matt J on 10 Jul 2020
Edited: Matt J on 10 Jul 2020
Here is a method using the Image Processing Toolbox (treating the signal as a 1D image, in other words).
load('dataTest.mat')
w=5000;
b=movmedian(a,w);
cmax=movmax(b,[w,1]);
cmin=movmin(b,[1,w]);
q=cmax-cmin;
lmap=bwlabel(q<0.01);
result = regionprops('table',lmap,b,'MeanIntensity')
plot(a)
hold on
for i=1:size(result,1)
plot(xlim,result{i,1}*[1,1],'--')
end
hold off

More Answers (1)

Matt J
Matt J on 10 Jul 2020
Edited: Matt J on 10 Jul 2020
Yes, you can use splitapply(@mean,data,G)
with G identifying the regions you wwant grouped together,
  4 Comments
noMathWiz
noMathWiz on 10 Jul 2020
Thanks for your help, but I think this is too much manual work (i.e. it requires me to manually find where the changes in data are, delete the data, then create G vectors of the same length).
Matt J
Matt J on 10 Jul 2020
Edited: Matt J on 10 Jul 2020
If you were to construct G manually, it would indeed be a lot of work, but the philosophy behind splitapply is that you would find some automated way to construct the group labels. Below is a method that makes use of group1s from the File Exchange.
load('dataTest.mat')
w=5000;
b=movmedian(a,w);
cmax=movmax(b,[w,1]);
cmin=movmin(b,[1,w]);
q=cmax-cmin;
lmap=group1s(q<0.01);
result = splitapply(@mean, b, findgroups(lmap))
result(1)=[];
plot(a)
hold on
for i=1:size(result,1)
plot(xlim,result(i)*[1,1],'--')
end
hold off

Sign in to comment.

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!