detrend a signal with break point but without jumps

Question

Yesterday I asked a question and got the answer to detrend my time series which I think is really better than applying a highpass filter. So I read the description of scipy.signal.detrend and there I learned that there is a parameter bp:

bp : array_like of ints, optional
A sequence of break points. If given, an individual linear fit is performed for each part of data between two break points. Break points are specified as indices into data. ...

I tested this feature with this program:

import math
import numpy as np
from scipy.signal import detrend
from matplotlib import pyplot as plt

low  = np.zeros((1500), dtype = float)
high = np.zeros((1500), dtype = float)
for i in range(500):
    p1 = i * math.tau / 1000
    p2 = i * math.tau /  500
    p3 = i * math.tau /   50
    c1 = -math.cos(p1)
    c2 =  math.cos(p2) / 2 + 1/2
    c3 =  math.cos(p1)
    c4 =  math.cos(p3) / 2
    low [i +    0] = c1
    low [i +  500] = c2
    low [i + 1000] = c3
    high[i +    0] = c4
    high[i +  500] = c4
    high[i + 1000] = c4

raw = low + high
detrended = detrend(raw, bp=[500,1000])
trend = raw - detrended

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(nrows = 5, ncols = 1, sharex = True)

ax1.plot(low, linestyle="solid")
ax1.set_title('low frequency data')

ax2.plot(high, linestyle="solid")
ax2.set_title('high frequency data')

ax3.plot(raw, linestyle="solid")
ax3.set_title('to be detrended = low + high')

ax4.plot(detrended, linestyle="solid")
ax4.set_title('detrended data')

ax5.plot(trend, linestyle="solid")
ax5.set_title('trend')

plt.tight_layout()
plt.show()

And I got this diagrams as output:

I put the break points at exactly that two points where the low frequency curves has its maxima.

As you can see, there are jumps in the trend line at the break points and therefore also the detrended signal has jumps that didn't exist in the raw data. This is exactly what I feared and expected, but not what I hoped to get.

I hoped to get something like this:

The trendline still has sharp kinks (which is fine), but it no longer has jumps. Therefore, there are no jumps in the detrended data either.

Is there a simple way to detrend a signal with breakpoint that do not cause jumps?

Answer 1

I'm answering my own question

I found a solution for my detrending problem that fits well to me, and maybe other researchers or developers who stumbled across this question might find this solution helpful for their own purposes. That's why I post it as an answer to my own question. But be aware, that other solutions exist and might fit better in other situations.

---

My solution:

Detrending means to calculate a certain simple trend function (usually a constant or a linear function) from the input data and then to subtract this trend from the original function. The difference is the detrended version of the original data.

But what if a simple linear function seems to be much too simple like in my example? (I have a time series with a few million samples recorded over many hours. Having only one linear function as a trend over the entire record does not seem very likely.) My idea is that the trend is a 'blurred' version of the input data, like a sliding average, calculated from a short section of the input data (a "sliding window") with a fixed length, but moving across the whole input data.

But the hard edges of such a window might cause some unwanted side effects, so I used a spline function with soft edges instead. There are many possible forms that such a spline might have. I have chosen this spline:

import math
import numpy as np
from matplotlib import pyplot as plt

half_splinewidth = 100
spline =  np.zeros((half_splinewidth * 2 + 1), dtype = float)
for n in range(half_splinewidth + 1):
    x = (n / half_splinewidth) ** 2 * math.pi
    w = (math.cos(x) + 1) / 2
    spline[half_splinewidth + n] = w
    spline[half_splinewidth - n] = w

# normalize such that the area under the curve == 1
integral = np.trapz(spline)
spline /= integral

plt.plot(spline)
plt.title('spline')
plt.show()

I applied that spline to the raw data with an algorithm named convolution (wikipedia). Scipy can do this trick (scipy.signal.convolve):

from scipy.signal import convolve

trend = convolve(raw, spline, mode='same')
detrended = raw - trend

The parameter mode influences the behavior at the beginning and at the end of the input data and the length of the output. same means that the output has exactly the same length as the input and that the input is prefixed and appended with so many zeros that the convolution can be calculated.

Having the same length is what I wanted, but assuming leading and trailing zeros is not really optimal. As a result I got this:

The trend differs from the low frequency data at the beginning and the end of the output, but is very close to it in the main section of the data.

---

The behavior at the edges can be improved my manually adding a short section of a mirrored copy of the input data to the beginning and the end:

# Note that len(spline) == half_splinewidth * 2 + 1
padded = np.concatenate((
    np.flip(raw[:splinewidth]), 
    raw, 
    np.flip(raw[-splinewidth:])
))
trend = convolve(padded, spline, mode='valid')

Note, that here I used mode='valid'.

And here is the result: