Linear segmentation of noisy data

20 views (last 30 days)
Omar Alahmad
Omar Alahmad on 2 Sep 2020
Edited: Bruno Luong on 3 Sep 2020
Hi,
I would like to divide my noisy data into a set of linear segments.
As seen in the above figure. I can measure the red noisy signal. The signal in its nature should be composed of linear segments.
Is there a robust way to perform this segmentation? I.e. convert the red noisy data into the black lined segments.
Thanks!
Regards,
Omar

Accepted Answer

Star Strider
Star Strider on 2 Sep 2020
This would be easier with your data. Lacking them, I created my own.
Try this:
x = 0:32;
y = [-0.5*(0:11)+10+randi([-1 1], 1, 12), 4.5*(12:14)-50+randi([-1 1], 1, 3), 0*(15:28)+13+randi([-1 1], 1, 14), -4*(29:32)+129+randi([-1 1], 1, 4)];
Lv = ischange(y, 'linear', 'Threshold',5);
Idx = [1 find(Lv), numel(y)];
for k = 1:numel(Idx)-1
P(:,k) = [x(Idx(k)) 1; x(Idx(k+1)) 1] \ y(Idx(k:k+1)).';
end
Slopes = P(1,:)
Intercepts = P(2,:)
figure
plot(x, y, '-o')
hold on
plot(x(Idx), y(Idx), 'p')
hold off
grid
axis([0 35 0 15])
In this run, that produced this plot:
and these data:
Slopes =
-0.6667 3.6667 0 -3.0000
Intercepts =
11.0000 -41.0000 14.0000 98.0000
You will likely need to adjust the 'Threshold' value (here 5) to work with your data. It should produce the appropriate slopes and intercepts for the regression lines. The ischange function was introduced in R2017b. A similar function, findchangepts in the Signal Processing Toolbox (with significantly different arguments and outputs) was introduced in R2016a.
  5 Comments

Sign in to comment.

More Answers (1)

Bruno Luong
Bruno Luong on 3 Sep 2020
Edited: Bruno Luong on 3 Sep 2020
Similar topic discussed here
Using my solution of BSFK
% Generate random data
N = 5; % number of linear segments + 1
breaks = cumsum([0, 1+rand(1,N)]);
yb = rand(size(breaks));
coefs = zeros(N,2);
for k=1:N
coefs(k,:)= polyfit(breaks([k,k+1])-breaks(k),yb([k,k+1]),1);
end
pp = struct('form', 'pp',...
'breaks', breaks, ...
'pieces', N, ...
'coefs', coefs, ...
'order', 2,...
'dim', 1);
n = 1000; % number of data point
sigma = 0.01; % Gaussian noise std
x = linspace(min(breaks),max(breaks),n);
y = ppval(pp,x) + 0.05*randn(size(x));
%%
close all
BSFK(x,y,2,[],[],struct('annimation',1)); % FEX

Community Treasure Hunt

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

Start Hunting!