Introduction

The onboard atomic frequency standard (AFS) is a fundamental element of a GNSS satellite, it plays an important role in the determination of the user position and characterizing the timing performance of the GNSS satellites, thus the satellite clock stability directly affects the range measurement and limits the navigation accuracy (Senior et al. 2008; Cernigliaro et al. 2013a, b; Hauschild et al. 2013). Thus, it is meaningful to establish the GNSS clock models with high precision for the GNSS positioning and timing capabilities.

In GNSS, each satellite of the constellation carries more than one atomic clock but one controls all the timing operations, such as signal generation and broadcasting, GNSS satellites employ different clock technologies: the GPS uses Cesium clocks and Rubidium Atomic clocks, GLONASS uses Cesium clocks, whereas Galileo and the BeiDou use RAFS and passive hydrogen masers (PHM) (Cernigliaro et al. 2013a, b; Hauschild et al. 2013). The current GPS constellations are running four families of clocks as shown in Table 1: the Block IIR and Block IIR-M satellites carry rubidium clocks, exclusively (Hauschild et al. 2013). The newest generation of Block IIF satellites uses improved Rubidium clocks, which has demonstrated superior stability in ground tests, and improved Cesium clocks (Vannicola et al. 2010; Weiss et al. 2010; Montenbruck et al. 2012). Finally, the next generation of Block IIIA satellites will be launched in near future (Marquis and Reigh 2015).

Table 1 Satellite information of the GPS system (as of July 2, 2016)
Full size table

The International Global Navigation Satellite Systems Service (IGS) operates a global network of Global Navigation Satellite System (GNSS) ground stations, data centers, and data analysis centers to provide data and derived data products to satisfy the objectives of a wide range of applications and experimentation, such as earth science research, multi-disciplinary positioning, navigation, and timing (PNT) applications; and education (IGS 2013). The IGS published the Broadcast, Ultra-Rapid, Rapid and Final products. Among them, the IGS Final products have the highest quality and internal consistency of all IGS products, which provide an opportunity for detecting the GPS clock characteristic with higher accuracy and finer temporal resolution (Senior et al. 2008).

Periodic fluctuations are apparently present in almost all GPS signal clocks (Gonzalez et al. 2010; Gonzalez Martinez 2013; Galleani and Tavella 2016). Some studies attribute these periodic fluctuations to the residuals of imperfect orbit estimation or the influence of periodic factors such as temperature and environment (Heo et al. 2010; Sesia et al. 2011; Zhou et al. 2011; Cernigliaro et al. 2013a, b). The periodic variations are very important because they can be used to improve the clock modeling, including interpolation of tabulated IGS products for higher-rate GPS positioning and predictions in real-time applications. This is especially true for high-accuracy uses but could also benefit the standard GPS operational products (Senior et al. 2008). Some tools are used to detect these periodic fluctuations, such as the dynamic allan variance, power spectrum and time–frequency spectrum (Galleani 2008; Sesia et al. 2011; Amin et al. 2017).

Conventionally, the periodic and time-varying periodic signals also called the harmonics and time-varying harmonics, are modeled by least-squares fitting harmonic terms with constant amplitudes and frequency. However, the time-varying harmonics have time-variable amplitudes and phases (Chen et al. 2013). Band-pass filtering, such as Fourier transform (Popinski and Kosek 1995), wavelet transform (Liu et al. 2007; Chao et al. 2014) and zero-phase digital filter (ZPDF) are often used to extract the periodic or time-varying periodic with time-variable amplitude and frequency (Höpfner 2004). But the time-varying periodic obtained by these methods will be corrupted easily by noise components, especially in the case when a large amounts of noise exist in the mixed signals because these methods choose a frequency band range in which may exist a lot of noises. In addition, other spatio-temporal analysis methods are also used to extracting the time-varying periodic signals, such as the singular spectrum analysis (SSA) (Chen et al. 2013), and empirical mode decomposition (EMD) (Smylie et al. 2015). However, when implementing the SSA, a key problem is choosing an appropriate lag-window size that determines the spectral resolution of the algorithm (Chen et al. 2013); EMD has the problem of mode mixing (Wu and Huang 2009). To overcome these problems, we present the inaction method (IM) and its principle. The inaction method, which is essentially line-pass filtering, can avoid polluting by noises to the maximum extent when extracting the time-varying periodic signals.

In GPS satellites, harmonics are a well-known feature since the early estimations of GPS clock (Swift and Hermann 1988; Senior et al. 2008; Gonzalez Martinez 2013). Four harmonic frequencies have been detected in all GPS clock types, n × (2.0029 ± 0.0005) cycles per day (24-h coordinated universal time or UTC), for n = 1, 2, 3, and 4 (Senior et al. 2008).

The GPS satellite clock data used in this study comes from the IGS Final 300 s satellite clocks (ftp://cddis.gsfc.nasa.gov/pub/gps/products/), active during 8 weeks from May 8 to July 2, 2016. The satellite information can be seen in Table 1 for this time range. PRN04 is not included in this table because it contains lots of missing data or abnormal values during this time range.

In the next section, we introduce the concept of the normal time–frequency transform and present the inaction principle which is the theoretical basis of the IM. Subsequent sections compare the extraction capacity of IM, SSA, and ZPDF from simulation time series data, and then we focus particular attention on the amplitude variations of 6- and 12-h terms extracted by the IM from the GPS satellite clock data over time. Finally, we compare the Hadamard deviation of the satellite clocks before and after exacting the 6- and 12-h terms.

Normal time–frequency transform

Time–frequency analysis can identify the instantaneous frequency, phase and amplitude of the undulation, which has been used in detecting periodic fluctuations. For unbiased measurement of the instantaneous frequency, phase and amplitude of a time series, the normal time–frequency transform (NTFT) has been proposed by Liu and Hsu (2009, 2011). Refer to these studies for detailed proof and properties of the NTFT explained below.

Definition

For time function \(f(t)\), where \(f{\text{:}}\;R \to C\), its time–frequency transform

$$\Psi f(\tau ,\varpi )=\int_{R} {f(t)\bar {\psi }} (t - \tau ,\varpi ){\text{d}}t, \tau , \varpi \in R$$
(1)

is called a NTFT, if the Fourier Transform of the transform kernel \(\psi (t,\varpi )\)

$$\hat {\psi }(\omega ,\varpi )=\int\limits_{R} {\psi (t,\varpi )\exp ( - i} \omega t)dt$$
(2)

satisfies

  1. (i)
    $$\hat {\psi }(\omega ,\varpi )=1,\;{\text{if }}\omega =\varpi$$
    (3)
  2. (ii)
    $$|\hat {\psi }(\omega ,\varpi )|<1,\;{\text{if }}\omega \ne \varpi$$
    (4)

where \(\tau\) is the time index and \(\varpi\) the frequency index, the line ‘–’ denotes the conjugate operator, C is the complex field and R the real field; and “||” is the modulus operator.

Satisfying (1), a typical NTFT kernel can be constructed as:

$$\psi (t,\varpi )=\left| {\mu (\varpi )} \right|w(\mu (\varpi )t)\exp (i\varpi t),w(t) \in \Omega (R), (\mu (\varpi ) \in R) \ne 0$$
(5)

where \(\Omega (R)\) is the function space of all the normal windows \(w(t)\). The rescaling index \(\mu (\varpi )\) can be almost anything other than zero, for instance, \(\varpi\), \({\varpi ^{2/3}}\), \({\varpi ^5}\) and so on. Particularly, if letting \(\mu (\varpi )=1\), Eq. (2) yields a phase-updated Gabor transform (GT) and letting \(\mu (\varpi )=\varpi\) yields a normal wavelet transform (NWT).

$$\Psi f(\tau ,\varpi )=\left| \varpi \right|\int_{{ - \infty }}^{{+\infty }} {f(t)w(\varpi (t - \tau ))\exp ( - i\varpi (t - \tau ))dt } \tau ,\varpi \in R$$
(6)

In (5) and (6), the local time function \(w(t) \in {L^1}(R)\) is called a normal window whose Fourier Transform satisfies

$$\left\{ {\begin{array}{*{20}{c}} {\hat {w}(\omega )=1,}&{if \omega =0} \\ {\left| {\hat {w}(\omega )} \right|<1, }&{ if \omega \ne 0 } \end{array}} \right.$$
(7)

The local time function \(w(t)\) often adopts a normal Gauss window

$$w(t)=\frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( { - \frac{{{t^2}}}{{2{\sigma ^2}}}} \right)$$
(8)

where \(\sigma>0\). The index \(\sigma\) is a vital parameter, which affects the time–frequency resolution of NTFT. We adopt the normal Gauss window and set \(\sigma =36\) when plotting the normal time–frequency spectrum of the GPS satellite clock data (Figs. 8, 9, 10, 11).

Inaction method

To obtain the harmonics and time-varying harmonics directly from the NTFT, Liu and Hsu (2009, 2011) and Liu et al. (2016) proposed the inaction method, which is essentially a forthright estimation of the harmonics and time-varying harmonics. Its theoretical derivation is as follows:

Inaction principle A time-harmonic \(h(t)\) can be written as

$$h(t)=A\exp (i\beta t)=\left| A \right|\exp (i(\beta t+\Delta ))$$
(9)

Applying an NTFT \(\Psi\) to \(h(t)\) yields

  1. (i)
    $$|\Psi h(\tau ,\varpi )|={\text{Maximum}}=|h(\tau )| \Leftrightarrow \varpi =\beta ,{\text{ }}\forall \tau \in R$$
    (10)
  2. (ii)
    $$\Psi h(\tau ,\beta )=h(\tau ),{\text{ }}\forall \tau \in R$$
    (11)

where A is the complex amplitude, the constant \(\beta \in R\) is the angular frequency, constant, \(\Delta \in [0,2\pi )\)is the initial phase, \(\beta t+\Delta\) is the immediate (or local) phase, “\(\forall\)” means “for any”, and “\(\Leftrightarrow\)” means “if and only if”.

Proof

The kernel of the NTFT satisfies

  1. 1.
    $$|\hat {\psi }(\omega ,\varpi )|={\text{Maximum}}=1 \Leftrightarrow \omega =\varpi$$
    (12)

Applying an NTFT to \(h(t)\), we have

$$\Psi h(\tau ,\varpi )=\int_{R} {A\exp (i\beta t)\bar {\psi }} (t - \tau ,\varpi ){\text{d}}t$$
(13)

If \(t - \tau ={t^\prime }\), \(t={t^\prime }+\tau\), then

$$\begin{gathered} \Psi h(\tau ,\varpi )=\int_{R} {A\exp (i\beta ({t^\prime }+\tau ))\bar {\psi }} ({t^\prime },\varpi ){\text{d}}{t^\prime } \hfill \\ =A\exp (i\beta \tau )\int_{R} {\exp (i\beta {t^\prime })\bar {\psi }} ({t^\prime },\varpi ){\text{d}}{t^\prime } \hfill \\ =h(\tau )\int_{R} {\bar {\psi }({t^\prime },\varpi )} \exp (i\beta {t^\prime }){\text{d}}{t^\prime } \hfill \\ =h(\tau )\bar {\hat {\psi }}(\beta ,\varpi ) \hfill \\ \end{gathered}$$
(14)

Thus

$$|\Psi h(\tau ,\varpi )|=|h(\tau )||\bar {\hat {\Psi }}(\beta ,\varpi )|=|h(\tau )|\hat {\Psi }(\beta ,\varpi )|$$
(15)

Considering (12) and (15), we have that, if and only if \(\varpi =\beta\), \(|\Psi h(\tau ,\varpi )|={\text{Maximum}}=\;|h(\tau )|,\;\forall \tau \in R\). Thus, the fact (10) is proven. Noting (3) and (14), we have that \(\Psi h(\tau ,\beta )=h(\tau ),{\text{ }}\forall \tau \in R\). Thus, the fact (11) is proven.

The result as expressed by (10) and (11) is called the inaction principle of the NTFT. It means the harmonic becomes available automatically from its NTFT without the inverse transform. The IM estimates the harmonics and time-varying harmonics in a time signal according to the inaction principle. Equations (10) and (11) mean that for a complex harmonic, its NTFT along the spectral ridge is itself. For a time-varying harmonic, its NTFT along the spectral ridge is an estimation or approximation of itself. As a matter of course, for real harmonic or a real time-varying harmonic, the twofold of the real part of its NTFT along the spectral ridge is itself or estimation of itself. On the time–frequency plane of the NTFT spectrum, the minimum zones are applied in the estimation of the harmonics and time-varying harmonics. Relating the band-pass filtering, the IM is essentially line-pass filtering. Thus, the IM can avoid being polluted by noises to the maximum extent and then have the ability to extract the time-varying harmonics with robustness; thus, IM is very suitable to extract the time-varying harmonics from the mixed signals, especially when it includes a lot of noise components. In this case, extracting time-varying harmonic terms using any type of band-pass filtering will be inevitably polluted by noises.

Simulation test

The methods used most often in extracting the harmonics and time-varying harmonics are SSA and ZPDF (Höpfner 2004). The details of SSA and ZPDF can be found in Chen et al. (2013) and Gustafsson (1996), respectively. To compare the extraction accuracy of these two filtering methods with IM, we simulate a set of time series which contain two time-varying harmonics and various noises. The sampling interval is 1/12,

$$y(t)=5{A_1}(t)\cos \left( {\frac{{2\pi t}}{{12}}+\frac{\pi }{3}} \right)+5{A_2}(t)\cos \left( {\frac{{2\pi t(1+0.00006(t - 600))}}{6}+\frac{\pi }{6}} \right)+{\text{noises}} 0 \leq t \leq 1200$$
(16)

where the first and second terms on the right-hand side of the equation are two original sub-signals, namely 12- and 6-h terms, respectively, which are shown in Fig. 1 (middle and bottom panels); the time-varying amplitudes of these two original sub-signals are \({A_1}(t)\) and \({A_2}(t)\) and can be seen in the top panel. It should be noted that the 6-h term has time-varying frequencies. The third term contains various noises, such as White phase modulation (WPM) noise, Flicker phase modulation (FPM) noise, White frequency modulation (WFM) noise, Flicker frequency modulation (FFM) noise, Random Walk frequency modulation (RWFM) noise, Flicker Walk frequency modulation (FWFM) noise, Random Run frequency modulation (RRFM) noise, and mixture noise, which are shown in Fig. 2, respectively.

Fig. 1
figure 1

Simulated 12- and 6-h terms and their amplitude

Full size image
Fig. 2
figure 2

Simulated noises

Full size image

To compare the robustness of these three methods regarding extracting the harmonics and time-varying harmonics, we simulated the synthesized signals with different signal-to-noise ratio (SNR). As shown in Fig. 3, the synthesized signals contain 12- and 6-h terms, and various types of noises.

Fig. 3
figure 3

Synthesized signals containing the same 12- and 6-h terms and various noises, respectively. The component in the bracket is SNR of the synthesized signals with a different value

Full size image

As shown in Fig. 3a, b, c, d, h, the simulated 12- and 6-h terms are almost covered by the noise. The extracted terms by these three methods can be seen in Figs. 4, 5, 6 and 7. Rows from top to bottom correspond to different types of noise, and the method adopted is labeled on each figure. In these figures, the blue and black lines indicate the extracted terms and the error between the true harmonics and the extracted terms, respectively.

Fig. 4
figure 4

Extracted 12-h terms (blue) and errors (red) in the case of the signals polluted by different types of noises. Rows from top to bottom correspond to the White PM noise, the Flicker PM noise, the White FM noise and the Flicker FM noise

Full size image
Fig. 5
figure 5

Extracted 12-h terms (blue) and errors (red) in the case of the signals polluted by different types of noises. Rows from top to bottom correspond to the Random Walk FM noise, the Flicker Walk FM noise, the Random FM noise and the mixed noise

Full size image
Fig. 6
figure 6

Extracted 6-h terms (blue) and errors (red) in the case of the signals polluted by different types of noises. Rows from top to bottom correspond to the White PM noise, the Flicker PM noise, the White FM noise and the Flicker FM noise

Full size image
Fig. 7
figure 7

Extracted 6-h terms (blue) and errors (red) in the case of the signals polluted by different types of noises. Rows from top to bottom correspond to Random Walk FM noise, Flicker Walk FM noise, Random FM noise, and the mixed noise

Full size image

As shown in Fig. 4a–i, the graphics of 12-h terms extracted by these three methods all seem to be similar when processing the synthesized signals, which includes White PM, Flicker PM and White FM, respectively. In fact, in this case, the extraction precision of IM is slightly better than SSA and ZPDF. The results of statistical error and correlation coefficient can be seen in Tables 2 and 3, respectively. As shown in Fig. 4j–l, the graphics of the 12-h terms extracted by SSA and ZPDF are apparently distorted when processing the Flicker FM noise. In this case, we see that the 12-h term extracted by IM is less affected by the Flicker PM noise compared to the methods SSA and ZPDF, which can be proven to be the results of statistical error (Table 2) and correlation coefficient (Table 3).

Table 2 Correlation coefficient between the estimation/extraction or the simulation harmonics
Full size table
Table 3 Relative errors of the estimation or extraction using the three methods
Full size table

As shown in Fig. 5a, j, the burr appeared in the extracted 12-h terms by SSA, which includes the Random Walk FM and the mixed noise, respectively. The 12-h terms extracted by ZPDF are apparently irregular as shown in Fig. 5b, k; however, in this case, the 12-h term extracted by IM is the smoothest among these three methods. In this case, it indicates that the 12-h harmonic is more accurately recovered by the IM than both the SSA and ZPDF when the signal includes Random Walk FM and mixed noise (as shown in Tables 2, 3).

As shown in the first row of Fig. 6, the errors of the extracted 6-h terms by ZPDF and IM are significantly less than by SSA when processing the synthesized signals, which include the White PM noise. The results are shown in Tables 2 and 3. The 6-h terms extracted by SSA and ZPDF are slightly distorted as shown in Fig. 6g, h, j, k. In the case, the 6-h term extracted by IM is the smoothest among these three methods. It indicates that the 6-h term extracted by IM is less affected by the White FM noise and the Flicker FM noise than those by SSA and ZPDF.

As shown in Fig. 7, top and bottom rows, the 6-h terms extracted by SSA and ZPDF are slightly distorted. In the case, the 6-h term extracted by IM is smoothest among these three methods. It indicates that the 6-h term extracted by IM is less affected by the Random Walk FM noise and the mixed noise than the terms extracted by SSA and ZPDF.

For evaluating the extraction accuracy of these three methods, the relative error was introduced as follow:

$${\text{Relative error}}=\frac{{{\text{RMSE}}}}{{{\text{signal's}}\;{\text{RMS}}}}=\frac{{\sqrt {\sum\limits_{{i=1}}^{N} {{{({x_i} - {{\hat {x}}_i})}^2}/N} } }}{{\sqrt {\sum\limits_{{i=1}}^{N} {x_{i}^{2}/N} } }},$$
(17)

where \({x_i}\) and \({\hat {x}_i}\) indicate the ith simulation value and ith estimation value, respectively. \(N\) is the length of the simulation data.

The correlation coefficient and relative errors can be seen in Tables 2 and 3, respectively. Because of the edge effect present in all of these three methods, some values at both ends of extraction time series are removed when evaluating the extraction accuracy.

As shown in Table 2, the correlation coefficient between the 12-h terms extracted by IM and the simulation harmonics are all greater than 0.99. As shown in Tables 2 and 3, the correlation coefficient between the extracted terms by IM and the simulation harmonics, and the accuracy of extraction by IM are all the best, except in the case of the Random Run FM only. This simulation demonstrates that the IM can robustly extract time-varying harmonics from a noisy time signal, even if the SNR is −21.

Time–frequency analysis of the GPS satellite clocks

To study the time–frequency characteristic of the GPS satellite clocks, the IGS Final 300 s satellite clocks as shown in Table 1 are adopted for the time May 8–July 2, 2016. The GPS satellite clocks with significant periodicity are chosen to analysis their time–frequency characteristic. In addition, the two GPS cesium clock, prn08 and 24, are also chosen for comparison. To reduce the effect of the trend term on the time–frequency spectrum as far as possible, quartic polynomials are used for removing the trend term of the GPS satellite clocks. According to the normal wavelet transform (8), the NTFT spectrum of the GPS satellite clocks can be obtained.

As can be seen in Figs. 8, 9 and 10, significant 6- or 12-h periodic signals exist in GPS rubidium clocks; most of the rubidium atomic clock contains both of the two periodic signals, especially, the block IIF rubidium atomic clock. Even some satellite clocks contain 3- and 4-h periodic signals, such as prn06 and prn21. As can be seen in Fig. 10, the GPS cesium clocks, prn08 and 24, do not contain these two periodic signals. It is possible that the short term of the GPS cesium clocks is worse than that of the rubidium clock (Cernigliaro et al. 2013a, b); in that case the periodic signal may be drowned by the noise of the GPS cesium clocks. These conclusions from above agree with Senior (2008).

Fig. 8
figure 8

NTFT spectrum of GPS satellite clocks (prn01, 03, 05, 06, 09, 12, 13 and 15)

Full size image
Fig. 9
figure 9

NTFT spectrum of GPS satellite clocks (prn16, 17, 18, 21, 22, 25, 26 and 27)

Full size image
Fig. 10
figure 10

NTFT spectrum of GPS satellite clocks (prn29, 30, 31, 32, 08, and 24)

Full size image

Amplitude variations of the harmonic

As can been seen in Figs. 8, 9 and 10, there are significant 6- and 12-h periodic terms are both of these satellite clocks, such as prn01, 03, 05, 06, 09, 12, 13, 15, 16, 17, 18, 21, 22, 25, 26, 27, 29, 30, 31 and 32. Then, according to the inaction method, equations 7 and 8, the 6- and 12-h periodic signals of these satellite clocks can be extracted from their NTFT spectrum, respectively.

As can be seen in Fig. 11, the amplitudes of the 6-h periodic terms of these satellite clocks show different variations. For example, within this time span, prn12, 13, 16 have a gradually decreased fluctuation, whereas prn21, 25, 26, 31 have a gradually increased fluctuation. In addition, the amplitudes of the 6-h terms of prn16, 17, 29 and 31 are relatively large, being greater than 0.2 ns.

Fig. 11
figure 11

6-h periodic terms extracted by IM from GPS satellite clock bias

Full size image

As can be seen in Fig. 12, the amplitudes of the 12-h periodic terms of these satellite clocks show significant variations, and amplitudes of the 12-h terms of prn05, 12, 13, 16, 17, 18, 27, 29, 30, and 31 are relatively large, being greater than 0.2 ns.

Fig. 12
figure 12

12-h periodic terms extracted by IM from GPS satellite clock bias

Full size image

As shown in Figs. 11 and 12, the amplitudes of the 12-h terms are generally greater than the 6-h terms. Figures 11 and 12 both show that their amplitudes of the 6- and 12-h terms vary over time, which should be taken into full consideration when modeling the GPS satellite clocks. In addition, we note that the amplitudes at the beginning and end of the 6- and 12-h terms may be less accurate because of an edge effect.

Comparison of the Hadamard deviation

To verify whether the 6- and 12-h terms are extracted completely and exactly from GPS satellite clock using the IM, the Hadamard deviation of the satellite clocks before and after extracting these periodic terms are compared. More information on the Hadamard deviation can be found in Hutsell (1996) and references therein.

As shown in Fig. 13, after removing the 6- and 12-h periodic terms, the oscillation phenomenon in Hadamard deviation virtually disappears, which demonstrate that these periodic signals have been extracted effectively. In addition, with respect to the red line, the blue line indicates that the Hadamard deviation of the satellite clocks within the averaging interval 300–3000s is the response of satellite clock noise, for averaging intervals beyond 3000s, it is the response of periodic signals existing in satellite clocks. The red lines demonstrate that the Hadamard deviation of the remaining signal after removing the periodic terms from satellite clocks is the real response of the clock noise within the averaging interval 90000s. Therefore, to evaluate the clock noise characteristic within the long-term averaging interval, it is possible to extract the periodic terms using the IM.

Fig. 13
figure 13

Comparison of the Hadamard deviation. AI is an abbreviation of averaging interval. The blue and red lines indicate the Hadamard deviation of GPS satellite clocks and the Hadamard deviation of the residual after removing the 6- and 12-h periodic terms, respectively

Full size image

Conclusions

We demonstrated the performance and robustness of the IM, which essentially is line-pass filtering, on extracting the harmonics and time-varying harmonics from a noisy signal. The 12- and 6-h periodic terms in most GPS rubidium clocks have been detected in NTFT spectrums, which can provide an unbiased measurement of the instantaneous frequency, phase and amplitude of the periodic signals. Through analyzing the 6- and 12-h periodic terms, we find that their amplitudes show significant variations between May 8 and July 2, 2016, which are very important to model satellite clocks. After removing 6- and 12-h periodic terms, the oscillation phenomenon in Hadamard deviation virtually disappears, this also shows that these terms have been extracted effectively by IM. Therefore, when evaluating the response of noise within long-term averaging interval using Hadamard deviation, one can detect the periodic signals by NTFT and remove these by the IM.