Detrend Flux Time Series with Non-Linear Trend

Question

I am needing to detrend flux time series data (light curves), but I'm running into a problem when the time series data doesn't have a simple linear trend.

I've been using scipy.signal.detrend() for the detrending of linear cases, but that isn't sufficient here.

I've used numpy.polyfit() to attempt polynomial detrending, but I'm not sure what to do with the polynomial coefficients it returns.

Can someone advise me as to the next intelligent step to take? Or, if someone has a better method for detrending non-linear data, I'd be delighted to hear that as well.

Answer 1

In a nutshell, you take the coefficients that polyfit returns and pass them to polyval to evaluate the polynomial at the observed "x" locations.

As a stand-alone example, let's say we have something similar to the following:

import numpy as np
import matplotlib.pyplot as plt

num = 1000
x = np.linspace(0, 10, num)
y = np.exp(x)

# Add some non-stationary noise that's hard to see without de-trending
noise = 100 * np.exp(0.2 * x) * np.random.normal(0, 1, num)
y += noise

fig, ax = plt.subplots()
ax.plot(x, y, 'ro')
plt.show()

Note that I haven't used a polynomial function here to create y. That's deliberate. Otherwise, we'd get an exact fit and wouldn't need to "play around" with the order of the polynomial.

Now let's try detrending it with a 2nd order polynomial function (note the 2 in the line model = np.polyfit(x, y, 2)):

import numpy as np
import matplotlib.pyplot as plt

num = 1000
x = np.linspace(0, 10, num)
y = np.exp(x)

# Add some non-stationary noise that's hard to see without de-trending
noise = 100 * np.exp(0.2 * x) * np.random.normal(0, 1, num)
y += noise

# Detrend with a 2d order polynomial
model = np.polyfit(x, y, 2)
predicted = np.polyval(model, x)

fig, axes = plt.subplots(nrows=2, sharex=True)
axes[0].plot(x, y, 'ro')
axes[0].plot(x, predicted, 'k-')
axes[0].set(title='Original Data and 2nd Order Polynomial Trend')

axes[1].plot(x, y - predicted, 'ro')
axes[1].set(title='Detrended Residual')

plt.show()

---

Notice that we didn't fit the data exactly. It's an exponential function and we're using a polynomial. However, as we increase the order of the polynomial, we'll fit the function more precisely (at the risk of starting to fit noise):