The use of the empirical mode decomposition for the identification of mean field aligned reference frames

The magnetic field satellite data are usually referred to geocentric coordinate reference frame. Conversely, the magnetohydrodynamic waves modes in magnetized plasma depend on the ambient magnetic field, and is then useful to rotate the magnetic field measurements into the mean field aligned (MFA) coordinate system. This reference frame is useful to study the ultra low frequency magnetic field variations along the direction of the mean field and perpendicularly to it. In order to identify the mean magnetic field the classical moving average (MAVG) approach is usually adopted but, under particular conditions, this procedure induces undesired features, such as spectral alteration in the rotated components. We discuss these aspects promoting an alternative and more efficient method for mean field aligned projection, based on the empirical mode decomposition (EMD).


Introduction
Ultra-low frequency (ULF; 1 mHz-10 Hz) magnetohydrodynamic (MHD) plasma waves received particular attention in the past years [Eriksson et al. 2006, Pilipenko et al. 2008, Agapitov and Cheremnykh 2015, Belakhovsky et al. 2016, Balasis et al. 2015]. Generated by a variety of instabilities, ULF waves transport energy throughout the magnetosphere, and may play important roles in the energization and loss of radiation belt particles (see Menk [2011] for a review). In particular, ULF waves provide a convenient probe of the magnetosphere, by means of ground [Lichtenberger et al. 2013] and/or satellites magnetic field measurements [Glassmeier et al. 2001, Clausen et al. 2009.
Satellite data are usually referred to a geocentric reference frame, while MHD waves propagation and properties can be established, using a more convenient reference system. In this work we refer to measurements from fluxgate magnetometers, which provide three components of the magnetic field, for both magnetospheric and interplanetary measurements. We also assume, without losing generality, that the geocentric solar ecliptic (GSE) is the original spacecraft reference frame.
Usually, for satellite mission that continuously explore both the magnetosphere and upstream (foreshock) regions, the mean field aligned (MFA) coordinate system is largely used for the rotation procedure [e.g., Clausen et al. 2009, Sarris et al. 2009, Francia et al. 2012]. In the solar wind, a technique widely used to study parallel and perpendicular magnetic field components of MHD waves is represented by the minimum variance analysis (MVA) [e.g., Tu et al. 1989, Klein et al. 1991, Bruno and Carbone 2013. However, the minimum variance direction does not necessarily coincide with that of the ambient magnetic field. In particular, Bruno et al. [1985] found that the angle between the minimum variance axis and the mean magnetic field direction ranges in the interval 8-30 degrees, for heliocentric distances from 0.29 to 0.88 AU. Other methods compute the average magnetic field direction by means of the moving average (MAVG) procedure, that will be described in Section 2, and then apply the MVA on the projected components into the orthogonal plane with respect to the mean direction [Pulkkinen and Rastätter 2009]. In these cases, the common key point concerns the computation of the main magnetic field direction.
ULF waves property can be associated with several anisotropy conditions. An important issue for many types of anisotropy is the scale over which the mean magnetic field is defined. Time-scales of the ambient magnetic field and of the perturbation are crucial quantities. For example, considering the solar wind measurements, if time scales of the fluctuations are comparable with the scale of the background flow, interactions between the background flow and phase variations of the fluctuations will influence the background flow [Tu et al. 1989]. For a local mean field, it might be the same scale as that of the considered fluctuations, or perhaps 2, 5 or 10 times that [Oughton et al. 2015]. On this regard, the magnetic variance anisotropy (see TenBarge et al. [2012] for a review) of the solar wind, as well as the ratio between the ambient magnetic field and the magnetic field fluctuations [Tu et al. 1989, Bruno andCarbone 2013], can be regarded as useful quantities to identify the nature of solar wind turbulent fluctuations.
If a time scale separation exists within the magnetic field measurements, these time series can be thought as a superposition of a slowly varying (amplitude) ambient field B 0 (t), an higher frequency signal b(t) (or perturbation) and an incoherent noise n(t): The MFA coordinates system, showed in Figure 1 for an assigned position in the inner magnetosphere, is established by means of the unit vectors defined as ( indicates the norm): where µ, z and o are usually associated with compressional, toroidal and poloidal ULF waves modes respectively, while r(t) represents the position vector of the spacecraft [Sarris et al. 2009.
Although this definition is usually used for magnetospheric field measurements, it may also be extended in upstream regions (see for example Francia et al. [2012]), but in this case the z and o components are simply related with transverse oscillations in the interplanetary region (e.g. foreshock upstream waves).
Using the Expressions (2) we define the instantaneous rotation matrix from geocentric to MFA reference frame as that allows us to project the instantaneous magnetic field vector from the original geocentric reference frame into the MFA one: It is clear from Equations (2) and (1) that, in order to obtain the time series B MFA (t) in the MFA reference system, the slowly varying mean field vector B 0 (t) must be found through an appropriate filtering procedure. A suitable filter should be applied to the data to remove longer timescale variability. Although filtering procedure has the benefit to potentially reduce the noise, it should be applied with caution, as it may also introduce artefacts affecting the signals altogether. In particular, inverse fast Fourier (IFFT) based filters should be used only for stationary and linear phenomena. On this regard, in the MFA rotation procedure, the moving average (MAVG) is a good choice with respect to any low-pass filter, and hence it is used to remove noise and higher frequency fluctuations in a given time series (see for example Clausen et al. [2009], Francia et al. [2012, Regi et al. [2013]).
In applying the MAVG method it is assumed that the characteristic fluctuation time T 0 related to B 0 (t) is much greater than the period T b of the perturbation b(t). Moreover, T 0 depends on both satellite motion and natural magnetic field variation (e.g. high velocity stream, coronal mass ejections, corotating interaction regions), and could be related to non linear and non stationary phenomena.
Under these conditions the MAVG might be unsuitable in the rotation procedure, while a method such as the empirical mode decomposition (EMD), is useful to identify non linear and non stationary processes [Huang et al. 1998].
In Section 2 we discuss on the intrinsic problems involved in the MAVG procedure on a discrete time series in the time domain; in Section 3 we introduce a new method to identify B 0 (t) using the EMD, comparing it to the classical ones (Section 4); finally, we discuss on the spectral modification induced by MAVG proce- r q q p q q q q q q q q q q q q q q q q q q q q q r q q q q r  [Tsyganenko 1995] during solar quiet conditions. dure in the MFA coordinate system, showing an experimental example in Section 5.

The MAVG effects
The simplest procedure to identify B 0 (t) is the estimation of the average over a time window T, centered at a given time t. By considering a periodic perturbation, Equation (1) assumes, for each component, the following form where b 0d (d = x, y, z) are the amplitudes, a d are the phases, and ~d = 2r/T b,d are the angular frequencies.
Assuming that T b = max{T b,x , T b,y , T b,z } and T 0 = min{T 0,x , T 0,y , T 0,z }, and selecting a time length T greater than T b and smaller than T 0 since the average of the perturbation in a time interval T ≥ T b is negligible and the average noise rapidly approaches to zero, because its auto correlation time is short compared to that of the physical phenomena involved [Regi et al. 2015[Regi et al. , 2016. We now adapt the procedure to discrete time series, defining the time t i = ih, where h is the sampling period and i = 1,2,…,n, with n the total number of measurements. In order to compute the average field around the i th -sample, we use a time-window where N is the integer number of data on the left and on the right with respect to the i th -central data, and hence N w is the window size (i.e. the total number of data involved in the average). Under this condition the MAVG, in discrete form, for the i th term of the magnetic field, can be written as that represents the procedure for the symmetric MAVG computation.
The main goal of the MAVG is to suppress both 〈b〉 i and 〈n〉 i terms in Equation (11), so that 〈B 0 〉 i ≃ 〈B〉 i , and Equation (2) can be rewritten as follows (here 〈r〉 i represent the averaged position vector of the spacecraft), obtaining the MAVG rotation matrix below However, the assumption 〈b〉 i + 〈n〉 i ≃ 0 due to the MAVG procedure is not always verified, and each term in Equation (13) (generally) depends on the perturbation 〈b〉 i , and lesser on the noise. Figure 2 illustrates qualitatively the MAVG effects on the perturbation and noise, i.e. 〈b〉 i + 〈n〉 i terms in Equation (11 r q q p q q q q q q q q q q q q q q q q q q q q q q q q r q q q q r more sophisticated test on the MAVG effects will be illustrated in Section 5.1. In this test an artificial signal with h =1 s, obtained as the sum of a 7 mHz sine wave of unitary amplitude, that represents the perturbation (along an assigned direction), and a normally distributed white noise with zero mean and variance v 2 ≃ 0.25 (top panel), are considered. We can see that the shape of the MAVG (bottom panel) strongly depends to the choice of the windows size N w , with amplitudes decreasing with increasing N w . In particular, the noise amplitude decreases more rapidly with respect to that of the perturbation, however both quantities do not completely vanishes. A quantitative evaluation of MAVG effects on this kind of signals can be found analytically solving Equation (11). Considering the discrete version of Equation (5), and that the harmonic perturbation at instant t i is represented by the amplitude b i and the angular frequency ~d = 2r/T b,d = 2r/N b,d h, Equation (11) can be written as In this case, we assumed that 〈n〉 i~0 , and the moving-averaged 〈b〉 i is computed by means of the relationship (29) (Appendix A). In Equation (14) b 0d sin (2ri/N b,d +a d ) is essentially the perturbation b d,i , while g is the time-independent (dimensionless) factor (see Equation (29) in Appendix A) resulting in a superposition of cosine functions related only with N w and N b .
In general, for assigned N w , |g| increases with N b , while for assigned N b it decreases with N w . More in detail, local minima (close to but not equal to 0) are obtained for N w = cN b where c = 1, 2, 3, ..., c max , and hence when the MAVG window length is exactly equal to an integer multiple of N b . Here, c max is computed as [N s /N b ], where N s is the total number of samples, and [] is the integer part operator. The choice of N w in the MAVG procedure strongly affects the estimation of the rotation matrix coefficients. In the satellite reference frame, the time fluctuations of the ambient field B 0 (t) can be due to its intrinsic variation, but also to satellite motion. In the magnetosphere, the prevailing time variations are typically due to the satellite motion through the magnetic field lines (except transient features observed for example in the geomagnetic tail) while, in the interplanetary medium, the ambient field variations depend on natural phenomena ranging from coronal mass ejections (CME) and plasma cloud passages (that are more frequently observed during the higher solar activity period), as well as co-rotating interaction regions (CIR) linked with coronal holes (mainly observed during the declining phase of the solar cycles). The prob-lem becomes more complicated if we consider the effect of several perturbations with different N b , embedded in the ambient magnetic field. In this case, the higher frequency components (with lower N b ) are removed from the time series, while the lower frequency ones (or higher periodicity) can affect the procedure. The MAVG effects on rotation matrix are discussed in Section 4.

An EMD based procedure for low-frequency components identification
An alternative procedure to extract the main magnetic field from a signal can be obtained using the empirical mode decomposition (EMD) method [Huang et al. 1998]. The EMD technique decomposes a time series into roughly zero-mean (mutually orthogonal) components called intrinsic mode functions (IMFs), and a residual term (Res). Each IMF is found using an appropriate convergent procedure called sifting algorithm.
The sifting process uses the upper e u and lower e l envelopes computed respectively by means of cubic spline interpolation of the local maxima and minima separately identified (see Huang et al. [1998] for detail). At the first step, the envelopes from the original data B(t) are computed. By defining 〈e〉 1 as the mean of the envelopes, the first extracted component h 1 is the difference between B(t) and 〈e〉 1 [Huang et al. 1998], (i.e. B(t) −〈e〉 1 = h 1 ). This procedure has to be repeated several times, until preassigned stopping criteria are not satisfied (e.g., Huang et al. [1998], Rilling et al. [2003], Flandrin et al. [2004], Wu and Huang [2004]). Let assume that criteria are satisfied at the k th -iteration, the first intrinsic mode function (IMF) B 1 is represented by h k . At this step, the sifting algorithm is applied to B(t) − B 1 , to extract B 2 , and so on. The final mode (hereafter residue, Res) can still be different from zero; for data with a trend, then the final Res should be that trend, characterized at most by a unique zero crossing.
After the EMD decomposition, the original time series can be regarded as superposition of all IMFs and the residue: where B k (k = 1,2,...,N m ) are the IMFs.
In particular, the EMD decomposition can be regarded as a dyadic filter bank [Flandrin et al. 2004], and hence it is expected that the number of extracted IMFs (or modes) N m on a time series represented by N s samples, should be N m < log 2 (N s ) or, at least, of the order of log 2 (N s ). If the inequality condition just above is not satisfied, the studied time se- ries possesses a more complex structure and a higher information content than one of a purely stochastic noise [De Michelis et al. 2012].
In principle, sifting algorithm extracts IMFs with energies and time-scales (or periodicities) that should be compared with that obtained with pure randomized time series (see Huang [2009] andDe Michelis et al. [2012]). In order to retain physical meaning in the extractedmodes, different sifting algorithm, with different stopping criteria are proposed by several authors (e.g., Huang et al. [1998], Rilling et al. [2003], De Michelis et al. [2012]).
In this work we used the EMD procedure described in Rilling et al. [2003] (see also Huang et al. [1998]), that utilizes a stopping criterion for the sifting algorithm based on two thresholds i 1 and i 2 . This criterion ensures globally small fluctuations, but at the same time, takes into account locally slow large excursions. Following the authors we introduce the mode amplitude a(t) = (e(t) max − e(t) min )/2, where e(t) max and e(t) min are the envelopes, and the evaluation function v(t) = |m(t)/a(t)|, where m(t) is the mean between these envelopes. The sifting algorithm is iterated until v(t) < i 1 , for some assigned fraction (1 − a) of the total duration, while v(t) < i 2 for the fraction a of the total duration. Regarding the boundary conditions, the algorithm in Rilling et al. [2003] minimizes the error propagations due to finite observation lengths by mirroring the extrema close to the edges. We performed in a separate analysis, the EMD based rotation algorithm for the real case of upstream waves event observed by Cluster satellite showed in Section 5.2, for several stopping parameters i 1 and i 2 , while for simplicity we chose a constant tolerance of a = 0.05. The results (not shown here) confirmed that our algorithm does not introduce significant spectral modification in the MFA reference frame assuming i 1 close to 0.05 and i 2 close to 0.5, as suggested in Rilling et al. [2003]. Then, according with the authors in Rilling et al. [2003], we chose as default values a = 0.05, i 1 = 0.05 and i 2 = 10i 1 = 0.50. For simplicity, in the following discussion, we will refer to a single component of the magnetic field B i = B 0,i +b i + n i . As explained in the introduction, we hypothesize that the B 0 (t) time scale is higher than that of b(t), and hence we found useful to define some quantity in order to reconstruct B 0 (t) after the EMD procedure. Referring to the perturbation field b(t), we search for the maximum time length T max unaffected by typical perturbation phenomena, defining the equivalent maximum number of samples as M max = T max /h.
Each k th -IMF can be subdivided in a finite number of subintervals, identified by two consecutive zero crossing. Let assume that each subinterval I encloses q I samples. We compute, for each of them, the period T I = 2q I h, and an energy-linked quantity P I = ∑|B i∈I |. By means of the previous quantities, we evaluate: (a) the average periodicity of an IMF; (b) the ratio between the sum of P I of subintervals characterized by a periodicity T I > T max (namely P 0 ), and the total energy of all subintervals (namely P = ∑ P I ) Finally, we reconstruct B 0 (t) EMD simply summing the residual Res of the EMD method and IMFs characterized by 〈T I 〉> T max or R > R min , while the remaining IMFs, with lower periodicities, reconstruct [b(t) + n(t)] EMD . The last condition ensures that the IMFs characterized by a high energy even for a short time interval can be identified as part of the main field. The lower is R min the more the method accepts the IMFs having high mean frequency.
We performed a Monte Carlo test in order to evaluate the best value of R min , generating 10,000 surrogates of Gaussian white noise (each one with v n =0.3), while B 0 and b are fixed and are represented by sinusoids with period greater/lower than T max . We computed, for different R min , the correlation coefficients r between B 0 and B 0 EMD . The results showed that our algorithm reconstructed B 0 with r > 0.99 in over 99% of the cases, for R min = 0.1.
An example is showed in Figure 3, using artificial signals with h = 1 s. Panel (a) shows B 0 (t) and b(t) + n(t), while in the panel (b) the resulting B(t) = B 0 (t) + b(t) + n(t) is showed. The perturbation is obtained windowing a 23 mHz sine wave with a 1000 points Hanning window centered at t = 1000 s. We assume in this simulation that B 0 (t) has a period greater than 1000 s. Following our procedure, the EMD method decomposes B(t) into 7 IMFs and a Res (panels c); the average period M of each IMF is also showed. By choosing T max = 1000 s (i.e. M max = 1000), identifying IMF7 and the Res (red lines) as B 0 (t). Using our method, both B 0 (t) EMD and [b(t) + n(t)] EMD are then reconstructed (panel d). It can be noted that the original B 0 (t) is well reconstructed, while b(t) + n(t) can be obtained by superimposing the remaining IMFs. Although we used the EMD method described in Rilling et al. [2003], any kind of EMD procedure can be used following our methods, such as the ensemble empirical mode decomposition (EEMD; Wu and Huang [2009]) or the complete ensemble empirical mode decomposition (CEEMD; Torres et al. [2011]). However, we stress that different sifting algorithms work with different stopping criteria, and hence, the R min should be evaluated by means of a dedicated Monte Carlo test for other EMD algorithms. We conclude this section with some details on the selection criteria examined here. They are based on the observed numerical test results, as well as on the direct application on several magnetospheric and upstream region events. Regarding T max , it should be chosen accordingly with what assumed in Section 1. Indeed, if the periodicity T b is much smaller than T B , a time scale separation exists within the magnetic field measurements, and T max will be chosen within Finally we remark also that, because of the adaptive nature of the basis, EMD is ideally suited for analyzing data from nonstationary and nonlinear processes. However, the reader must bearing in mind that EMD still cannot resolve the most complicated cases, when the processes are nonlinear and the noises also have the same time-scale as the signal. In these cases their separation becomes impossible (see Wu and Huang [2004] for details).

The comparison between EMD and MAVG methods
We compared the EMD and MAVG methods, basing on the difference between the original signal B 0 (t) and that, B 0 th (t), obtained using different methods where, for simplicity, we used the same artificial signal B(t) showed in Figure 3b. Figure 4 shows the original signal B(t) (panel a), and DB 0 (t) computed using the EMD (green) and the MAVG (blue) procedures (panel b). In the MAVG method a window size of 75 points is used. As we can see the MAVG method is affected by the presence of a residual oscillation at the same frequency of the original perturbation, but with a phase shift of 180°, and both method show irregular deviation from the expected signal due to the noise. The results are more clear by looking at the panel (c) (see also the magnified time interval showed in panel d), where we averaged the results for 10,000 signals differing only for the white noise, highlighting then the method dependent effects. In particular, the residual signal (red) arising from the MAVG procedure reproduces the pattern of the original perturbation, but with a phase shift of 180°, while the EMD residual results almost insignificant. The amplitude of the residual oscillation and the phase shift depend on the window size used for the moving-average. We computed the maximum amplitude difference and phase lag between the original B 0 (t) and moving-averaged B MAVG (t) signals, for several values of N w . The maximum difference ( Figure 5, top panel) generally decreases with N w , this trend is affected by oscillations with relative minimum corresponding to integer multiples of the perturbation period N b (accordingly with Figure 12).
Moreover, the phase lag ( Figure 5, bottom panel), calculated via a cross-phase analysis at the frequency of the perturbation f b = 23 mHz, close to zero in correspondence of the intervals 0 − N b , 2N b − 3N b , etc. The dashed line marks N w = 75. Conversely, in a separate analysis, we do not found any significant phase lag between B(t) and B EMD (t).

Energy conservation test
Basing on the previous results it is now clear that the MAVG affects the rotation matrix coefficients. In principle, the total energy should be conserved under rotation procedure. Numerically, it should correspond to a very low difference between power spectra of the original and that of the rotated signals. In order to evaluate the energy difference induced in the MFA aligned components by the MAVG procedure, we computed a Monte Carlo test comparing the power spectra of the time series in the original coordinate to that of the MFA reference frame. We assumed that the B x and B z components of the original signal are constant values (5 and 6 a.u. respectively) with added red noise, with AR(1)=.95 (autocorrelation coefficient at lag=1) and the standard deviation v red = 0.1 respectively, that are generated at each run, while the B y component is the same signal of Figure 4 except for the red noise generated at each run with properties just above described. In this test, the satellite position is fixed at the coordinates x = 1, y = 0 and z = 1 a.u. Figure 6 shows the absolute values of the averaged power difference |〈DP〉|, computed over 10 surrogates, as a function of N b and N w . The results show behaviors similar to |g| (see Figure 12) that affects the averaged perturbation 〈b d 〉, and hence |〈DP〉|. Indeed, |g| rapidly increases for N w < N b (see Appendix A for details), that corresponds to a violation of the inequality condition in the left side of relation (6).  Figure 3; (b) the difference DB 0 (t) between the original B 0 (t) and that estimated from the EMD (green) and MAVG (blue), where we used a window size of 75 points for the MAVG method; (c) the resulting average 〈DB 0 (t)〉 for 10,000 signals, differing only for a different noise; (d) detail of the phase shift between original perturbation and the averaged 10,000 signals [ b(t)+n(t)] estimated from different methods. In order to evaluate the effects of the noise, we computed a Monte Carlo test using 10,000 different noise keeping fixed the frequency f b = 23 mHz and the moving-average windows size N w = 75; the resulting spectra are averaged to obtain the total PSD ratio between rotated and original signals. The results are shown in Figure 7. EMD method gives excellent results, showing the same pattern of original signal through the frequencies. Conversely MAVG method presents several undesired issues: lowering energy at low frequencies, rising energy around the perturbation frequency, a pronounced peak in correspondence of the first harmonic of the perturbation and a regular oscillation at high frequencies. In particular, the last behavior depend on the window size of the moving average, with a characteristic frequency that increase with decreasing N w .

The MFA rotation matrix instability induced by MAVG procedure: a numerical example
We now discuss how the rotation matrix coefficients depend on the MAVG and EMD. We remark that, generally, these coefficients depend also on the mean field variation. In order to simplify our discussions, we assume a constant ambient B 0 = (0, B 0 , 0). In addition, we assume that the measuring point r = (x, 0, 0) is fixed. The chosen conditions will be clear in the next discussions. With these assumptions and in absence of any perturbing signal, each component of the rotation matrix is time invariant. If we assume b = (0, b 0 sin(~hi), 0) as the perturbation (t = hi) using Equation (14) the rotation matrix (13) is given by where We can see that, generally, G i switches from −1 to 1, so that the rotation matrix could be unstable inducing undesired features. To avoid this problem, the following condition must be satisfied: The resulting MFA components are We can see that, with our assumption, the MFA p q r q q p q q q q q q q q q q q q q q q q q q q q q q r q q $ t p q r q q p q q q q q q q q q q q q q q q q q q q q q r q q q q r p q r q q p q q q q q q q q q q q q q q q q q q r q q q q r t u v u u t u u u u u u u u u u u u u u u u u u u u u v u u u u v $ t u v u u t u u u u u u u u u u u u u u u u u u v u u u u v   As discussed in Section 2 (see also Appendix A for details), |g| decreases with N w , with local minima at N w = cN b (c = 1, 2,...,c max ). In the following test we used B 0 =0.3 nT, b 0 =1 nT, ~= 33 mHz, and h = 1 s, so that N b~3 0 and B 0 /b 0 = 0.3. The time series has a duration of 3600 s (i.e. i = 1, 2, …, 3600), and is analyzed using FFT algorithm. In the frequency domain we used the relationships (21) and (29) in order to select two values of N w equal to 11 and 35: those values correspond to |g| N w =11 = 0.78 > 0.3 (i.e. unstable condition) and |g| N w =35 = 0.03 < 0.3 (stable condition), as it can be seen in the top panel of Figure 8.
While the stable rotation matrix satisfies the energy conservation (i.e. G i = 1 = const, lower panel), the unstable condition does not preserve energy and G i is represented by an asymmetric square function and the negative values of B y (GSE) become positive in the resulting B µ , since G i switches from −1 to +1 with time.
Regarding the unstable condition, we show in Figure 9 the total power spectra of the original and rotated components, together with the G i spectra. We observe several features: (a) energy attenuation in the rotated components (at the signal frequency of 33 mHz and at lower frequencies); (b) the appearance of new peaks, in correspondence to multiples of the fundamental frequency (i.e. 66 mHz, 99 mHz,…). In a separate analysis (not shown here) we observed that the energy of the even harmonics decreases with decreasing ratio B 0 /b 0 , while the energy peaks related with the odd ones increase.
In particular, assuming that the mean field is zero (B 0 = 0), G i becomes a symmetric square wave function, and the fundamental frequency disappears, together with its even harmonics, while only the odd harmonics survive. Although this corresponds to an extreme case study, it is interesting to analyze because it allows us to highlight some of the issues covered here. Indeed, in this case the Fourier expansion of B µ = b 0 sin(~hi) is (see Appendix C) that clearly contains only the odd harmonics of the fundamental frequency of the signal.
Regarding the energy attenuation at lower frequencies, it can be tentatively explained by means of the Fourier expansion of B µ = G i [B 0 + b 0 sin(~hi)], obtained using the Fourier expansion of a generic square wave G i (see Appendices B and C) with amplitude A = 1 For a practical purpose, we truncate the Fourier expansion at order k = 0, obtaining c 0 B 0 + c 0 b 0 sin(~t) = c 0 B y . Since A = 1, for each coefficient it results |c k | 2 < 1, and hence, indicating with P 0 = P By (t) the power of the original signal, the resulting rotated component becomes P B µ = P c 0 B y = |c 0 | 2 P 0 ≤ P 0 . This implies that only for c 0 = 1 the power is equal (i.e. P B µ = P 0 ), and this is possible only if x + = 1 and x − = 0, and hence G(t) = 1 = constant (see Appendix B, Equations (33), (34), (36)), that corresponds to the validity of inequality condition (21). It is possible to study the Fourier series of B µ following Appendix B for any k-order, but it is out of the scope of this work.
The energy alteration, due to the MAVG application, can be regarded as a fundamental element of diagnostic in the rotation procedure, as we will discuss further in the next sections, by means of an experimental example.

A case study of the upstream waves event on May 6, 2005, observed by Cluster satellites
In order to compare EMD and MAVG effects on real data, we present a case study. Basing on the theoretical results showed in the previous section, we expect that major differences occur during time intervals characterized by large amplitude oscillations of the geomagnetic field (i.e. during disturbed geomagnetic conditions).
However, it is extremely difficult to find a magnetospheric event characterized by a large amplitude oscillation with respect to the mean magnetic field. Conversely, in the foreshock region, i.e. regions upstream the Earth's bow shock, the ULF upstream waves (20-100 mHz frequency range) amplitude is comparable with the interplanetary magnetic field strength. In that region, when the interplanetary magnetic field makes an angle with the bow shock normal direction i n,B < 45° (see Regi et al. [2014aRegi et al. [ , 2014b), the existing waves in the solar wind are amplified by means of the ion-cyclotron resonance mechanism.
We searched for a time interval characterized by quasi-monochromatic ULF fluctuations in the transversal components with respect to the ambient interplanetarymagnetic field, inside a database of previously examined events in Regi et al. [2014a]: in that work each event was classified by the signal-to-noise ratio, in each GSE component, using Cluster satellites fluxgate magnetometer data at 4 s sampling period [Balogh et al. 2001]. We selected the time interval 9-11 UT on May 6 (DoY 126) of 2005, characterized by a high signal-tonoise ratio and a low average ambient field in the same component. This condition should emphasize the effects induced by MAVG procedure, theoretically predicted and described above. Figure 10 shows the magnetic field components in the GSE reference frame (top panel), with oscillation amplitudes up to ~3 nT, while the average magnetic field was B 0 = (−4.8,−0.2, 2.1) nT. We computed the MFA components by means of both EMD and MAVG methods, also evaluating the absolute values of the difference between homologous components µ, z, o), using N w = 5 (unstable condition), as showed in the bottom panel of Figure 10. It can be seen that, during time intervals characterized by higher oscillations of B y (i.e. component with lower average 〈B 0,y 〉 = 0.21 nT), the difference becomes more evident. Similar, although less clear, results are obtained using greater N w value in the MAVG procedure. In Figure 11 we show a comparison between total power spectra (i.e. S GSE ( f ) = S x ( f ) + S y ( f ) + S z ( f )) of original magnetic field in the GSE coordinate system and that obtained with EMD and MAVG (i.e. It can be seen, for S GSE (marked in both panels of Figure 11 with black curves), a power peak at f uw~2 5.8 mHz (T b~3 8.7 s, corresponding approximately at N b = 9 samples). In this test, we used a number of weights N w = 5, 9, 19 for the MAVG procedure, while for the EMD we fixed a cut-off frequency of 1 mHz. Remarkable correspondence is found between S EMD (green line) and S GSE (black line). Conversely, S MAVG departs in the lower frequency range (i.e. f < f uw ) for any value of N w , while for higher frequency and for N w = N b , 2N b the correspondence is good. Remarkable difference is found around f~2 f uw , 3 f uw using N w = 5 (i.e. lower than N b ), accordingly to the theoretical discussions in previous section.
Further evidences of energy alteration as well as phase difference in the MAVG procedure are from multiple coherence analysis. Since no external noise is introduced in the rotation procedure, a multiple coherence c 2 close to 1 is theoretically expected. We computed c 2 between original GSE time series (3 inputs) and that in the MFA reference frame (3 outputs), obtained with both EMD and MAVG techniques (Figure 11, bottom panel). It can be seen that c 2 EMD is close to 1 at any frequency, while clear differences emerge in the c 2 MAVG at several frequency bands, further demonstrating that MAVG alters the power spectra of the rotated components with respect to that of the original ones.

Discussions and conclusions
In this work we studied the effect of the moving-average procedure applied to a discrete time series, comparing the results with those obtained using a new method based on the empirical mode decomposition. Theoretical investigations on a simple periodic and monochromatic perturbation, show that the MAVG procedure cannot completely remove the high-frequency perturbation from the mean field, affecting the rotated MFA components, although they are time-scaled separated.
The choice of the window size N w in the MAVG procedure strongly affects the estimation of the rotation matrix coefficients. We showed that the amplitude of the perturbation b(t), in the estimated B 0 (t), are close to but not equal to zero only if the window size N w corresponds to a integer multiple of the number of sample N b , that corresponds to the fundamental period of b(t) (i.e. N w = kN B , k = 1, 2, …). Moreover, a phase shift of 180° is introduced by the MAVG procedure when using a window size N w in the range between an odd and even integer multiple of the perturbation fundamental period N b (i.e. 2(m − 1) N b ≤ N w ≤ 2m N b , m = 1, 2, 3...). These results are also supported by means of energy conservation test, using a Monte Carlo simulation.
We also showed that theoretically, when the perturbation amplitude is large with respect to the mean magnetic field, several power peaks, at frequencies multiple of the fundamental one, appear. However it is extremely difficult to found such conditions in the magnetosphere, where the ambient field is higher than ULF waves amplitudes (see for example Sergeev et al. [2003]) and the stability condition is satisfied.
Finally, we studied the MAVG and EMD effects on real magnetic field data measured by Cluster satellite on May 6, 2005, in the upstream region. In this case, we found significant differences between the power spectra of the original and that of rotated components, using the MAVG. In other words, the MAVG operates as a low-pass filter with (generally) a non-zero phase lag. The choice of appropriate N w is crucial; it should not be too large, in order to be comparable with the ambient field fluctuation period, and not too small, in order to ensure that the power spectra of the resulting rotated MFA components are not modified, especially in correspondence to frequency ranges of our interest. Moreover, in some cases, the involved phenomena are non-linear and non-stationary, and hence the MAVG cannot be applied to identify the long-periodicities components of the discrete time series. A new technique based on the EMD is presented; this technique is applied to artificial (but almost realistic) signals, as well as to real data: The results show that it allows us to identify the long periodicities component B 0 without phase shift, and energy modifications.

MFA COORDINATES USING EMD
We further remark that, in theMAVG procedure, an ideal choice of N w equal to (or multiple of ) N b is not possible because: (a) the periodicity of the perturbation T b generally does not correspond to an integer number of samples, i.e. T b /h is not an integer number (h is the sampling period); (b) real data, in the frequency domain, generally may contain several spectral peaks, each characterized by different N b , and then N w cannot be uniquely determined in order to suppress, with the same efficiency, all the oscillations.
Conversely, the EMD procedure does not introduce any significant distortion in the MFA components, probably due to the adaptive nature of the EMD. However, it cannot resolve the most complicated cases when, for example, the noise and signals time-scale separation does not exist (see Wu and Huang [2004] for details). The complete codes can be requested via email at the following addresses: mauro.regi@aquila.infn.it; alfredo.delcorpo@aquila.infn.it.

A. Moving-average applied to a sinusoidal discrete time series
Let s i be a generic data series, the balanced moving average 〈s〉 at i th -central interval is defined as where N is the number of data on the left and on the right with respect to the central value, used to compute the average; hence the total number of data values used in the average is the odd integer N w = 2N + 1. Assume now a periodic and monochromatic form for the series s k = s 0 sin(~hk + a), where ~is the pulsation, a is the phase and h is the sampling period. Under these assumptions, Equation (24) at i th -time can be written as follows In order to compute the sum (25) we used the trigonometric superposition relationships (prosthaphaeresis) formula (26) obtaining the following terms: Therefore Equation (25) becomes: In this expression one term is independent and another (containing i-index) is dependent on time. More-over, assuming that the discrete time series is sampled at h resolution, we can define the integer number of samples N b that corresponds to the signal periodicity T s such that T s = hN b , allowing us to define the angular frequency as ~= 2r/T s = 2r/hN b . After the substitution of this expression in Equation (27), we obtain: where we have conveniently defined the time independent term Equation (28) has a time-dependent term and a time-independent dimensionless term g. The former is essentially the original time series (or perturbation) s i , and the latter is a superposition of cosine functions depending only on N w and N b . Figure 12 shows the absolute value of |g| (in a color scale) as a function of N w and N b , for a sampling period h = 1 s. |g| is always different from zero except for N w = kN b (k = 1,2,3,…,k max ), is indicated by white dashed lines. Here, k max is computed as [N s /N b ], where N s is the total number of samples, and [] is the integer part operator. In particular, |g| rapidly increases for N w < N b , exceeding 0.5; this condition corresponds to a violation of the inequality condition in the left side of Equation (6) (see also Section 5).

B. Fourier expansion of a generic square wave
Assuming an asymmetric square wave g(t) as showed in Figure 13 with a fundamental period T, resulting from a superposition of a positive pulse, with amplitude +A and duration x + T, and a negative pulse with amplitude −A and duration x − T (it can be noted that x + + x − = 1). Defining ~= 2r/T, the Fourier coefficients can be calculated as follows: where c + k are the coefficients for the positive pulse, calculated as while for the negative pulse, the coefficients are calculated as and k = 0, 1, 2, 3, … For Equation (31) we easy obtain (~= 2r/T) where we used the relation 2jsinx = e jx − e −jx . In order to calculate the contribution of the negative pulse, we use the properties of the Fourier transform of a time shift in the time domain F[x(t ± x)] = X(j~)e ±jk~d , where X(j~) = is the Fourier transform of the non shifted, negative pulse. In our case the time shift is d = (x + + x − )T/2 = T/2. Finally, for a time shifted negative pulse we obtain Using the equality e −jkr = (−1) k , and assuming A = 1, the resulting coefficients for g(t), where k = ±1, ±2, ±3,... are and hence the Fourier expansion of g(t) is From Equations (33)-(36) the Fourier expansion has both odd (k = 1,3,5,…) and even (k = 2,4,6,…) frequencies.

C. Fourier expansion of the product of the sine function with its associated square wave
Let assume a sinusoidal signal in discrete form s(t i ) = s 0 sin(~hi) (a = 0) at time t i = hi, we define the associated square wave q s (t i ) as follows where n = 1,2,3,... In order to compute the product s(t i )q s (t i ), we used the Fourier series expansions of the q s (t i ) odd function: The product ŝ(t) = s(t)q s (t) assumes the following expression Grouping the terms with the same argument in Equation (38)