How to find period/envelope that is larger than sample size?

Question

Mercury's distance from the solar system's barycenter (in kilometers) can be approximated as:

5.9115960588705115e7 + 1.1573608483842954e7*Cos(2.086944367456105 - 0.0029751428085695573*t) 

based on 8 years worth of hourly "samples" from HORIZONS. (t is in hours; I forget what value corresponds to t=0, but it should be irrelevant to my question).

The residuals after this approximation look like this:

I could add additional Fourier coefficients, but there appears to be a period or envelope that is larger than the sample size.

Standard Fourier analysis will never find this period, and continuous Fourier analysis is fairly inefficient.

How can I find this period/envelope?

Perhaps more to the point: what's the simplest function that does a "fairly good" job of approximating this data?

The graph shows a pattern, so there "must be" such a function?

Answer 1

If you can guess the number of Fourier coefficients, then a decomposition method such as MUSIC or ESPIRIT might provide a parametric estimate or eigenvector decomposition, given that guess.

Answer 2

Just like you said, this is an envelope, so decomposing the equation to a sum of Fourier coefficients would require a VERY high order, which is not what you want.

You want to get past the "underlying" signal and analyse the envelope, which can be done in several ways like:

1. Low-Pass filter the squared signal to leave the envelope which is of low frequency,then FFT the square root of the low-passed signal.
2. FFT the analytic signal. Use Hilbert Discrete transform to get the analytic signal and FFT its absolute value.

Answer 3

The envelope appears to be just ,,beating'', that is, the result of frequencies almost, but not exactly, multiple of each other. If you can measure each component finely enough, this will also give you the correct ,,beating''. The question would be how to find all components.

In case of Mercury, there will probably be several components. I would try windowing the available data and looking for peaks in the Frequency domain. I wouldn't be surprised, were they the inverse of the difference of Mercuries orbital period and those of other planets. Of course, a window will smear your peaks, and you could have bad luck and not separate two components which actually are distinct.

Answer 4

Here's what I would do to figure out the component cycles of the envelopes of your residuals:

1. find local extrema.
   • construct the upper envelope signal from maxima,
   • construct the lower envelope signal from minima.
2. 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:

1. take the FFT of the envelope,
2. 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.