The XRF data from our sediment cores from the Chew Bahir show dramatic jumps. If we analyze them statistically, the results will be affected by these changes in the mean. That’s why I was looking for ways to remove them. Here are two of these ways to remove jumps.

First, let us generate synthetic data, a simple sine wave without noise, and add two jumps up and down.

clear, clc, close all t = 0 : 0.1 : 10*pi; x1 = sin(t); x1(1,131:315) = x1(1,131:315) + 5; x1(1,170:315) = x1(1,170:315) - 10;

The first way of removing jumps is to calculate the 1st derivative with diff, with jumps indicated by outliers to be removed e.g. by using the function filloutliers, introduced in R2017a. Then, we integrate the time series and add one single NaN at the end of the time series that got lost while using diff.

x2 = diff(x1); x3 = filloutliers(x2,'pchip'); x4 = cumsum(x3); x4(end+1) = NaN; figure('Position',[100 100 800 1000]) axes('Position',[0.1 0.8 0.8 0.15]) plot(t,x1), grid title('Original Data with Jumps') axes('Position',[0.1 0.6 0.8 0.15]) plot(t(1:end-1),x2), grid title('1st Derivative') axes('Position',[0.1 0.4 0.8 0.15]) plot(t(1:end-1),x3), grid title('1st Derivative, Outliers removed') axes('Position',[0.1 0.2 0.8 0.15]) plot(t,x1,t,x4), grid title('Integrated Signal')

This works fine for the noisefree sine wave but fails to remove jumps from noise, created with x1 = randn(size(t));

Therefore, I tried an alternative approach, using the function findchangepts, contained in the Signal Processing Toolbox and introduced in R2016. The function has been discussed to find change points in an earlier post. We use the function, define the number of change points we expect, and remove the mean (alternatively the median) from the segments between the change points. This approach seems to work fine for both sine waves and noise.

clear x2 x3 x4 maxnumchg = 2; [ipoint,residual] = findchangepts(x1,... 'MaxNumChanges',maxnumchg); m(1) = mean(x1(1:ipoint(1))); for i = 1 : length(ipoint)-1 m(i+1) = mean(x1(ipoint(i):ipoint(i+1))); end m(length(ipoint)+1) = mean(x1(ipoint(i+1):end)); x2 = zeros(size(x1)); x2(1:ipoint(1)) = x1(1:ipoint(1)) - m(1); for i = 1 : length(ipoint)-1 x2(ipoint(i):ipoint(i+1)) = ... x1(ipoint(i):ipoint(i+1)) - m(2); end x2(ipoint(i+1):end) = x1(ipoint(i+1):end) - m(end); figure('Position',[100 100 800 1000]) axes('Position',[0.1 0.8 0.8 0.15]) plot(t,x1), grid, title('Original Data with Jumps') axes('Position',[0.1 0.6 0.8 0.15],... 'XLim',[0 35],'Box','On'), hold on line(t,x1) line([0 t(ipoint(1))],[m(1) m(1)]) for i = 1 : length(ipoint)-1 line([t(ipoint(i)) t(ipoint(i+1))],... [m(i+1) m(i+1)]) end line([t(ipoint(length(ipoint))) t(end)],... [m(end) m(end)]) stem(t(ipoint),6*ones(size(t(ipoint))),... 'Marker','none',... 'Color',[0.4667 0.6745 0.1882],... 'BaseValue',-6) title('Original Data with Change Points') axes('Position',[0.1 0.4 0.8 0.15]) plot(t,x1,t,x2) title('Data without Jumps')

This method also works fine, of course, for the sine wave.

You can download the full code of both experiments using this link.

Very interesting results!!!

Could you add final code with findchangepts function?

Many thanks, Michal! I accidentally forgot to include the code, thanks for note! The code is now included. Cheers, Martin

Thanks a lot for the code…!!!

In my case the signal has typically three signal components:

1. trend

2. global amplitude modulation

3. noise

Something like this, for example:

x1 = t + sin(t) + randn(size(t))

I am looking for any methods which will be able to remove jumps in this kind of signals.

Do you have any idea how to modified findchangepts approach in this case? I think the ‘linear’ statistics option should be used now instead of ‘mean’, but I am not sure.

In this case I would determine and remove the trend, possibly filter it for any other disturbancies, then apply the findchangepts approach to the signal. The “linear” option searches for change points in the trend, if that’s what you are interested in? You find lots of little segments with different trends and jumps in the trend between those.

for signal:

t = 0 : 0.1 : 10*pi;

x1 = t + sin(t) + randn(size(t));

x1(1,131:315) = x1(1,131:315) + 5;

x1(1,170:315) = x1(1,170:315) – 10;

just try the command:

findchangepts(x1,’Statistic’,’linear’,’MaxNumChanges’,2)

or

findchangepts(x1,’Statistic’,’linear’,’MinThreshold’,50)

You get nice trend segments

Yes, I know … we have a paper in review using this technique (mean, std, trend) to find change points in paleoclimate records.

OK … but I am not sure how to correctly apply jumps removing in this case.

In your paper is solved jumps removing in a case of model:

x1 = const + trend + modulation + noise;

e.g. for example:

x1 = 300 + t + sin(t) + randn(size(t));

If yes, could send me this paper for info?

Constant term is not solved in your code, too.

In general I am looking for effective and reliable method to eliminate jumps in signal with constant offset + trend + noise + signal

Well, the paper is about using findchangepts, not about removing jumps. We just started with correcting the jumps in our micro-XRF data, which is a nightmare and will take weeks and months to be completed. I’ll keep you posted if we make some progress. I am sure we will come across very similar problems like you in near future.