Here's what I would do to figure out the component cycles of the envelopes of your residuals:
- find local extrema.
- construct the upper envelope signal from maxima,
- construct the lower envelope signal from minima.
- perform cosine curve fitting.
In Python step 1. can be easily done with scipy.signal.argrelmax(), scipy.signal.argrelmin() functions which return the indexes of local extrema found in the input signal. If the extrema in your residual series are fairly equidistant you can use the resulting slices for step 2., otherwise interpolate the minima and maxima over the original time index to get the approximate curves of the envelopes.
Step 2. is more complex. In Python I would use the method proposed by @hwlau here, that is:
- take the FFT of the envelope,
- use the frequency and phase of the strongest component and use it as
the starting point for the least-squares cosine fitting.
I acknowledge the concern raised that FFT does not find the cycles that have periods longer than the length of the input signal, but the curve fitting algorithm can take care of that. In Python for example scipy.optimize.curve_fit() performs up to 1000 fitting attempts and is capable of finding larger cycles in the input data even though it uses FFT results as the starting point. You can see such cases in @hwlau's answer 1.
If you expect the envelopes to be sums of multiple cosines, you can perform the curve fitting a number of times each time subtracting the previously found waves from the input. Here's some code to demonstrate this solution:
import pandas as pd
import numpy as np
from scipy.signal import argrelmax, argrelmin
from scipy.optimize import curve_fit
def get_envelope(s, type_='high'): # or low
e = pd.Series(np.zeros(len(s))*np.nan, index=s.index)
exi = argrelmax(s.values) if type_ == 'high' else argrelmin(s.values)
e.iloc[exi] = s.iloc[exi]
return e.interpolate().iloc[exi[0][0]:exi[0][-1]]
def cosine(t, A, w, ph, c):
return A*np.sin(w*t + ph) + c
def cosfit(s):
F = np.fft.rfft(s.values)
ff = np.fft.rfftfreq(len(s))
guess = [
s.std()*np.sqrt(2), # amp
2*np.pi*ff[np.argmax(F[1:]+1)], # angular frequency
np.angle(F[np.argmax(F[1:]+1)]),# phase
s.mean() # offset
]
return curve_fit(cosine, range(len(s)), s.values, p0=guess) # returns popt and pcov, see scipy docs for info
# Input: pd.Series of residuals
# Step 1.
ehi = get_envelope(res, 'high')
elo = get_envelope(res, 'low')
# Step 2.
N = 8
fithi = ehi.apply(lambda x: 0)
fitlo = elo.apply(lambda x: 0)
for n in range(N):
popt, pcov = cosfit(res[ehi.index].sub(fithi))
fithi = fithi.add([cosine(t, popt[0], popt[1], popt[2], popt[3]) for t in range(len(ehi))])
popt, pcov = cosfit(res[elo.index].sub(fitlo))
fitlo = fitlo.add([cosine(t, popt[0], popt[1], popt[2], popt[3]) for t in range(len(elo))])
# examine the results, ie. plot fithi and fitlo etc.
Notes:
- N is the arbitrary number of fitting operations to perform, could be derived empirically.
- It could be a good idea to save popt from each iteration in a data
frame or a 2D array to have the parameters if the cosine components
of the fitted curves.