Power Spectral Density Background Estimate and Signal Detection via the Multitaper Method

We present a new spectral analysis method for the identification of periodic signals in geophysical time series. We evaluate the power spectral density with the adaptive multitaper method, a nonparametric spectral analysis technique suitable for time series characterized by colored power spectral density. Our method provides a maximum likelihood estimation of the power spectral density background according to four different models. It includes the option for the models to be fitted on four smoothed versions of the power spectral density when there is a need to reduce the influence of power enhancements due to periodic signals. We use a statistical criterion to select the best background representation among the different smoothing + model pairs. Then, we define the confidence thresholds to identify the power spectral density enhancements related to the occurrence of periodic fluctuations (γ test). We combine the results with those obtained with the multitaper harmonic F test, an additional complex‐valued regression analysis from which it is possible to estimate the amplitude and phase of the signals. We demonstrate the algorithm on Monte Carlo simulations of synthetic time series and a case study of magnetospheric field fluctuations directly driven by periodic density structures in the solar wind. The method is robust and flexible. Our procedure is freely available as a stand‐alone IDL code at https://zenodo.org/record/3703168. The modular structure of our methodology allows the introduction of new smoothing methods and models to cover additional types of time series. The flexibility and extensibility of the technique makes it broadly suitable to any discipline.

(1) where f j = j/(NΔt) are the Fourier frequencies defined over the frequency interval [−f Ny , f Ny ], limited by the Nyquist frequency f Ny = 1/(2Δt), with the frequency resolution determined by the Rayleigh frequency f Ray = 1/(NΔt). The periodogram is defined as the product of the time series sampling rate over the number of points and discrete Fourier transform square modulus, that is S (p) (f j ) = (Δt/N)|X j | 2 . For real-valued processes, the spectral density function is two-sided, i.e., symmetric about the zero frequency so that S(−f ) = S( f ). In this case, we can define the one-sided spectral density function, hereafter referred to as power spectral density (PSD), with S( f ) doubled for 0 < f < f Ny and set to zero for f < 0. As a consequence, the PSD is defined on N f Fourier frequencies f j with j = 0, 1, …, (N f − 1). Note that N f = N/2 + 1 for even N and N f = (N + 1)/2 for odd N. The major issues of this estimator are well known (Percival & Walden, 1993): (i) the leakage of power into adjacent bins, due to the finite frequency resolution, (ii) a bias in the estimate not known a priori, depending on the time series itself, and (iii) the associated variance, that is equal to the estimate itself. These effects can be reduced by tapering the time series with appropriate weights w n , satisfying   2 1 n n w , and/or by averaging the PSD over adjacent frequency bins (Percival & Walden, 1993). Another procedure consists of averaging the PSD estimated on different weighted subintervals (possibly overlapped) of the original time series (Welch, 1967); since the intervals are shorter, the frequency resolution is reduced. Additionally, many parametric spectral analysis procedures exist, including minimum prediction error (Samson, 1983), maximum entropy (Vellante & Villante, 1984), and CARMA (Kelly et al., 2014) method. Among the nonparametric methods, the singular-spectrum analysis (Ghil, 1997) and the adaptive multitaper method (Thomson, 1982) have been extensively used for the identification of periodic signals in time series. The singular-spectrum analysis is able to reconstruct the original data in terms of oscillatory components, based on a data-adaptive basis-set, obtained with the eigen-decomposition of the lagged covariance matrix on M lagged copies of a time series. This method is particularly useful for the analysis of nonlinear systems, owing to the absence of assumptions on the basis-set. On the other hand, it is difficult to recover the frequency of a reconstructed oscillation as the singular-spectrum analysis searches for frequency bands containing a relevant amount of the time series variance, rather than discrete PSD enhancements.
Purely periodic or quasiperiodic signals appear in the PSD as enhancements relative to the continuous PSD background whose properties depends on the physical system. While harmonic analysis to identify the occurrence of periodic variations compared to a flat PSD (i.e., white noise) are well established (Percival & Walden, 1993), there is no standard techniques to assess the significance of a periodicity against a colored noise, such as the red PSD typically found in astrophysical and geophysical time series. The identification of the continuous part of the PSD constitutes a great challenge since sharp peaks can be created by completely different processes, such as random/stochastic processes with signals or deterministic chaotic systems (Kantz & Schreiber, 2003). Vaughan (2010) addressed this issue in an analysis of the occurrence of quasiperiodic oscillations in X-ray observations of Seyfert galaxies (Vaughan, 2005). They introduced a significance test for periodicity assuming red noise PSD with an approximately power law or bending power law spectral background shape. Using the statistical properties of the periodogram, Vaughan (2010) applied a Bayesian approach to estimate the posterior distribution of the PSD model parameters. After selecting the best representation of the PSD background via the sum of the squared standard errors and the likelihood ratio test (Vaughan, 2010;Vaughan et al., 2011), periodic signals appear as periodogram outliers.
More recently, Inglis et al. (2015) adapted the Vaughan (2010) method to the identification of quasiperiodic pulsations typically observed during the impulsive phase of solar and stellar flares over a wide range of wavelengths. From radio waves and microwaves to hard X-rays and gamma-rays (Nakariakov & Melnikov, 2009), the characteristic timescales of these fluctuations range from 1 s up to several minutes. The Automated Flare Inference of Oscillations (AFINO; Inglis et al., 2015Inglis et al., , 2016 technique probes the PSD of the time series for a single power law-plus-constant model, a broken-power law model, and power lawplus-constant combined with a Gaussian component in log-frequency space, representing the excess power due to the occurrence of a periodic oscillation. The most appropriate background model is selected via the Bayesian information criterion (Burnham & Anderson, 2004) and a modified χ 2 statistic for exponentially distributed data (Nita et al., 2014). The AFINO technique has been applied also to magnetometer data from the Magnetospheric Multiscale mission to study the role of the ULF waves in the dynamics of the inner magnetosphere and outer radiation belt (Murphy et al., 2018). This technique has been proven to be effective in the identification of strong PSD enhancements, but it is limited to the selection of a single wave mode (Murphy et al., 2020). Mann and Lees (1996) proposed a procedure for distinguishing between PSD background and peaks, based on the spectral and harmonic analyses of Thomson (1982). Briefly, the PSD background is estimated fitting a lag-one autoregressive model to the median smoothed PSD of the time series. Then, a periodicity is identified at locations where the PSD enhancements are concurrent with enhanced harmonic F test values, both above a defined confidence threshold (Thomson, 1982). This method has been applied to many studies of remote and in situ observations of the solar wind and the magnetosphere. A similar approach was developed by Di Matteo and Villante (2017) who combined the identification of narrow peaks in the PSD, estimated with the Welch method, with the harmonic F test of the MTM method.
Here, we combine and improve some of these approaches. Following a brief description of the spectral and harmonic analysis of the multitaper method, we discuss the extension of the maximum likelihood approach, developed for the periodogram, to the multitaper estimates of the PSD. We introduce various combinations of PSD models and smoothing approaches. We discuss robust statistical criteria to determine the best representation of the PSD background. Finally, we describe different options for the identification of periodic fluctuations. The method is validated with Monte Carlo simulations and demonstrates its application with real observations.

The Multitaper Method
Given a time series of length N with sampling time Δt, the multitaper method (MTM) uses a set of K orthogonal tapers to obtain K independent estimates of the PSD. The tapers result from the Fourier transform of the eigenfunctions of the Dirichlet kernel, namely the Slepian functions (Slepian, 1978). These functions minimize the spectral leakage outside a frequency band 2W/Δt with 0 < 2W < 1. Ordering the Slepian sequences with the corresponding eigenvalues in decreasing order, the first K ≤ 2NW − 1 eigensequences have eigenvalues close to 1 (Slepian, 1978) and provide, in the case of a white noise process, unbiased and uncorrelated estimates of the spectral density function at the Fourier frequencies f j (Thomson, 1982), ( ) in which the weights d k (f j ) are derived from: where σ 2 is the variance of the time series. The weights are obtained at the Fourier frequencies f j by recursively substituting S( f ) with the spectral density function estimated by (2). In particular, starting from the average of the spectral estimates ( ) ( ) mt k j S f calculated using the first two ordered Slepian sequences, we obtain a set of weights from (3), that, when substituted into (2) gives a new estimate of the spectral density function to be used for the evaluation of d k (f j ). As with the periodogram, the PSD estimator is the one-sided spectral density function. The main advantage of this procedure is the attenuation of the average broadband bias, i.e., the amount of power leakage outside a frequency band of 2W/Δt (Percival & Walden, 1993;Thomson, 1982).
A powerful tool that we use in conjunction with the MTM spectral analysis is the harmonic F test. The assumption is that a time series can be expressed as a superposition of sinusoidal components and a DI MATTEO ET AL.

10.1029/2020JA028748
3 of 23 background process with a continuous PSD (Ghil et al., 2002;Thomson, 1982). The MTM yields a complex-valued regression model (Di Matteo & Villante, 2017;Thomson, 1982) from which it is possible to estimate amplitude and phase of the sinusoidal components. The null hypothesis, that an estimated amplitude is zero, is tested with the harmonic F test according to a Fisher distribution, which provides the confidence interval of the least squares fit. If the initial assumption is not valid and the PSD background is not locally white, false positives can be identified. Protassov et al. (2002) cautioned that the F test deviates from the nominal Fisher distribution when the null value of the tested parameters is on the boundary of the possible parameter value, as for the MTM harmonic F test. As a consequence, this test is likely to identify false positives, especially at low confidence levels. Therefore, this method should never be used alone. Mann and Lees (1996) considered only values of the harmonic F test above a defined confidence level that were concurrent with PSD enhancements with respect to a PSD background. This combined test is more robust than either test alone. In the following section, we extend this methodology through additional smoothing approaches and background models, the latter fitted via an appropriate maximum likelihood approach.

Maximum Likelihood and Confidence Bounds
While fitting a model to an estimated PSD, we have to consider the probability density function of these estimates since they are not Gaussian distributed. The periodogram estimates follow an exponential distribution, that is is the expectation value at the Fourier frequencies f j ≠0, f Ny (Anderson et al., 1990;Bevington & Robinson, 2003;Vaughan, 2005). The adaptive MTM estimates instead follow a gamma distribution (Thomson & Haley, 2014), such that S (amt) ( f j ) ∼ Gamma(α j , B j /α j ), where α j is related at each Fourier frequency to the number of degrees of freedom, ν j , defined as (Percival & Walden, 1993): where d k (f j ) are the final weights obtained from equation (3).
We extend the approach already adopted for periodograms (Vaughan, 2005(Vaughan, , 2010Vaughan et al., 2011) to MTM PDSs. Given a time series of length N, the joint probability density function, which characterizes the distribution of the PSD estimates at f j ≠0, f Ny , is L = ∏ j p(S j ). When used as a function of the model parameters, this corresponds to the likelihood function that can be more easily managed considering the log-likelihood, namely M = −2 ln L. Table 1 summarizes the types of random variables, the probability density functions, and the log-likelihoods for the periodogram and MTM estimates. Note that the two approaches match each other for α j = 1, corresponding to one direct PSD estimate among the ones obtained from the different tapered data instances.
Once the PSD background has been estimated, we define confidence thresholds in order to identify statistically significant PSD enhancements. In previous work, the ratio between the estimated PSD and the DI MATTEO ET AL.
10.1029/2020JA028748 4 of 23 If we consider the ensemble of γ j as possible representations of a single random variable γ, the corresponding probability distribution function is: where p(α) is the probability distribution function of half the number of degrees of freedom that we estimate via a simple histogram of the α j values over the range [0, K] with a fixed step of Δα = 0.2. The use of more sophisticated methods for the estimation of p(α), like the nearest neighbor or the kernel methods (Silverman, 1986), determine differences lower than the 1.0% on the final confidence level. To define a confidence threshold z, we need the cumulative distribution function. Considering that 0 < α < K by definition and that z > 0, since the PSD is always positive, we obtain: Introducing the normalized lower incomplete gamma function: the cumulative distribution function for the random variable γ is: At a given confidence level ϵ, a threshold z ϵ can be evaluated by searching for the zero of the function g(z) = C K (z) −ϵ.

Practical Procedure
Our procedure is freely available as a stand-alone IDL code at https://zenodo.org/record/3703168 (Di Matteo et al., 2020). We assume a time series x n regularly sampled with no data gaps. By default, we subtract the average value < x n >. Note that data trends, due to long term variations on the same timescale of the length of the interval, might affect the results. In this case, the user should consider prewhitening of the time series if necessary. In the following sections, we carefully describe our new procedure for the characterization of the PSD background and the identification of signals.

Smoothing
Accurate estimation of the PSD background can be strongly influenced by embedded signals that create large local enhancements in the PSD (signal-to-noise ratio of several units or more). The major consequence of this energy excess is to increase the estimated background level, possibly along the entire frequency range, leading to selection of PSD peaks at lower confidence levels. The smoothing of the PSD is a way to reduce this effect (Percival & Walden, 1993). Here, we describe the four different approaches we offer as options in our algorithm. The italicized abbreviation used below to refer to each approach corresponds to the keyword for calling that version of smoothing in our code.
Our first smoothing approach is the running median (med; Mann & Lees, 1996) over frequency intervals of 2w + 1 points.
Near the edges of the frequency interval the window is truncated to fewer points. The number of points, determined by w, are evaluated from a percentage value p of the available frequency interval. For example, given the complete interval [0, f Ny ] and the percentage value p (such that 0 < p < 1), the width of the smoothing window is 2w + 1 ≈ (pf Ny )/f Ray . Note that the running median strongly distorts portions of the PSD that exhibit steep variations.
Our second approach is based on a running median on windows with uniform width with respect to the central frequencies in the logarithmic frequency space (mlog; Stella et al., 1997), namely: For geophysical signals, which are typically red noise spectra, the critical range is at low frequencies (Mann & Lees, 1996). This approach includes only a few points at low frequencies enabling the recovery of the steep PSD at low frequencies. However, at high frequencies, where a large portion of the frequency range is included, the smoothed PSD tends to flatten.
The third smoothing approach associates the running average of the logarithmic PSD over 2w + 1 data points to the geometric mean of the corresponding frequencies (bin; Papadakis & Lawrence, 1993), namely: with k = j − w, …, j + w. At the edges of the frequency interval, we neglect intervals of length less than 2w + 1, so that j = w, …, N f − w − 1. Papadakis and Lawrence (1993) showed that this is an unbiased estimator of the true PSD at the set of frequencies f bin,j in the case of a power law. They also note that the bias is small as long as the logarithm of the PSD varies smoothly with the logarithm of frequency. Note that the bin smoothed PSD can be significantly distorted if the raw PSD exhibit strong local spikes.
For our fourth smoothing approach we apply a Butterworth low pass filter to the PSD as if it were a time series (but). The Butterworth gain function is given by: where f′ are the "frequencies," f′ c is the cutoff frequency, and Ω is the order of the filter. Typically, the filtered series exhibits problems at the boundaries of the interval. To overcome this issue, we first extend the data by introducing a mirrored replica of itself at both ends of the PSD. Then, we apply the zero-phase forward and reverse Butterworth filter providing no phase distortion. Finally, the central part of the inverse Fourier transform provides the smoothed PSD. The percentage of smoothing p regulates the value of the cutoff frequency f ′ c = pf ′ Ny while the order is set to Ω = 8. The choice of the filter order is arbitrary, but it provides reasonable results in various synthetic data representations (white and colored noise). For PSD with steep variations, this procedure shows limitations similar to the med approach. Note that the but smoothed PSD can be affected by strong local spikes in the raw PSD if they occur.
The parameter p of the smoothing procedure must be chosen carefully. The width of the window must be greater than the width expected for the PSD enhancements, but not too large to distort the PSD. For the MTM, the width of the peaks in the PSD is typically greater than 2W/Δt. This set the minimum size of p, such that p > 2W/Δtf Ny = 4W. To avoid strong distortions of colored PSD, we assume an upper limit of p = 0.5. To have an estimate of the optimum window in different scenarios, we can use the information on the probability density function of the PSD. Stella et al. (1997) showed that a Kolmogorov-Smirnov (KS) test (Press et al., 2007), can be applied to the ratio between the periodogram and its smoothing. In a similar way, for the MTM, we can apply the same concept to γ, as defined in equation (5). The data points can be converted to an unbiased estimator of the cumulative distribution function C γ (z) with z > 0 providing the fraction of data points less than a certain value z. The theoretical cumulative distribution function for the ratio γ is C K (z) as defined in equation (9). The KS test probes the similarity between these two cumulative distribution functions evaluating their maximum distance: The optimal percentage of smoothing is the one in which C γ (z) minimizes the D KS value. First, we probe all the p values between 4W and 0.5, then we select the p corresponding to the minimum D KS value among all the local minima. In the supporting information, we provide the distribution of the optimal p values obtained from a Monte Carlo simulation of synthetic time series. This procedure is robust for peaks with a signal-to-noise ratio on the order of unity, even in the case of multiple PSD peaks. When stronger signals occur, the smoothed PSD results are distorted, especially for the bin and but approaches. In this case, additional steps are required that will be discussed in section 4.4.

Background Models
Once the adaptive MTM PSD, hereafter referred to as the raw PSD, and the different smoothed versions med, mlog, bin, and but have been evaluated, we can test the PSD background against simple parametric models representative of a wide range of geophysical systems. The best parameters are determined by minimizing the log-likelihood as outlined in section 3.
A common representation of colored PSD is the power law (PL) model: with constant factor c and spectral index β. For β = 0, (15) reduces to a simple white noise process, where the power is evenly distributed among the Fourier frequencies f j . In this case, we can analytically determine the maximum likelihood, such that the PSD background is the weighted average of the adaptive MTM PSD estimates at the Fourier frequencies f j with weights equal to half of the corresponding number of degrees of freedom, namely For the power law model we use a numerical procedure to minimize the log-likelihood. We start from a rough estimate of the spectral index, as the slope of the logarithmic PSD, and the corresponding analytical solution for the constant factor, namely: where the indices j′ and j″ refer to the lower and upper limit of the frequency range of interest. To find the solution, the minimization procedure needs the definition of the parameter space boundaries. For the PL model, we need only the lower and upper limit of β. The default interval in our code is 0 < β < 10, but it can be modified by the user.
When considering discrete finite red noise time series, the simplest statistical process one can assume is the lag-one autoregressive process (AR(1)) represented by x n = ρx n−1 + w n . The present value of a time series x n depends on the past values x n−1 by the degree of serial correlation (the lag-one autocorrelation coefficient 0 ≤ ρ < 1) together with some random effect w n (white process with variance σ 2 ). It is representative of many geophysical systems (Mann & Lees, 1996). The autocorrelation of a AR(1) process decays exponentially with a characteristic time determined by τ = −Δt/log(ρ); therefore, on time scales longer than τ, it behaves as a white noise process. The corresponding PSD is given by (Mann & Lees, 1996;Vaughan et al., 2011): Note that for ρ = 0, (17) reduces to a white process. For the numerical minimization procedure, we define the starting values for ρ, using the Yule-Walker equation, and for c, using its analytical solution for the log-likelihood minimization, namely: In our code, the default interval for the lag-one autocorrelation coefficient is 0 < ρ < 1, but it can be modified by the user.
A more flexible approach is the adoption of analytical functions able to reproduce the general behavior of geophysical PSDs, even though they are not related to a particular stochastic process. An example is the bending power law (BPL; McHardy et al., 2004;Vaughan et al., 2011) defined as: There are four parameters: the constant factor c, the spectral indices β and γ dominating, respectively, the frequency intervals below and above the frequency break f b at which the model bends. This model is particularly helpful when analyzing time series of turbulent systems that exhibit different spectral indices at frequencies below and above a frequency break, corresponding to different regimes of the energy cascade. As in the previous models, we provide a starting value for the model parameters. We initialize the estimate with the frequency break at the center of the interval in analysis, the spectral indices as the slopes of the logarithmic PSD at frequencies below and above the frequency break, and the constant factor from its analytical solution for the log-likelihood minimization, namely: where the indices j′ and j″ refer to the lower and upper limit of the frequency range of interest. When β > γ > 0, the spectral indices β and γ dominate in the opposite frequency interval, that is above and below the frequency break, respectively. We can recover our original definition, with the parameter transformation of In our code, the default parameter space intervals for the BPL model are − 5 < β < 10, 0 < γ < 15, and 0 < f b < f Ny , but they can be modified by the user.

Best PSD Background Choosing Criteria
The combination of the possible smoothing and models creates an array of PSD background estimates that in some cases are very similar. Here, using the stochastic properties of the adaptive MTM PSD estimates, we introduce three tests providing objective criteria to identify the best representation of the PSD background. In the following, B j indicates a possible PSD background and S j the raw unsmoothed PSD.
Based on the likelihood and the number of free parameters, N θ , of each model, a useful method of comparison is the Akaike Information Criterion (AIC; Akaike, 1973): It corresponds to the sum of the log-likelihood with a penalty value for including more free parameters. This is a standard tool in maximum likelihood analysis allowing the comparison of non-nested models (Vaughan, 2005), that is, models in which parameter values are not a subset of those of another model. The best PSD background corresponds to the model that minimizes the AIC. Anderson et al. (1990) defined a fit acceptable when a MERIT value, defined as the ratio between the weighted sum of squared errors and the number of degrees of freedom (difference between the number of points and the number of the model free parameters), was lower than 1. For the adaptive MTM PSD, the MERIT value is where N S is the number of PSD values considered. We use the adaptive MTM expected value and variance (Thomson & Haley, 2014), that are, respectively, E{S j } = B j and The MERIT value represents the goodness of fit for least squares problems (Bevington & Robinson, 2003), but in our case, since the distribution of our data differs from a Gaussian distribution, it represents only a comparison tool. As with the AIC, the lower the MERIT value is, the better the representation is of the PSD background.
A Kolmogorov-Smirnov test can be applied to the ratio between the adaptive MTM PSD and the background model (see equation (5)). After evaluating D KS , the significance level can be approximately estimated by (Press et al., 2007): A confidence level for the fit is defined as C KS = 1 − P(D > D KS ) so that, in a similar way to the previous approaches, the minimum value identifies the best PSD background representation. The performance of the three criteria is discussed in section 5.2.

Selection of Signals and PSD Reshaping
Once the PSD background has been identified, we provide four options for signal identification. In the first option, we identify every portion of the PSD above a threshold defined as the product of the PSD background and the value z ϵ obtained by equation (9) at the confidence level ϵ. We refer to this procedure as the γ test. Among the frequency intervals corresponding to the PSD portions passing the γ test, we select only those whose width is greater than W, the half-bandwidth of the MTM spectral window. In the MTM approach, spurious peaks exhibit a triangular shape, while enhancements due to real periodicities show a rectangular shape of width ≈ 2W (Thomson & Haley, 2014). Due to the distortion of the enhancement's shape caused by noise, especially for low signal-to-noise ratio, a lower limit of W on the width is more appropriate. No upper limit is imposed in order to include the possibility of broad enhancements related to the occurrence of multiple signals at close frequencies (Di Matteo & Villante, 2017) or quasiperiodic signals whose frequency varies in time. For each portion of the PSD that passes the test, we identify the central frequency and the half-width. In the second option, we provide all of the harmonic F test local maxima above the defined confidence level (see section 2). We combine the results of the two tests in the third option: the selected frequencies are the ones identified in the harmonic F test that are within a PSD enhancement passing the γ test. Finally, for the last option, we impose the more stringent criterion allowing only the harmonic F test absolute maximum within each PSD frequency band selected by the γ test.
When strong periodic fluctuations, with a signal-to-noise ratio of several units or more, occur in the time series, the PSD background level provided by our procedure can increase above the true value and/or be distorted. In the case of narrow-band PSD peaks, the smoothing step reduces this effect, but better results DI MATTEO ET AL.
10.1029/2020JA028748 9 of 23 can be obtained by reshaping the PSD (Percival & Walden, 1993;Thomson, 1982). In our code, we implement an option to remove from the adaptive PSD estimate the contribution of strong signals identified by the combined γ and F tests at a given confidence level. Once these spectral peaks are removed, we apply our procedure again starting from this reshaped PSD. In the case of strong broadband PSD enhancements, where these expedients are ineffective, we suggest applying our procedure to a portion of the PSD unaffected by the strong signal to recover the global PSD background. Another solution is to implement a new PSD background model and/or the inclusion of a proper parameterization of the features of interest (e.g., power law-plus-constant combined with a Gaussian component; Inglis et al., 2015). This task is relatively simple given the modular structure of our code. In this scenario, the PSD model will provide information on the PSD background, the parameterized features, and the possible additional signals.

Examples with Synthetic Data
We discuss the performance of our procedure using synthetic time series representing lag-one autoregressive, power law, and bending power law processes. There are many methods to generate synthetic data with a specific PSD shape (Anderson et al., 1990;Timmer & Koenig, 1995;Vaughan et al., 2011). We use the approach of Timmer and Koenig (1995). Briefly, the square root of half the desired PSD is multiplied for two different series of Gaussian distributed random numbers. These vectors constitute the real and imaginary parts of a complex variable that, when extended with its complex conjugate, retrieve the double-sided Fourier transform of the desired data such that the synthetic data are obtained as its inverse Fourier transform. We generate synthetic time series in which we vary the N lengths, but we hold Δt = 1 s. We set the parameters: (i) c = 1 a.u. 2 Hz −1 and ρ = 0.90 for the AR(1); (ii) c = 1 a.u. 2 Hz β−1 and β = 1.5 for the power law; (iii) c = 1 a.u. 2 Hz β−1 , β = 1, γ = 3, and f b = 0.25 Hz for the bending power law. In the following, the PSD are evaluated using NW = 3 and K = 5 tapers.
The following discussion covers common scenarios for time series frequently observed in space physics environments, but it is not exhaustive. We note that the choice of the analysis parameters such as the time series length, N, the time-half-bandwidth product, NW, the number of tapers, K, the width of the smoothing window, and the parameter space of the models might influence the results. Therefore, we always recommend a preliminary investigation on specific sets of measurements to determine the best parameters for a robust spectral analysis.

Smoothing
The primary purpose of the smoothing procedure is to reduce the fluctuations of the estimated PSD around the true value in order to recover the shape of the PSD background, even when enhancements due to periodic signals may be present. Figure 1 shows results of a Monte Carlo simulation of 10 4 repetitions of time series with N = 512 points. From the left, each column shows results for the AR(1), PL, and BPL processes. From the top, we report an example of time series and the average of the raw and the med, mlog, bin, and but smoothed PSD (black thick line) with the 90% percentiles bounds (black thin lines). The red lines are the true PSD used to generate the synthetic time series. The smoothing windows have been identified automatically with the KS test. The distribution of the values p is available in the supporting information. Each procedure provides a different background approximation, primarily due to the different behavior at the edges of the frequency interval. Note that, compared to the raw PSD, all the smoothing procedures significantly reduce the 90% percentiles bounds.
The average raw PSD shows at low and high frequencies the effect of the convolution of the spectral window with the true PSD. At low frequencies, it underestimates the PSD for an AR(1) process at f j < W. On the other hand, for the PL and BPL processes, the average raw PSD flattens for f j < 2W and results in a spurious peak at f j ≈ W with respect to the true PSD. In all three processes, the true PSD is underestimated for f j > f Ny − 2W. The average med smoothed PSD provides a good representation of the PSD except at low frequencies, where it flattens due to the rapid rise of power. Therefore, it systematically underestimates the true PSD. The average mlog smoothed PSD, on the other hand, follows exactly the raw PSD at lower frequencies and flattens at high frequencies. This approach is particularly well-suited for AR(1) processes with ρ values that lead to a flattening toward a white noise PSD at high frequencies. In contrast, for PL and BPL processes, the average mlog smoothed PSD overestimates the PSD background. The bin smoothed PSD, known to be an unbiased estimator of power law PSD (Papadakis & Lawrence, 1993), gives a good representation of the PSD background in all three cases. Unfortunately, since the corresponding frequency range is reduced due to the binning, the values at low and high frequencies are extrapolated. Unlike the other smoothing procedures, the bin smoothed PSD is unaffected by the flattening at low frequencies, even though it slightly underestimates the PSD for the AR(1) process, and overestimates the PSD for the PL and BPL processes. The average but smoothed PSD is similar to the med one with a better representation of the true PSD at low frequencies for the AR (1)

(a) (b) (c)
high frequencies. Decreasing the percentage of smoothing, that is, increasing the pass band of the low pass filter, the but smoothed PSD could give a better representation of the PSD at low frequencies, but it would rapidly reduce to the raw PSD.

Background Estimate
Once the raw and/or the smoothed PSDs have been evaluated, the next step is to fit the different models with the maximum likelihood method. For each of the three simulated cases (AR(1), PL, and BPL), Figure 1 shows the average of the PSD backgrounds (green dashed lines) estimated from the raw PSD and its four different smoothed versions, while the red lines represent the true PSD. All of the PSD background estimation techniques provide a good representation of the true PSD for some portion of the frequency interval, but there are some differences. For the AR(1) process, the true PSD is underestimated at low frequencies by the med, bin, and but smoothed PSD, showing a clear flattening. At high frequencies, the true PSD is overestimated using the raw PSD and underestimated using the mlog, bin, and but smoothed PSD. For the PL case, the fits tend to overestimate the true PSD at all frequencies with the exception of the med approach, which lies below the true PSD at low frequencies. For the BPL case, instead, the true PSD is underestimated at low frequencies, while it tends to be overestimated at mid frequencies near the frequency break. At high frequencies, the true PSD is slightly overestimated using the raw PSD and underestimated using the mlog, bin, and but smoothed PSD.
Starting from the same synthetic time series used in Figure 1, in Figure 2 we show the distribution of the model parameters estimated via each smoothing + model combination. Since there are no signals in the simulated data, we expect the raw PSD to give the best results. However, this analysis helps to understand the biases that each smoothing procedure might introduce. For the AR(1) time series, the raw/AR(1) combination provides unbiased estimates for both the c and ρ parameters as the length N of the time series increases, as expected. The constant factor c, estimated using the smoothed PSDs, deviates from the true value as N increases, while the bin and but approaches provide a good estimate of the lag-one autocorrelation coefficient. For the PL time series, the raw PSD provides good estimates only for the slope β. The mlog, bin, and but smoothed PSD correctly estimate the constant factor c. We also note that the spread of the c distribution remains almost constant for time series longer than 1,024 points. For the BPL time series, both the raw/BPL and bin/BPL combinations provide unbiased estimates of all the model parameters. In addition, the mlog smoothed PSD determines a good approximation of the constant factor c and the slope β. The frequency break f b shows an uncertainty of about half of the entire frequency range even at the longest probed time series (2,048 points). However, the identification of the bending is strictly related to the gap between spectral indices. For large gaps, the uncertainty of the frequency break will be significantly reduced.
Finally, we investigated the performance of the AIC, MERIT, and CKS criteria in the selection of the PSD background. We report in Table 2 the rate of identification of a model by each criteria in the case of AR(1), PL, and BPL processes. Note that the BPL model, due to its flexibility, is able to approximate both the AR(1) and PL PSD retaining a relevant selection rate even in these cases. The MERIT criterion provides the best performance in all three of the scenarios with the correct selection rate above ≈ 57%. In addition, it properly excludes PL for AR(1) time series (false positive 0.22%), AR(1) for PL time series (false positive 2.53%), and both AR(1) and PL for BPL time series (false positive, respectively, 0.12% and 0.22%). The CKS criterion shows almost uniform rates with values between ≈ 30% and ≈ 37%. However, the maximum rates occur at the correct models, reaching ≈ 56% for BPL time series. The AIC criterion shows results similar to the MER-IT criterion for PL and BPL time series, but it completely excludes the AR(1) model. For AR(1) time series, the AIC criterion almost always selects the BPL model. Following these considerations, the MERIT criterion is the default in our code, while the values of the other criteria are provided.

Signal Identification
The last step is the identification of periodic signals according to the γ test, harmonic F test, and their combination. In the following discussion we do not impose constraints on the minimum width of the frequency intervals corresponding to the PSD portions selected by the γ test. Therefore, the frequency distributions shown here represent an upper bound for those obtained from signals identified by the γ test and its combination with the F test when imposing a minimum frequency width on PSD enhancements.

(c)
First, we evaluated the performance of our tests on the same synthetic time series used in the previous section. We evaluated the occurrence rate of false positives when imposing confidence thresholds at the 90% level. Figure 3 shows the distribution of false positives for AR(1) (panel a), PL (panel b), and BPL (panel c) time series. Each row corresponds to the results obtained using a specific smoothing procedure for the identification of the PSD background. For the AR(1) time series, the occurrence of false positives identified by the harmonic F test (green lines) is around the 10% level, as expected for a 90% confidence threshold, except at the edges of the frequency interval f j < 2W and f j > f NY − 2W (limits identified by the red vertical lines) that show higher rates. The γ test (blue lines) shows significantly fewer false positives. Unlike the F test, where the maximum F value above the confidence level is a single isolated outlier, the MTM windowing creates a group of consecutive outliers for periodic signals. For the γ test, we assign a single central frequency to this entire group. If we consider the distribution of every outlier according to the γ test, we would obtain a level of ≈10%. For comparison, we show the 10% level divided by the average width of the PSD enhancements above the confidence threshold (horizontal black lines). We found good agreement with the raw, mlog, and but approaches except for frequencies below 6W (left black vertical line) and above f NY − 2W where the occurrence rate decreases. For the med and bin approaches the distribution of false positives increases toward lower frequencies with a peak at f j ≲ 2W. The γ + F test (black lines) and γ + maximum F test (red lines) results are similar, but produce a flatter distribution than the single γ test with the number of false positives almost halved. For the PL time series, the distribution of false positives identified by the harmonic F test is similar to the AR(1) case with an additional small decrease toward low frequencies in the interval 2W < f j < 6W. The γ test determines a peak in the distributions at frequencies lower than 2W for all the approaches. The distributions obtained via the raw and mlog smoothed PSD provide good agreements with the expected level of false positives for the rest of the frequency interval. For the bin and but approaches, the distribution shows a slight increase toward high frequencies, while for the med approach, the distribution far exceeds the expected level in the first half of the frequency interval. The γ + F test and γ + maximum F test produce flatter distributions with the absence of the peak at low frequencies, except for the med smoothing results. For the BPL, the distribution of false positives identified by the harmonic F test is similar to the PL case. The γ test finds a peak in the distributions at frequencies lower than 2W, less pronounced for the raw and the mlog approaches. Other than with the raw PSD, the distribution of false positives manifests local enhancements between ≈0.2 and ≈0. We also estimated the rate of identification at the 90% confidence level of a monochromatic signal with frequencies spanning the entire frequency range. Each dot, at a specific frequency, in panel d-f of Figure 3, represents the results for 10 4 repetitions of a synthetic time series plus a signal at the corresponding frequency with signal-to-noise ratio equal to 0.8. Given the frequency of the signal f 0 , we considered the ratio of the power of a monochromatic signal A 2 /2, where A is the amplitude of the sinusoid, and the noise level, estimated integrating the theoretical PSD generating the time series over the interval Bold numbers indicate the highest rate in each scenario. show the identification rate at the correct frequencies of a monochromatic signal with signal-to-noise ratio equal to 0.8 and frequencies spanning the entire frequency range. The red (black) vertical lines correspond to the width (three times the width) of the main lobe of the multitaper spectral window, that is 2W ≈ 0.012 Hz (6W ≈ 0.035 Hz), from the limits of the frequency interval. PSD, power spectral density; AR, autoregressive process; PL, power law; BPL, bending power law. and models, with the exception of a jump to higher values for f < 2W. The γ test (blue dots) determines occurrence rates of ≈ 35-40% in all scenarios, except for f ≲ 6W where we observe a decrease for lower frequencies.
Only the but/BPL combination exhibits an opposite behavior. The γ + F test (black dots) and γ + maximum F test (red dots), in contrast with the analysis of false positives, find higher occurrence rates with respect to the single γ test. The reason is that for the latter, the chance of identifying nearby frequencies is high, with a probability to select f 0 − f Ray and f 0 + f Ray of ≈ 20% each, but the PSD enhancement contains the correct frequencies, f 0 , that are then identified by the harmonic F test, which has an almost null rate of false positives at nearby frequencies (Di Matteo & Villante, 2017). For the AR(1) time series (panel d), the rate of true positives is ≈60% for f > 6W and slowly decreases for f < 6W except for the results obtained via the med approach, which determines rates of up to 70%. For the PL time series (panel e), the identification rate is between 50% and 55% except for f < 2W, where values up to 75% occur, and for the med related results, showing an almost linear decrease from 70% to 45%. For the BPL time series, the occurrence rates of true positives closely follow the shape of the distribution of false positives (panel c) with values ranging between 40% and 65% for f > 2W. Note that the raw PSD and the bin approaches provide almost flat distributions.
The results provide insight into the performance of our method for the types of spectra commonly observed in geophysical environments. Depending on the circumstances, each smoothing + model combination provides good results in specific frequency ranges. In addition, we can use the biases quantified with these simulations to make a reasonable conclusion about the physics of the actual system.

Demonstration of the Technique Applied to Observations
In this section, we apply our approach to a case study using data taken in the solar wind, magnetosphere, and ground observations to demonstrate the performance of our methodology. We consider the periodic fluctuations previously identified by Viall et al. (2009) in the solar wind proton density and magnetic field measurements at geostationary orbit on January 15, 1997, extending the analysis to a longer time interval and to ground observatories. Viall et al. (2009) showed that the magnetospheric field fluctuations observed by GOES9 can be a consequence of the quasi-static modulation of the magnetospheric cavity size by the solar wind dynamic pressure in turn related to the solar wind density variations. However, the frequency of these oscillations is in the same range of magnetospheric ULF waves that can be triggered by numerous processes, including solar wind pressure pulses, flow shear instabilities at the magnetopause, and wave-particle interactions in the inner magnetosphere. The ability to identify these waves is the first step in distinguishing between the different possible formation mechanisms and in furthering our understanding of them. Identification of these waves, especially in the case of low signal-to-noise ratio, is often affected by the limitation of the adopted spectral analysis techniques. Here, we show that even though the PSD background in the solar wind, magnetosphere, and ground observations exhibit considerably different shapes, our technique exhibits great flexibility and is able to provide good background estimates and identify a common periodicity among all of the PSDs.

Solar Wind
Periodic structures in the solar wind proton density were observed by the Wind spacecraft on January 15, 1997, between 12:40 and 19:10 UT. We used proton density data derived from the Wind-Solar Wind Experiment (Ogilvie et al., 1995) measurements. The time interval of 6.5 h determines a Rayleigh frequency of f Ray ≈ 43 µHz, while the average sampling rate of Δt ≈ 83 s corresponds to a Nyquist frequency of f Ny ≈ 6 mHz. We choose NW = 3 and K = 5 as parameters for the MTM analysis, therefore the bandwidth for the spectral window is 2W/Δt ≈ 0.26 mHz corresponding to the minimum separation needed to distinguish two signals with close frequencies. Figure 4a shows the proton density observations, n p , while Figure 4b shows the corresponding PSD, the γ and harmonic F tests. Applying our spectral analysis procedure, the best fit PSD background identified (red line) is the bin/PL pair with parameters c ≈ 0.033 cm −6 Hz β−1 and β ≈ 1.39. Then, we tested the occurrence of periodic signals at the 90% confidence levels (red dashed lines). We placed circles above the PSD enhancements passing the γ test and crosses at the frequencies that also passed the harmonic F test within the same frequency range. We identified three clear signals passing both tests at ≈ 0.88, 2.25, and 3.89 mHz corresponding, respectively, to ≈19, 7.4, and 4.3 min. An additional periodicity at f ≈ 0.17 mHz (≈100 min) was identified only by the γ test. Viall et al. (2009), using the Mann and Lees (1996) DI MATTEO ET AL.

10.1029/2020JA028748
approach over part of the same time interval analyzed here, identified periodic fluctuations passing both the narrow band and the harmonic F test at f ≈ 0.2, 0.8 and 2.8 mHz.

Magnetosphere
The solar wind described in the previous section was measured near L1, and impacted the magnetosphere ≈45 min later, corresponding to the time range from 13:25 to 19:55 UT. We investigated the magnetospheric response considering the 60 s (f Ny ≈ 8.3 mHz) averaged magnetic field components derived from the triaxial  (Singer et al., 1996) on the GOES9 geostationary satellite (LT = UT-9) located in the dawn-morning sector (between 4:25 and 10:55 LT). The data have been rotated in the Mean Field Aligned (MFA) coordinate system at each point along the spacecraft trajectory. In MFA coordinates (Takahashi et al., 1990),  is along the average field, as defined by a vector running average;  is perpendicular to  and the spacecraft position vector, positive eastward; ˆ completes the orthogonal system. To avoid the introduction of spurious periodicity due to the rotation procedure, the average magnetic field is evaluated on a running window of 6.5 h (Di Matteo & Villante, 2018). Figure 4c shows the three components of the magnetospheric field, while Figure 4d shows the corresponding PSDs, the γ and harmonic F tests. The similarity of the compressive component B μ with the solar wind density fluctuations is clear, even though the higher frequencies component seems to be filtered out in the magnetosphere at the GOES9 location. Next, we investigate the occurrence and properties of the magnetospheric field fluctuations with our spectral analysis approach. Our method selects the raw/BPL PSD background with parameters c ≈ 0.72 nT 2 Hz β−1 , β ≈ 1.17, γ ≈ 3.57, and f b ≈ 0.33 mHz for the compressive component (B μ ), the raw/PL PSD background with parameters c ≈ 1.05 × 10 −6 nT 2 Hz β−1 and β ≈ 2.31 for the toroidal component (B φ ), and the raw/BPL PSD background with parameters c ≈ 6.18 × 10 −3 nT 2 Hz β−1 , β ≈ 1.56, γ ≈ 3.15, and f b ≈ 0.36 mHz for the poloidal component (B ν  an enhancement at f ≈ 0.9 mHz (≈20 min) clearly observed also in the solar wind proton density. In the toroidal component, the signal at f ≈ 0.43 mHz, corresponding to oscillations of about ≈39 min, is probably related to the first three oscillations observed at the beginning of the time interval (≈26, 32, and 36 min). Similar fluctuations appear also in the solar wind proton density (≈26, 40, and 36 min), even though there is no clear enhancement in the PSD. In fact, other stronger fluctuations at nearby frequencies dominate the low-frequency range of the solar wind density PSD making it difficult to distinguish additional quasiperiodic signals.

Ground Observatories
We extended the analysis to ground magnetic field observations from two stations located near the GOES9 magnetic field line footpoint: Yellowknife (YKC, λ = 62.48° and ϕ = 245.52°) and Fort McMurray (FMC, λ = 56.66° and ϕ = 248.79°), where λ and ϕ are the geographic latitude and longitude, respectively. For these examples, we used the 60 s data from the SuperMAG collaboration providing the three components of the magnetic field in the NEZ coordinate system where B N and B E are directed toward the locally magnetic north and east, respectively, and B Z is vertically down. We analyzed the B N component after the removal of the daily variations and yearly trend determined by the Gjerloev (2012) algorithm. Figure 4e shows the magnetic field observations from the two stations, while Figure 4f shows the corresponding PSD, the γ and harmonic F tests. Applying our procedure, we obtain the raw/BPL PSD background with c ≈ 0.05 nT 2 Hz β−1 , β ≈ 1.82, γ ≈ 9.25, and f b ≈ 6.64 mHz at YKC, and the raw/PL PSD background with c ≈ 0.02 nT 2 Hz β−1 and β ≈ 1.57 at FMC. As in the previous section, we classified the signals identified at the 90% confidence level. At YKC, we observed three PSD peaks at f ≈ 0.86, 4.88, and 5.14 mHz passing both the γ and the harmonic F test. The spectral analysis at FMC identified one signal satisfying both the γ and the harmonic F test at f ≈ 0.86 mHz and two signals at f ≈ 6.04 and 7.84 mHz selected only by the γ test. The two ground observatories observed magnetic field oscillations at f ≈ 0.9 mHz as in the magnetospheric field at geostationary orbit, in turn, driven by the solar wind density fluctuations. Viall et al. (2009), using the Mann and Lees (1996) approach, identified, during part of the same time interval, periodic fluctuations passing both a narrow band and harmonic F test, with the 95% confidence level, at f ≈ 0.2, 0.8, and 2.8 mHz in both the solar wind proton density and B z magnetospheric field component at the geostationary orbit. We find similarity with our results at f ≈ 0.17, 0.89, and 2.26 mHz in the solar wind and at f ≈ 0.9 mHz at the geostationary orbit and in the two ground observatories. The time series of B μ at GOES9 suggests that the longer timescales are directly driven by the solar wind density fluctuations.

Additional Remarks
In addition, we note that in the low-frequency range of the γ statistic, three enhancements centered at f ≈ 0.2, 0.4, and 0.9 mHz occurred in all the observations, but our procedure identified only the strongest component at f ≈ 0.9 mHz. The difficulty in the identification of PSD peaks at nearby frequencies and at the edges of the frequency interval are two known limitations of spectral analysis methods. The simplest way to overcome these issues is to increase the frequency resolution, either by increasing the length of the time interval or decreasing the width of the spectral window main lobe (reducing the NW parameter for the MTM). Another alternative is the analysis of overlapping time intervals to construct dynamic γ and F tests. While the γ test will always show PSD enhancements with a width equal to or greater than 2W, the F test might distinguish simultaneous signals at close frequencies, that is |f i − f j |≲ 2W, depending on the characteristics of the signal itself (Di Matteo & Villante, 2017). The occurrence of multiple signals might be revealed by the distribution of the frequencies identified with the F test in each patch of the dynamic γ test above the confidence thresholds. This approach can also be extremely useful when the signal frequency changes in time. A possible improvement to our procedure, especially when facing multiple signals, is the implementation of the approach developed by Denison et al. (1999), who provide an alternative significance test to the simple harmonic F test when facing time series with embedded signals at close frequencies. Another alternative is the extension of our approach to multivariate spectral analysis (Walden, 2000), simultaneously analyzing the time series of interest.

Discussion
We presented a new spectral analysis procedure, based on the adaptive MTM method, for the robust modeling of the PSD background and identification of signals at discrete frequencies. The adaptive MTM was specifically introduced by Thomson (1982) to investigate colored PSD when common spectral analysis techniques might suffer from strong energy leakage, especially for short time series. One major challenge in the analysis of the PSD of space physics time series is their wide range of variability. In general, even when the physical process at work in the creation of the PSD is well known, any individual instances may not produce a fully developed PSD of that type, so the flexibility provided by our algorithm may still be needed. For example, the PSD spectral slope of the solar wind parameters in the inertial range evolves with increasing distance from the Sun, steepening from −3/2 to −5/3 for the velocity (Roberts, 2010) and magnetic field (Chen et al., 2020), or may tend toward −2 in the presence of discontinuities (Roberts, 2010) or anisotropies (Horbury et al., 2012). We use the statistical properties of the adaptive MTM to develop a maximum likelihood determination of the PSD background as in Vaughan (2010). In addition, we extended the Mann and Lees (1996) approach, combining different smoothing methods (raw, med, mlog, bin, and but) and spectral models (WHT, PL, AR(1), BPL). Finally, we defined objective criteria to select the best representation of the PSD background and to identify spectral peaks in the PSD and F values at defined confidence levels.
We examined the characteristic features of PSD background identification via Monte Carlo simulations of synthetic time series representing lag-one autoregressive, power law, and bending power law processes. The first step is the smoothing of the raw PSD, useful when large PSD enhancements due to geophysical periodic signals are present. The user can choose from four different raw PSD smoothing approaches, each of which has its own advantages and disadvantages for fitting colored PSDs. The med approach systematically underestimates steep PSD at low frequencies on an interval comparable to the width of the running window. However, it might give a better representation of the PSD background when strong clear peaks occur at very low frequencies. The mlog approach, instead, reproduces the raw PSD at low frequencies, while at high frequencies, due to the running window covering a large portion of the frequency interval, returns almost constant values. This behavior is optimal for a AR(1) process, when the PSD flattens at high frequencies, but is not well-suited for very steep PSD where it will overestimate the PSD background. The bin approach defines the smoothed PSD on a limited range of frequencies, therefore the PSD background at the edge of the frequency interval is extrapolated. However, this procedure provides a good representation of the PSD background in all three processes studied. The but approach provides results similar to med but with better estimates in the low-frequency range.
When we fit the different models to the smoothed PSDs, we obtain a good representation of the true PSD for the majority of the smoothing + model combinations in the interval 2W < f j < f Ny − 2W. In absence of signals, the use of the raw PSD ensures good results in all of the scenarios, as expected.
In the examples with synthetic time series, we show that for steep PSDs, especially for power law processes, the low frequencies portion, which is not well represented even by the raw PSD for f j < 2W, plays a critical role in the identification of a reliable background model. This is mainly a concern for short time series that might have few points in the low-frequency range. For bending power laws, additional complications might arise when the frequency break is too close to the edges of the frequency interval or when the two spectral indices are similar; in these scenarios the BPL can be easily mistaken for a PL. Therefore, a necessary condition for the BPL is to have enough points in each of the two frequency intervals that exhibit different spectral slopes. Both problems might be resolved by considering time series long enough to ensure adequate coverage for both regimes of the PSD. When there is a lack of information about the properties of the background model, our technique allows for the smoothing + model combinations to be calculated and the best representation selected according to objective statistical criteria. This is particularly helpful when PSD enhancements due to periodic fluctuations are present.
Using synthetic time series, we found that fitting the chosen model to the raw PSD provided the best results, except for a constant factor offset for PL time series. For the AR(1) process, ρ is best estimated with bin and but; for the PL process c can be estimated with mlog, bin, and but, and β with mlog and bin. For the BPL process c and β are best estimated with mlog and bin, while γ and f b via bin. Overall, among the smoothed PSD, we obtain the best performance with the bin approach followed in order by the mlog, but, and med approaches.
Regarding the identification of periodic signals, the distribution of false positives estimated with the Monte Carlo simulations agrees with the expected rate over the frequency interval 6W < f j < f Ny − 2W. In particular, we find that, in this frequency range, the identification rates of true signals are flat for the raw PSD, as expected. In addition, we obtain flat distributions via the mlog, bin, and but approaches for AR(1) and PL time series, and via the bin approach for BPL time series. Outside this interval, the false positive rate can significantly differ from this expected rate, and care should be taken when testing in this range.
We demonstrated our technique by analyzing observations of solar wind proton density, magnetospheric field at geostationary orbit, and magnetic field at two ground stations. We considered a previously studied time interval during which the solar wind density directly drove compressional fluctuations in the magnetospheric field at geostationary orbit and the magnetic northward component at ground observatories (Viall et al., 2009). The best PSD background representation identified by our procedure corresponded to a power law for Wind measurements of n p , B φ at GOES9, and for B N at FMC, and a bending power law for B μ and B ν at GOES9; and B N at YKC. AR(1) was not found to provide the best fit background model for any of the data. This demonstrates the need for utilizing different models for a correct evaluation of the PSD background, especially in cases like the YKC observatory where only the BPL provided reasonable results.

Conclusions
We have developed an automated method for identifying both the background and significant enhancements of PSDs. We start with the adaptive MTM, a sophisticated nonparametric spectral analysis tool suitable for the analysis of colored PSD. The knowledge of the statistical properties of the PSD allows a robust maximum likelihood fitting of four models on the raw PSD and four smoothed PSDs. The best representation of the PSD background, selected via a robust statistical criterion, determines the confidence thresholds used to identify statistically significant PSD enhancements and, when combined with the harmonic F test, robustly identifies the frequency of the periodic oscillations occurring in the time series.
The Monte Carlo simulations of synthetic time series demonstrates how different combinations of smoothings and models influence the determination of the PSD background, and hence the confidence levels of the PSD enhancements. Our method is not meant to be a black box to be applied to any time series, but rather a useful tool providing different paths from which a user can choose the best combinations for the data being analyzed. The Monte Carlo simulations of synthetic time series show clearly that not all paths provide good results. We highlight that a preliminary analysis on the data of interest is the best practice to assure a robust application of our method. For the specific case analyzed in our simulations, we can conclude that the recommended smoothing are bin and but for AR(1) time series, mlog and bin for PL time series, and bin for BPL time series. Note that the other smoothing approaches provide good results in a narrower frequency range.
We also demonstrated the inherent flexibility of our method by applying the analysis to real measurements in three different geophysical environments for the same event.
The approach developed here can be extended to a broad range of disciplines that need to distinguish between continuous PSD and discrete PSD enhancements. Such applications range from analyzing time series for statistically significant periodicities to robustly characterizing the PSD background. The present work also lays the foundation of a Bayesian approach for estimating the posterior distribution of the PSD model parameters using the MTM PSD. The modular structure of our methodology allows the introduction of new smoothing methods and models to cover additional types of time series. The flexibility and extensibility of the technique makes it broadly suitable to any discipline. Generally speaking, this technique provides a good representation of the PSD background thanks to the different smoothing + model pairs covering more scenarios than previous spectral analysis methods. When combined with an independent harmonic analysis, this allows the robust identification of PSD enhancements related to monochromatic fluctuations occurring in the time series.

Data Availability Statement
The solar wind observations from Wind and the GOES magnetic field data are available at the NASA's CDAWeb site (http://cdaweb.gsfc.nasa.gov/istppublic/). Ground observations rely on data collected at magnetic observatories and are available at http://supermag.jhuapl.edu. Code is available on the Zenodo platform at https://zenodo.org/record/3703168.