2.2 Digital filtering

The typical output of a long-term numerical integration of an orbit consists of the time series of orbital elements (and possibly other quantities) sampled at regular intervals, much longer than the time step of the integration itself. This decimation of the output is in our case combined with the digital filtering to cut off the short-period perturbing terms. The N-body problem is strongly degenerate, that is the spectrum of the orbital elements is divided in two distinct regions, at high and low frequencies. This allows filtering to be performed independently from the Fourier analysis of the secular system, which by definition contains only the low frequencies. This requires to design filters with performances specifically adapted to the problem; some of the requirements are conflicting, and in practice a suitable compromise must be reached for an optimum performance. The digital filtering procedure applied in the present paper is described in detail in Carpino et al. ([1987]). Thus we only briefly outline the algorithm and the particular choices of the main parameters. The output of the numerical integration can be represented as a time series of orbital elements sampled at regular intervals of length T_s; for a generic element x

$$x_j=x(t_j)\quad(t_j=t_0+jT_s,\ j=0,1,2,\ldots)\,.$$

In order to produce a correct spectral representation of the time series, it is necessary to take into account that the discrete Fourier transform of the signal

$$\hat x(\omega)=\sum_{k=-\infty}^{+\infty} x_k\exp(-ik\omega\,T_s) \quad(i=\sqrt{-1})$$

is periodic with period equal to the sampling frequency $\omega_s=2\pi/T_s$

$$\hat x(\omega)=\hat x(\omega\pm k\omega_s),\quad k=1,2,\ldots$$

and therefore any spectral line having frequency $\vert\omega\vert>\omega_s/2$ will generate a low frequency alias in the region $[-\omega_s/2,\omega_s/2]$, more precisely at the frequency

$$\omega_{\rm alias}=\omega-k_a\omega_s,\qquad {\rm where\ }k_a={\rm int}\,\left({\omega\over\omega_s}+{1\over 2}\right)\,.$$

Of course, this situation is highly undesirable, because the interpretation of a spectrum in which real features are mixed with aliased lines would be difficult. In order to avoid this, it is necessary to use a sampling periodicity T_s satisfying the so-called Nyquist criterion, according to which the sampling frequency $\omega_s$ must be larger than $2\omega_{\rm max}$, where $\omega_{\rm max}$ is the highest frequency present in the spectrum of the series. In the case of the output of the numerical integration of planetary orbits, this would mean to sample the orbital elements at least several times per revolution of the body having the shortest orbital period, in order to avoid aliasing from the short periodic terms produced by its perturbation. For a long-term integration, the adoption of such a short sampling periodicity would in turn produce a huge amount of output: this is not only demanding from the point of view of the computer resources, but would also make the subsequent steps of analysis of the results very cumbersome; on the other hand, if one is interested only in the long-periodic features of the spectrum, the information about short periodic terms is anyway not necessary. A possible solution of this problem consists in the cancellation of the high frequency components of the signal by appropriate digital filtering techniques; after the application of a suitable low-pass digital filter removing all short-periodic terms, the sampling periodicity of the output can be reduced substantially without generating aliasing. Basically the filtering procedure consists in a linear convolution of the time series of the orbital elements x_j with a suitable finite real succession d_k:

$$y_m=(d*x)_m=\sum_{k=-M}^{+M} d_k\,x_{m-k}\,;$$

in this way the spectrum of the convolved signal $\hat y(\omega)$ is equal to the product of the spectrum of the input signal $\hat x(\omega)$ and the frequency response of the filter $\hat d(\omega)$ (which is the discrete Fourier transform of the filter coefficients d_k); the filter can be designed in such a way that its frequency response is almost zero for frequencies larger than a prescribed cutoff frequency $\omega_{\rm cut}$

$$\vert\hat d(\omega)\vert\simeq 0\qquad ({\rm for\ }\vert\omega\vert>\omega_{\rm cut})$$

For the filtered signal y_j the Nyquist criterion becomes

$$\omega_s={2\pi\over T_s} > 2\omega_{\rm cut}\,;$$

of course, in the actual design of the filter, $\omega_{\rm cut}$ must be selected in such way as to preserve the spectral features of interest (in our case, the region of the spectrum corresponding to secular terms).

  

Figure: Frequency response of the digital filter, as used in the 2 Myr integrations.

For the computation of synthetic proper elements, we have used a filter with decimation 100; e.g., in our 2 Myr integrations, we have used an input sampling frequency of one data point every 2 yr, and the output contains one data point every 200 yr. Long frequency aliases can be produced, as a result of a beat with the input frequency, only if they have periods less than the input Nyquist period of 4 yr; thus, in the outer asteroid belt, the perturbations with a frequency close to the mean motion cannot be responsible for spurious long periodic terms. The frequency response of this filter is shown in Figure 1, in particular the dark band is for periods up to $\simeq 300$ yr with response $\vert\hat d(\omega)\vert<10^{-4}$, and the pass band is for periods above $\simeq 720$ yr with response $0.999 < \vert\hat d(\omega)\vert<1.001$. In this way the short periodic perturbations are removed from the filtered output, while the secular perturbations are not significantly distorted. Note that the response at zero frequency is $\hat d (0)=1$; this implies that the long term average is preserved. Moreover, the coefficients of the filter we are using are symmetric, that is d_m=d_-m. In case the filter is applied to an angle variable, provided the angle time series is transformed into a continuous real function by adding $2\pi$ times the number of revolutions, the filter output has the same long term slope and phase as the original data. This is essential to be able to use the filter output to fit frequencies and phases. For the 10 Myr integrations, to avoid excessive output, we have used the same filter with one input data point every 5 yr.