Articles | Volume 17, issue 2
Research article
31 Jan 2024
Research article |  | 31 Jan 2024

Assessing atmospheric gravity wave spectra in the presence of observational gaps

Mohamed Mossad, Irina Strelnikova, Robin Wing, and Gerd Baumgarten

We present a thorough investigation into the accuracy and reliability of gravity wave (GW) spectral estimation methods when dealing with observational gaps. GWs have a significant impact on atmospheric dynamics, exerting influence over weather and climate patterns. However, empirical atmospheric measurements often suffer from data gaps caused by various factors, leading to biased estimations of the spectral power-law exponent (slope) β. This exponent describes how the energy of GWs changes with frequency over a defined range of GW scales. In this study, we meticulously evaluate three commonly employed estimation methods: the fast Fourier transform (FFT), generalized Lomb–Scargle periodogram (GLS), and Haar structure function (HSF). We assess their performance using time series of synthetic observational data with varying levels of complexity, ranging from a signal with one frequency to a number of superposed sinusoids with randomly distributed wave parameters. By providing a comprehensive analysis of the advantages and limitations of these methods, our aim is to provide a valuable roadmap for selecting the most suitable approach for accurate estimations of β from sparse observational datasets.

1 Introduction

Gravity waves (GWs) are ubiquitous phenomena that play a crucial role in the dynamics of the Earth's atmosphere, where they impact weather and climate patterns (Hines1960; Ern et al.2018). Various sources, including convection, topography, and jet streams, generate these waves (Crowley and Williams1987; Fritts1989). As they propagate through the atmosphere, they can break and mix with the surrounding atmosphere, redistributing their energy and momentum. This leads to significant changes in the atmospheric thermodynamics and large-scale circulation patterns of the atmosphere, including wind speeds and temperature gradients (Lindzen1981; Holton1983; Fritts and Alexander2003). Observations of these meteorological variables reveal that GWs exist for the most part in the form of a spectrum of superposed waves within a wave packet and occasionally as quasi-monochromatic waves (Maekawa et al.1984; Eckermann and Hocking1989). To understand the physical processes that govern the generation, propagation, and dissipation of these wave packets, it is often useful to examine their spectral properties such as their frequencies, amplitudes, and scales (Axford1971; Fritts and VanZandt1993).

On that note, VanZandt first introduced the concept of a “universal atmospheric GW spectrum” (VanZandt1982). This spectrum facilitated efficient parameterizations of how GWs affect the mean atmospheric state (Narendra Babu et al.2008). For instance, the spectra of GWs are often used in model parameterization, including source spectra parameterization, Lagrangian spectral parameterization, and subgrid-scale parameterization, enabling the simulation of the dynamics of the middle and upper atmospheres (Beres et al.2005; Song and Chun2008; Houchi et al.2010). Overall, accurate predictions of GW activity can improve weather forecasting, while they contribute significantly to climate modelling in parameterizing physical processes like turbulence and mixing (Alexander et al.2002; Smith2012; Liu et al.2014).

This GW spectrum exhibits a power law scaled by an exponent (or slope) β, which describes the rate at which wave energy changes with its wavenumber (or frequency). The basis for this spectrum for atmospheric GWs is supported not only by a strong foundation in theoretical works (Dewan and Good1986; Weinstock1990; Hines1991; Dewan1994; Gardner1994), but also in observational studies (Smith et al.1987; Fritts et al.1988; Gardner et al.1995; Nastrom et al.1997; Zhang et al.2006, 2017). These values of β depend on not only the type of spectra (e.g. temporal or horizontal, or vertical wavenumber) but also the geophysical variables measured (e.g. temperature, horizontal or vertical wind); see Table A1 for a summary. Thus, an accurate estimation of β is essential to validate different theoretical predictions of GW power spectral densities (PSDs) (Dewan and Grossbard2000) and improve climate models and weather forecasts (Lindgren et al.2020).

Determining β from empirical atmospheric measurements is challenging due to various factors, such as the inevitable presence of data gaps, observational noise, and the finiteness of data length and resolution. Data gaps can occur for numerous reasons, including instrumental errors, data transmission issues (e.g. due to weather conditions like clouds in the case of lidar), and signal interference (in the case of radar). When gaps exist in multiscale time series, data points representing certain frequencies are lost, which distorts the spectra and introduces significant bias into the estimation of β (Brown and Christensen-Dalsgaard1990; Rigling2012). To minimize the effect of these gaps on the spectra, data-filling schemes are often applied. Though linear interpolation is usually used to fill in these gaps (Meisel1978; Lepot et al.2017), even adaptively implemented interpolators produce artefacts in the time series at low gap percentages (GPs), which contribute additional bias in the spectra (Schulz and Stattegger1997; Hall and Aso1999). Bias in spectral estimates can also be caused by other relevant sources, such as spectral leakage, steep spectra (β>2), and in-signal components with larger periods than the observed time span T (Klis1994).

In this paper, we systematically quantify the advantages and limitations of estimation methods of GW spectra in handling these error sources. We also propose a procedure for selecting unambiguously suitable methods for β estimation. Even though this study is motivated by the analysis of atmospheric GWs, the conclusions and the methods can be generalized to different fields with similar time-series characteristics in other branches of geophysics. Two commonly used methods are considered, namely the fast Fourier transform (FFT) (Cooley and Tukey1965) and the generalized Lomb–Scargle periodogram (GLS) (Zechmeister and Kürster2009), as well as the fairly recent Haar structure function (HSF) (Lovejoy and Schertzer2012). FFT is the standard method to analyse spectra of evenly sampled data. The Lomb–Scargle periodogram (LS) has been used in many studies as the main analysis method (or as a reference) of GW spectra (e.g. Hall and Aso1999; Zhang et al.2006; Guharay and Sekar2011; Qing et al.2014). As far as the authors are aware, HSF has never been used in atmospheric GW studies. Both GLS and HSF are specifically known to handle unevenly sampled data. In an effort to closely mimic real observations of GWs, we simulate time-series data with varied levels of complexity, beginning with a signal with one frequency and increasing in complexity to a superposition of sinusoids with randomly distributed frequencies (and phases).

Previous studies have investigated these spectral methods and others for estimating power-law spectra and compared their performance using synthetic and observed data. For instance, Zhan et al. (1996) found that using FFT of linearly interpolated signals is the best approach to analyse radar wind data at 50 % GP (only for the case of β=5/3) compared to the correlogram and LS. However, a quantitative analysis of the effect of changing β or the GP was not conducted. Similarly, Munteanu et al. (2016) showed that FFT outperforms LS, Z transform, and discrete Fourier transform in estimating β from Venus' magnetic field data, although the effect of changing β was not considered either, since the power-law spectra were not simulated. In contrast, Hébert et al. found that HSF consistently surpassed other methods in estimating β, without the need to interpolate the gapped (simulated palaeoclimate) data for β(0,3), except in the case of β(-1,0), where they concluded that LS would be the best option (Hébert et al.2021). Nonetheless, the impact of altering the GP was not quantitatively presented; instead, the skewness of the gaps' (gamma) distribution was used as a parameter to refer to the irregularity of the time series. These three studies showed that the LS method suffers from significant leakage in the case of power-law spectra, which persists even when multi-taper methods (MTMs) are used.

The rest of the paper is organized as follows: in Sect. 2, we describe the methods used in our study, including a description of FFT, GLS, and HSF. In Sect. 3 we introduce the data simulation procedures. In Sect. 4 we discuss data processing. In Sect. 5 we present the results of our simulations, comparing the performance of these methods in different scenarios. In Sect. 6 we discuss the implications of our findings and provide recommendations for spectral analysis of GW time series with data gaps. Finally, in Sect. 7, we present a summary of our relevant results and conclusions.

2 Spectral methods

2.1 The fast Fourier transform

FFT is the most commonly used method for estimating frequency spectra of evenly sampled data (Cooley and Tukey1965). It enables the approximation of a time series sampled from a continuous distribution over discrete time steps, through a series of complex sine and cosine waves with varying frequencies. Under the assumption of a unit sample interval, (forward) FFT converts a time series zn of length N from its original domain (time or space) into a set of coefficients Zk in the (temporal or spatial) frequency domain by employing the following relation:

(1) Z k = n = 0 N - 1 z n e - 2 π i k n / N , k = 0 , 1 , , N - 1 .

In our work, FFT will serve as the benchmark spectral estimation method. The expected Fourier transform of a discretized signal is given by the convolution of the true transform and the transform of a Dirac comb window function designating those measurement times (Vanderplas2018). In the case of gapped data, the symmetry in the Dirac comb is destroyed, causing the resulting transform to be noisy with incorrect peak positions and heights. Consequently, the true transform of gapped data will not be recoverable. This disadvantage can be bypassed by applying data reconstruction methods such as interpolation, sparse approximation, etc., to approximate the true Fourier transform (Babu and Stoica2010). Unfortunately, these reconstruction methods can introduce artefacts into the signal, which depend on the distribution of the gaps and their sizes (Munteanu et al.2016).

2.2 The generalized Lomb–Scargle periodogram

The GLS periodogram developed by Zechmeister and Kürster (2009) offers a method for estimation of the PSD of unevenly sampled time series. It is a generalization of Lomb's least-squares approach (Lomb1976), which is equivalent to the modified Schuster's periodogram (Schuster1898; Scargle1982) (based on FFT) in the case of evenly sampled data. GLS produces a spectrum by least-squares fitting a model of a weighted sinusoid given by

(2) y ( t ) = a cos ω t + b sin ω t + c

to the time series at each sampled frequency ω. The offset c compensates for the assumption that the mean of the time series z is equal to the mean of the fit y. This floating-mean approach is advantageous, considering that the mean of a periodic signal may change statistically, especially for small N (Ferraz-Mello1981). Furthermore, the purpose of using weighted sums is to account for the observational noise for which the original LS does not.

The LS method has often been used to seek dominant periodic frequencies or cycles (Zhang et al.1993; Pichon et al.2015; Rao et al.2017), analyse seasonal changes in significant modulations of GW fields (Beldon and Mitchell2010), and estimate the spectral indices β and amplitudes of GW power-law spectra (Hall and Aso1999; Zhang et al.2006; Guharay and Sekar2011; Qing et al.2014). In addition, LS is known as the most efficient method for estimating the variance in both gapped and non-gapped stationary time series with a single periodicity, without the need to fill in missing data (Marinna et al.2019). In contrast, Vio et al. (2010) found that LS is reliable for analysing neither semi-periodic nor aperiodic signals with non-stationary noise or signals made of more than one wave, without additional steps.

2.3 The Haar structure function

HSF is a mathematical tool used in conducting scaling analysis of signals (Lovejoy and Schertzer2012), which is based on the Haar wavelet (Haar1910). It is a simple yet powerful method for decomposing a signal x(t) whose power spectral density exhibits a power law, i.e. PSD∝τH with a scaling (Hurst) exponent H, over a scale (lag) τ=1/f, into fluctuations Δx=x(t+τ)-x(t). The first-order Haar fluctuations τ at a lag τ are defined by the following relation:

(3) H τ ( x ( t ) ) = 2 τ t + τ 2 < t < t + τ x ( t ) - t < t < t + τ 2 x ( t ) .

The qth-order structure function Sq(τ) is then obtained as an approximation by ensemble-averaging these fluctuations as follows:

(4) S H , q ( τ ) = H τ ( x ( t ) ) τ q H - K ( q ) .

In the quasi-Gaussian case, the moment scaling function is K(q)≈0, so for the first-order structure function (q=1) only H determines the scaling of mean fluctuations.

A power spectrum which follows a power-law PSD∝τβ is related to the Hurst exponent by β=1+qH-K(q); since the power spectral density is a second-order moment, we take q=2. Thus, under our quasi-Gaussian approximation, we re-scaled HSF to a comparable scale to the PSD using the following relation:


Despite the term K(2) being fairly small, it is nontrivial in the non-Gaussian case and for higher moments q, and HSF allows for its calculation (Lovejoy and Schertzer2012), but this is beyond the scope of this paper. HSF is particularly suitable for estimating the scaling exponent of time series with H(-1,1) or β(-1,3). This range of β values covers the vast majority of atmospheric processes from weather (where τ<10 d and 1<β<3) to macroweather (where 10d<τ<10–30 years and -1<β<1) systems. HSF also possesses the advantage of handling unevenly sampled data, which is a consequence of the fact that it is computed by taking the mean of absolute fluctuations (Lovejoy2014). Nevertheless, HSF is not employed to estimate the frequency of a wave or its amplitude, since it only measures how much frequency components contribute to the total variance. The Python code implementation of HSF is readily accessible (Mossad2023) and was derived from the R code originally developed by Raphaël Hébert (Hébert2021).

3 Data simulation

In this section, we present the simulation procedures used to generate time series similar to actual GWs measurements. In measurements, GWs can exhibit various behaviours, ranging from superposed waves within wave packets with multiple frequencies, amplitudes, and phases to more coherent quasi-monochromatic waves (Maekawa et al.1984; Eckermann and Hocking1989; Sica and Russell1999). In Sect. 3.1, we simulate signals with each having one frequency to mimic quasi-monochromatic waves, by which the goal is to accurately estimate the correct frequency and amplitude of each signal. In Sect. 3.2, however, we simulate time series composed of a superposition of waves with random frequencies and power-law amplitudes. The exponents/slopes β of their power-law spectra are used to assess the bias in the spectral-analysis methods. While our simulation adopts a simplified linear saturation theory approach through the superposition of sine waves (Dewan and Good1986; Smith et al.1987), we acknowledge that other explanations for the spectral character of GWs exist, including “nonlinear damping” (Weinstock1982, 1990; Gardner1994), “Doppler spreading” (Hines1991), “saturated-cascade similitude” (Dewan1994), and high-Reynolds-number “stratified turbulence” (Pinel and Lovejoy2014); see Table A1 for more details.

By analysing both simulations, we can determine the accuracy of each of the methods at different levels of signal complexity and identify potential limitations and sources of error in the analysis of GWs spectra. Random gaps are then introduced to resemble observational gaps for both simulations. The units and values of the variables used in this simulation have been chosen to represent average values or ranges characteristic of typical GW time series.

3.1 Single-frequency signal simulation

Quasi-monochromatic GWs can be observed under specific conditions where a single frequency dominates other components (Muraoka et al.1988; Swenson et al.1999). This kind of GW can be approximated as an evenly sampled single sinusoid x(t)=Asin(2πft+φ), with a known frequency f, phase shift φ, and amplitude A, at time t. For both simulations, the time resolution Δt is 5 min with a total span of T=6 h, as this resolution and this duration align with the average values of lidar measurements commonly used in atmospheric studies (Gardner et al.1995; Gerding et al.2008). As a result, the number of points N in each signal is equal to TΔt=72. Each simulated sinusoid has an amplitude of A=4 K and a randomly chosen frequency f from the set 1/{6,3,1.5,1,0.5,1/3} [h−1]. Changing the frequency serves as a test to examine whether the bias in the methods is frequency-dependent. The phase shift φ is also randomly chosen from a uniform distribution within the interval [0,2π].

The amplitude of the time series is equivalent to the estimated height of the main peak in the spectrum. It is computed from the FFT coefficients (Eq. 1) as AFFT=maxk2ZkN and from the GLS fit coefficients (Eq. 2) as AGLS=maxkak2+bk2. The frequency of that peak fk corresponds to the estimated frequency of the signal. As a metric for the accuracy of estimation of the true values of frequencies and amplitudes, we used the relative bias given by

(6) relative bias = value estimated - value expected value expected .

Since real data are susceptible to observational noise, it is crucial to consider the case where white noise is added to the simulated signal as a random variable r(t) from a standard normal distribution. Here, the signal-to-noise ratio is defined by SNR=A2/2σr2, where σr2 is the noise variance (Horne and Baliunas1986). To strike an appropriate balance between capturing meaningful noise characteristics and minimizing low SNR bias, an average SNR value of 8 is chosen for this simulation.

3.2 Spectral power-law simulation

As reported before (see Table A1), the spectra of GWs are characterized by a power law; i.e. PSD1/fβ (in the case of spatial data, PSD1/(k,l,m)β, where k, l, and m are horizontal and vertical wavenumbers). On that account, we are interested in estimating the value of β of simulated time series whose spectra would have β{-1,0,1,5/3,2,2.5,3} and comparing it with the true value. The simulated evenly sampled time series x(t) consists of a sum of M sinusoids, each with frequency fi and power-law amplitudes Ai=fi-β, as follows (Rice1944; Billah and Shinozuka1990; Kirchner2005):

(7) x ( t ) i M A i sin ( 2 π f i t + φ i ) = i M ( f i ) - β / 2 sin ( 2 π ( f i ) t + φ i ) .

We used for our simulation M=35 for each time series. It is driven by the objective of reconciling the inclusion of a minimum of 20 waves based on observational (Sica and Russell1999) and modelling (Dewan1994; Hamilton1997) suggestions while simultaneously incorporating a sufficiently large number of waves to mitigate random spectral errors (Keisler and Rhyne1976; Shinozuka2005). This approach enables the demonstration of the power-law spectrum without the need for excessive averaging. The phase shifts φi are also randomly chosen from a uniform distribution within the interval [0,2π]. The frequencies fi are statistically independent, uniformly distributed random values, selected within the range 1T,12Δt=16h,110min. Hence, the time series is composed of non-harmonic components and there are no favoured frequencies, which is a better approximation of atmospheric GWs than an idealistic case where frequencies are only integer multiples of a fundamental frequency. Here, x(t) is proportional to the square root of the amplitudes Ai, since β is estimated from PSD-normalized spectra, which are the squared modulus of the amplitudes.

The PSD is obtained by FFT using the relation 2ZkZk*ΔtN, where Zk* is the complex conjugate of Zk. This definition is equivalent to the GLS spectrum Nak2+bk2Δt2 of evenly sampled data. HSF is however normalized according to Eq. (5c) to estimate a comparable scale to the PSD. The bias in β estimation is defined for this simulation as

(8) β bias = value estimated - value expected .

3.3 Gap simulation

After creating a time series with the desired spectral properties, gaps are introduced by randomly removing data points (except both endpoints), assuming that all data points are equally likely to be removed (i.e. a uniform distribution). Based on the simulated GP p in the data, an integer number of random points NG=Np100 is removed. Thus, a 0 % GP means that no points were removed, while a 50 % GP means that 36 points were randomly removed, since N=72. To assess the dependence of bias in spectral-analysis methods on the gaps, we conducted simulation runs spanning GPs ranging from 0 % to 90 % in increments of 10 % for each time series analysed. Each of these simulation runs was repeated 1000 times at each GP increment to ensure the statistical significance of our results, since the frequencies, phases, and gaps are randomized. Then we computed the average values of the estimated amplitude and period (for the single-frequency signal simulation) and β (for the spectral power-law simulation).

4 Processing steps

Before applying spectral methods to the generated time series from the simulations in Sect. 3, the following steps are taken:

  • The time series is first interpolated using the original time step of 5 min; this is only necessary for FFT.

  • The mean of the signal z is subtracted to account for the zero-frequency component of the Fourier transform Z0.

When computing the spectrum for a time step Δt, the frequency grids of all methods are defined as follows:

  • The frequency range spans from 1T to 12Δt.

  • The frequency spacing is given by Δf=1T.

In the presence of gaps where Δt is not constant, the Nyquist frequency is then defined as follows:

  • fNy=12p, where p is the largest value that allows ti=t0+nip to be possible for all ti values, and ni comprises integers (Eyer and Bartholdi1999). This value corresponds to the same Nyquist frequency as that of our non-gapped data, which is 0.1 min−1.

  • This approach is more appropriate for GLS and HSF than applying a “pseudo-Nyquist” limit based on an average or a minimum value of Δt (Scargle1982; Vanderplas2018).

To estimate β, the spectra are fitted by taking the following steps:

  • A maximum likelihood estimator (MLE) is employed to determine the fit parameter β (Duvall and Harvey1986).

  • The MLE fit involves minimizing the negative log-likelihood function −ln L(O) of observations Oi at frequency fi using the equation

    (9) - ln L ( O ) = i = 1 n ln O i + O i O i ,

    where Oi refers to the power-law model c(1fβ) being fitted, with c as a normalization coefficient.

  • The MLE fit is recommended over least-squares regression because the latter assumes a Gaussian distribution of periodogram residuals, leading to a biased estimate of β (Clauset et al.2009).

5 Results

5.1 Single-frequency signal

First, we show an example of the time series generated by the simulation described in Sect. 3.1, which consists of a 0.5 h wave in a 6 h time series. As can be seen in Fig. 1a, in the absence of gaps, an accurate estimation of the 4 K amplitude and the 0.5 h period of this signal is acquired by both GLS and FFT. In addition, the spectra obtained by FFT and GLS are (as expected) equivalent in the case of evenly sampled data (Scargle1982). When the random gaps replace 50 % of the data points, a significant difference between the amplitude spectra is observed; see Fig. 1b. Both methods still provide an accurate estimate of the signal's period. Linear interpolation of the 50 % gapped signal preserves the structure of the wave but loses some of the high-frequency components, which leads to a significant underestimation of the amplitude by 43 % in the FFT spectrum. In contrast, the amplitude of the highest peak in the GLS spectrum is not affected by the gaps and has not changed from the expected value.

Figure 1Time series of a 0.5 h wave in a 6 h observation time generated according to Sect. 3.1. Panel (a) shows the time series (upper left) and its temporal amplitude spectrum (lower left) in the absence of gaps. Panel (b) shows the time series (upper right) and its temporal amplitude spectrum (lower right) after the addition of a number of random gaps equal to 50 % of the total data length.


When comparing average relative period bias for different simulation periods, Fig. 2a shows that GLS demonstrates no period bias below 80 % GP and a negligible bias beyond (within ±20 % deviation interval). Figure 2b shows, contrarily, that the smaller the simulated period of the signal, the more FFT overestimates it at GPs larger than 40 %. This is due to linear interpolation replacing the removed high-frequency components with lower-frequency ones, which eventually dominate as the GP gradually increases. To put these results in perspective, at 70 % GP, FFT inaccurately determines the period of a 0.5 h wave as 2.93 h (i.e. overestimates it by 486 %). Given that there are quite few data points left in the time series at 70 % GP, it is expected that FFT is not able to recover the correct period; however, GLS is still capable of obtaining the exact value of the period (i.e. an error of 0 %). Not only is GLS a much better estimator of the period on average, but also since its standard deviation remains trivially small until 80 % GP, it is a more reliable choice on a case-by-case basis.

Figure 2Comparison of relative period bias (Eq. 6) as a function of the gap percentage for GLS (a) and FFT (b), here shown for each simulated period (frequency). The bias is estimated from the average estimated values of the periods, and its standard deviation is scaled accordingly. Note that the y axis (relative period bias) is limited between [-0.2,0.2] for GLS and [-4,20] for FFT. This shows how extremely different the accuracy of each method is.


Similarly, FFT's amplitude bias experiences clear dependency on the frequency of the signal, while GLS demonstrates a negligible amplitude bias at GPs below 80 %; see Fig. 3. In particular, the results show that the mean estimated amplitude from the FFT spectrum deviates exponentially from the expected value as the frequency of the signal increases, even at gap percentages below 50 %. In addition, the standard deviation of amplitude estimation by GLS increases as the GP increases, although it remains within a ±10 % deviation interval up to 80 % GP. On the contrary, the standard deviation of amplitude estimation by FFT significantly increases as both the frequency of the signal and the GP increase, implying that FFT is more inconsistent and highly sensitive to missing data, especially for high-frequency signals.

Figure 3Comparison of relative amplitude bias (Eq. 6) as a function of the gap percentage for GLS (a) and FFT (b), here shown for each simulated period (frequency). According to Fig. 2b, past a certain GP threshold, the highest peak in the FFT spectrum does not belong to the true frequency but to meaningless interpolation noise. For this reason, the amplitude bias reported in Fig. 3 is limited by relative period bias<3, since amplitude values beyond this threshold should not be relied upon. Note that the upper and lower limits of the y axis are different for each method as well.


Overall, GLS provides a more robust estimation of the period and amplitude of gapped time series, while FFT's performance is simultaneously dependent on the GP and frequency of the signal.

5.2 Spectral power-law signal

A time-series example is shown in Fig. 4 to illustrate the complexity of a signal produced according to Sect. 3.2 for β=2, showcasing the signal before and after the introduction of gaps. As the percentage of gaps in the data increases, the impact on the spectral components varies depending on their frequency range. High-frequency components (rapid fluctuations over short periods) are most vulnerable when data are being removed. Subsequently, the lower-frequency components (slow variations over longer periods) follow, exhibiting greater resilience to data gaps. In essence, a considerably greater number of gaps are required to significantly affect the estimation of the latter. A distorted spectrum of this time series can be seen, which is a result of the limited sample length and resolution (Roberts et al.1987). In addition, as the signal comprises numerous waves with random frequencies, the presence of closely spaced frequencies leads to the emergence of complex and broad peaks in the spectrum (Horne and Baliunas1986; Dewan and Grossbard2000).

Figure 4A time-series example (upper left and right) generated by the spectral power-law simulation (according to Sect. 3.2) with a spectral exponent of β=2 within a 6 h observation time, before (a) and after (b) the addition of 50 % gaps. The estimated power spectral densities of both non-gapped and gapped time series are shown in the lower left and right. Here, both x and y axes in the spectra figures were log-scaled so that a linear function can be identified. The dotted lines (in red, green, and blue) represent the fits of the PSD of each method.


In the absence of gaps, the true value of β is accurately estimated from the spectra of all methods (see Fig. 4a). However, after removing 50 % of the data points (Fig. 4b), the estimated spectrum by HSF remains relatively unaltered, while the GLS and FFT spectra diverge. The overestimation of β (bias=0.6) by FFT is due to the amplitudes of high-frequency components being underestimated, which is a result of the interpolation (Schulz and Stattegger1997; Hall and Aso1999). This overestimation of β can also be explained by the fact that these interpolated (high-frequency) components contribute locally by β=3 to the overall slope of the spectrum, which results in positive bias when the true β is less than 3 (Lovejoy2014). In the case of HSF, where no interpolation of the signal is done, an even smaller bias (0.1) results. In contrast, the true power law can be seen for the first few low frequencies in the GLS spectrum; then it starts to flatten at intermediate and high frequencies with a substantial bias of −1.41. This occurs because the lack of data points constrains the least-squares fit by GLS, which leads to the interpolation of power at these frequencies (i.e. a flat line).

For a statistically significant picture, we show the distributions of estimated β values from 1000 simulation runs in Fig. 5. In the non-gapped case (see Fig. 5a), the distributions of all methods overlap within a small standard-deviation range of ±0.2 around β=2. In contrast, we see that gaps cause the estimated β values from the FFT spectra to shift their distribution to higher values around β=2.3 in Fig. 5b, while the distribution of GLS becomes narrower and is diverted far below the expected value of 2. The mean of estimated β values from HSF spectra is clearly the closest to both the true value and its mean in the absence of gaps, which shows consistency and less sensitivity to the gaps. It is worth noting that the distributions in Fig. 5a show that even in the absence of gaps, estimated β lies mostly within the range of [1.5,2.5] and not exactly at 2, as single spectra are distorted without averaging. It is also important to mention that the results of the bias are almost identical whether β is estimated from averaging power-law exponents of single spectra or it is estimated from fitting an averaged spectrum.

Figure 5Histograms of the estimated β values from 1000 spectra of time series generated by the spectral power-law simulation (according to Sect. 3.2) before (a) and after (b) the addition of 50 % gaps. The simulated spectral exponent was β=2, and each time series had a 6 h observation span. The vertical lines in the histograms of the gapped case refer to the mean estimated values of β by each method.


To further explore the behaviour of the bias under different conditions, we evaluated the effect of changing the simulated value of β on the estimation bias (see Fig. 6). At 0 % GP, all methods show no bias, except for β>2 and β=-1, where there is an apparent deviation. The first deviation is expected because the more β is larger than 2, the more the spectrum suffers from “low-frequency leakage” due to the finite length of the time series (Klis1994; Schulz and Mudelsee2002). The other case of small deviation takes place on the opposite end of the spectral slope range (β=-1), where high frequencies dominate. The power at these frequencies is quite easily aliased, hence underestimated, as they are the closest to the Nyquist limit. These deviations do not mean that spectra of β=-1 or β>2 cannot be obtained but that on average, they are very likely to be misestimated.

Figure 6Comparison of the bias in the mean β estimates obtained from FFT (a), GLS (b), and HSF (c) as a function of the gap percentage. The results of each method are shown for spectra with power-law exponents β{-1,0,1,5/3,2,2.5,3}.


In the instance of gapped time series, our results show that as the GP increases, the biases in the estimated exponent also became more pronounced for all three methods. Similarly to the non-gapped situation, GLS demonstrates an exceptional efficiency in estimating flat spectra where β=0 with a negligible bias. This is a consequence of the absence of frequency dependency in a flat spectrum, which renders the gaps irrelevant in terms of introducing bias, since the GLS spectrum is already flat. As β increases (indicating a steeper decline in power with increasing frequency) and the percentage of gaps in the data increases, the bias in the GLS spectrum becomes more prominent. For instance, in the case of β transitioning from 1→3 where power is skewed towards low frequencies, gaps cause GLS to mistakenly assign excessive power at the missing high frequencies, ultimately resulting in a steady underestimation of β. In contrast, when β=-1, high-frequency components dominate low-frequency ones. Consequently, there is a loss of power at these high frequencies as the gaps disrupt their sampling. This loss of power causes the GLS method to overestimate β, mirroring the GLS bias observed when β=1.

In similar fashion, both FFT and HSF demonstrate a relatively constant bias for β=3 of approximately −0.3 at all GPs. However, as β decreases from 2-1, their biases monotonically increase as the GP increases. Nevertheless, HSF shows substantially less bias than FFT when the GP exceeds 10 %. The FFT bias is attributed to the established interpolation effects. Therefore, as more data points are interpolated, the FFT spectrum progressively underestimates the amplitude of the high frequencies. This underestimation results in the bias being positive for all β values, except β=3, where leakage causes FFT to overestimate these frequencies. Overall, averaging β values from single spectra is a good measure of the expected value because of their low standard deviation except at very high GPs.

In light of the aforementioned considerations, it can be argued that the FFT technique demonstrates competence in generating accurate spectral estimations for non-gapped time series. Nevertheless, it encounters challenges as the data incorporate an increasing number of gaps, necessitating interpolation techniques which introduce inherent biases. Meanwhile, HSF is demonstrated to be a particularly reliable approach for analysing GW time series with spectral power-law exponents β{1,5/3,2,2.5,3} and in between. Its performance, however, exhibits limitations primarily in cases where the spectrum of a time series possesses a power-law exponent β<1. Notably, such occurrences have only been observed and predicted within measurements of vertical wind time series, as indicated in Table A1. One can reasonably anticipate that spectra falling within the range of β values between 1 and 3 will be prevalent across the majority of atmospheric time series.

Conversely, the GLS method yields similarly favourable outcomes, particularly for time series whose spectra are flat and high-frequency-dominated, so that it even surpasses the accuracy of HSF when β possesses values between −1 and 0. Nonetheless, the GLS method exhibits an increasing bias as the value of β increases beyond 0, rendering it a less certain choice for timescales extending beyond a few hours, which commonly occur in the context of atmospheric gravity waves. Clearly, the consistent and overall impeccable GLS performance in estimating periods and amplitudes of signals with single frequencies does not seem to translate to universally resolving the level of superposition of many random periodicities with power-law amplitudes.

6 Discussion

6.1 Low-frequency leakage

The problem of power leakage from low frequencies into higher ones arises as a result of the constrained frequency range, which itself is limited by the observed time span T. This leakage takes place not only in the case of spectra with β>2, but also in the spectra of time series with periods longer than T. GWs can often have these kinds of periods longer than the simulated time span of 6 h, and normally these periods are not resolved. Thus, we also quantitatively tested the effect of these long periods on the estimated spectra by adding three extra waves of 8, 10, and 12 h periods into the simulated time series. Two cases were examined, one where each of the three waves has an amplitude equivalent to the lowest-frequency component in each simulation. In another case, we scaled their amplitudes by the same power-law exponent β as all frequencies in the simulation.

In both cases, longer-than-T periods produced quite similar effects on the spectra, with a substantial positive bias observed for all signals with β<2 and a significant negative bias for β>2, even in the absence of gaps. For instance, Fig. 7a and b show an example of this leakage in averaged spectra for β=3, without and with the extra waves. A common feature of both cases is the spectral power being excessively concentrated at the lowest frequencies. When the extra long waves are added, the leakage becomes more drastic for GLS and FFT, and their absolute biases of β increase. In contrast, HSF is less affected by these long periods compared to FFT and GLS. This effect of longer-than-T waves in the time series resembles that of trends, which contribute to the spectral shape with a power-law exponent β=2 (Klis1994).

Figure 7Averaged temporal spectra of non-gapped time series generated by the spectral power-law simulation with a spectral exponent of β=3 within a 6 h observation time: (a) non-prewhitened, without extra long waves; (b) non-prewhitened, after the addition of three extra waves with frequencies lower than f=1/T (particularly 8, 10, and 12 h) to the simulation; and (c) the postdarkened spectra of the prewhitened time series with extra waves. Here, both axes were also log-scaled so that a linear function can be identified. The dotted lines (in red, green, and blue) represent the fits of the PSD of each method.


While a weighted fit of the spectra can reduce the bias, it does not fully rectify the problem of leakage; it also requires a smoothed spectrum and may be confounded by other biases from observational noise, gaps, or method inefficiencies. One approach to counteract this leakage from unresolved low-frequency power is based on “prewhitening” the time series and then “postdarkening” the spectra (Blackman and Tukey1958). Prewhitening is a technique to decorrelate the time series (reduce its autocorrelation close to that of a white noise signal) before calculating the PSD (Keisler and Rhyne1976; Dewan and Grossbard2000; Guharay and Sekar2011). This step serves to decrease the rate of change in PSD(f) with frequency f; hence, the distribution of the power is more even. The prewhitened time series is obtained by the means of first differencing, where each data point is subtracted from its successive one; i.e. Δx(t)=x(t)-x(t-Δt). By doing this, the contribution of the mean becomes negligible and trends are transformed into constants (Bieber et al.1993). The spectrum of this prewhitened series Δx(t) is related to the original spectrum of x(t) by

(10) PSD x ( t ) ( f ) = PSD Δ x ( t ) ( f ) 2 ( 1 - cos ( 2 π f Δ t ) ) .

This factor 2(1−cos (2πfnΔt)) is derived from the Fourier transform of the autocorrelation function of Δx(t) and is used to compensate for the prewhitening process (Houbolt et al.1964). This compensation postdarkens (recolours) the resulting spectrum. Note that the spectrum of the prewhitened series Δx(t) should be smoothed first using a Hann window to suppress the random fluctuations before postdarkening.

In Fig. 7c, we present the postdarkened spectra after prewhitening the time series for β=3. This approach completely cancels out the bias in all methods for both cases. This confirms the effectiveness of the prewhitening and postdarkening method in correcting the leakage problem. However, this approach is not a perfect solution, since it may introduce additional bias for less steep spectra (where β<1) which do not suffer from leakage. When the spectrum is constant in the low-frequency part, the prewhitening process might cause it to be distorted as well.

6.2 Method selection procedure

The spectral analysis of GW time-series data is a complex task that requires careful consideration of various factors. Based on our simulation results, we propose a flowchart (see Fig. 8) that outlines a practical guide for selecting appropriate spectral estimation methods for GW studies, taking into account the characteristics of the observed data such as their complexity and percentage of gaps in them. From the flowchart, it is clear that there are recognizable differences between the patterns of the time series of superposed waves and signals with one frequency. Even within the former classification, the anti-persistent time series whose spectra have β[0,-1] are still differentiable from those long-range-dependent ones where β[1,3]. It is also safe to say that theoretical predictions of GW spectra with β[0,-1] exist only for measurements of vertical wind time series (see Table A1). Otherwise, β[1,3] spectra should be expected for the vast majority of atmospheric time series. Note that, even in the absence of gaps in the signal, caution must be taken if the estimated β is approximately equal to or less than 2. If it is, then it can be an accurate estimation or caused by one of the systematic errors such as interpolation or longer-than-T variations.

Figure 8Recommended procedure for estimating power spectra of gravity wave time series.


7 Conclusions

The domain of observational analyses, especially in atmospheric physics, is vast and complex. The intrinsic nature of these measurements requires precision, reliability, and adaptability of analytical methods to extract meaningful insights. Observational gaps due to instrument failures or adverse conditions are also common in atmospheric physics. Our research on synthetic atmospheric gravity wave (GW) time series offers a comprehensive overview of different spectral estimation methods, shedding light on their strengths and limitations. The methods compared in this study are the fast Fourier transform (FFT) (which requires interpolation of gaps) and the generalized Lomb–Scargle method (GLS) and the Haar structure function (HSF) (both which can handle unevenly sampled data without interpolation). Their performance is assessed by evaluating whether the output spectra of simulated time series (with known a priori features and gap distributions) match the input parameters of the data, including the frequency, amplitude, and spectral slope. Building on our findings, we propose the following recommendations:

Determination of the time-series nature. Initially, the series is assessed to determine if it represents a pattern of a signal with one frequency. If it does, the analysis proceeds directly to using the GLS method, especially if the time series contains gaps and/or demonstrates rapid variations. Rather than interpolating these gaps to use FFT, which alters the data's characteristics, GLS offers far more accurate spectral estimates.

Evaluation of data gaps. For more complex series, particularly those characterized by a superposition of waves or the power-law spectrum, the data are first assessed for observational gaps. The presence of such gaps profoundly affects the accuracy of traditional spectral-analysis methods like FFT. If variations are rapid (i.e. the slope of that power-law spectrum β is in the range [-1,0]), the series is again subjected to the GLS method. However, for slower variations such as the vast majority of atmospheric measurements, represented by a β value in the range of [1,3], the HSF method is consistently the least biased.

Mitigate low-frequency leakage. When the spectra of the time series with slow variations are analysed (even in the absence of gaps), it should be investigated whether low-frequency leakage takes place. Finite time series, very steep spectra, and waves with periods longer than the observed time span cause this kind of leakage. All tested methods, to varying extents, struggled with this leakage from unresolved low-frequency power, leading to biases in estimating the spectral slope β. Thus we recommend cautiously applying the approach of prewhitening the time series and postdarkening the spectrum to address this problem. Prewhitening can be actualized by first-differencing the series, which renders its spectrum close to that of white noise. Subsequently, postdarkening the spectrum is required to compensate for the prewhitening process via division by a frequency-dependent factor. Nevertheless, depending on the characteristics of the analysed time series, further steps (e.g. windowing) might be necessary to help counter other causes of leakage.

These sequential decisions highlight the importance of tailoring analysis methods to the characteristics of the data at hand. Employing a one-size-fits-all method can result in biased spectra, which are especially critical in GW parameterizations in all atmospheric models. Our recommendations, grounded in rigorous research on synthetic gravity wave time series, aim to guide researchers in making informed decisions, ensuring the accuracy of spectral results and advancing our understanding of the dynamics of the atmosphere.

Appendix A
(VanZandt1982, and ref. within)(Dewan and Good1986; Smith et al.1987)(Dewan1994)(Dewan1994)(Dewan1994)(Dewan1994)(Shibata et al.1988)(Gardner1994)(Gardner1994)(Gardner1994)(Gardner1994)(Weinstock1990; Zhu1994)(Weinstock1990; Zhu1994)(Hines1991)(Hines1991)(Gardner et al.1995)(Gardner et al.1995)(Gardner et al.1995)(Gardner et al.1995)(Gardner et al.1995)(Gardner et al.1995)(Zhang et al.2017)(Zhang et al.2017)(Zhang et al.2017)(He et al.2020)(Hertzog and Vial2001)(Hertzog and Vial2001)(Podglajen et al.2016)

Table A1Comparison of theoretical predictions and selected observed values of the power-law exponent β of GW spectra. Here T refers to temperature, W is wind, and ρ is density.

Download Print Version | Download XLSX

Code availability

All codes used in this paper are available under Mossad (2023).

Author contributions

MM wrote the codes to conduct these simulations, analysed their results, and drafted the manuscript. IS, RW, and GB provided supervision and scientific insight and edited the text of the manuscript.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Measurement Techniques. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


This paper is a contribution to the project W1 (Gravity Wave Parameterization for the Atmosphere) of the Collaborative Research Centre TRR 181 “Energy Transfers in Atmosphere and Ocean” and the project Analyzing the Motion of the Middle Atmosphere Using Nighttime RMR-lidar Observations at the Midlatitude Station Kühlungsborn (AMUN).

Financial support

This research has been supported by the Deutsche Forschungsgemeinschaft (grant nos. 274762653 and 445400792).

The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Lars Hoffmann and reviewed by two anonymous referees.


Alexander, M., Tsuda, T., and Vincent, R. N.: Latitudinal Variations Observed in Gravity Waves with Short Vertical Wavelengths, J. Atmos. Sci., 59, 1394–1404,<1394:LVOIGW>2.0.CO;2, 2002. a

Axford, D.: Spectral analysis of an aircraft observation of gravity waves, Q. J. Roy. Meteor. Soc., 97, 313–321,, 1971. a

Babu, P. and Stoica, P.: Spectral analysis of nonuniformly sampled data – a review, Digit. Signal Process., 20, 359–378,, 2010. a

Beldon, C. L. and Mitchell, N. J.: Gravity wave–tidal interactions in the mesosphere and lower thermosphere over Rothera, Antarctica (68 S, 68 W), J. Geophys. Res., 115, D18101,, 2010. a

Beres, J., Garcia, R. R., Boville, B. A., and Sassi, F.: Implementation of a gravity wave source spectrum parameterization dependent on the properties of convection in the Whole Atmosphere Community Climate Model (WACCM), J. Geophys. Res., 110, D10108,, 2005. a

Bieber, J. W., Chen, J., Matthaeus, W. H., Smith, C. W., and Pomerantz, M. A.: Long-term variations of interplanetary magnetic field spectra with implications for cosmic ray modulation, J. Geophys. Res., 98, 3585–3603,, 1993. a

Billah, K. Y. R. and Shinozuka, M.: Numerical method for colored-noise generation and its application to a bistable system, Phys. Rev. A, 42, 7492–7495,, 1990. a

Blackman, R. B. and Tukey, J. W.: The measurement of power spectra from the point of view of communications engineering – Part I, Bell Syst. Tech. J., 37, 185–282,, 1958. a

Brown, T. M. and Christensen-Dalsgaard, J.: A technique for estimating complicated power spectra from time series with gaps, Astrophys. J., 349, 667,, 1990.  a

Clauset, A., Shalizi, C. R., and Newman, M.: Power-law distributions in empirical data, Siam Rev., 51, 661–703,, 2009. a

Cooley, J. and Tukey, J. W.: An algorithm for the machine calculation of complex Fourier series, Math. Comput., 19, 297–301,, 1965. a, b

Crowley, G. and Williams, P.: Observations of the source and propagation of atmospheric gravity waves, Nature, 328, 231–233,, 1987. a

Dewan, E.: The saturated-cascade model for atmospheric gravity wave spectra, and the wavelength-period (W-P) relations, Geophys. Res. Lett., 21, 817–820,, 1994. a, b, c, d, e, f, g

Dewan, E. M. and Good, R. E.: Saturation and the “universal” spectrum for vertical profiles of horizontal scalar winds in the atmosphere, J. Geophys. Res., 91, 2742,, 1986. a, b, c

Dewan, E. M. and Grossbard, N.: Power spectral artifacts in published balloon data and implications regarding saturated gravity wave theories, J. Geophys. Res., 105, 4667–4683,, 2000. a, b, c

Duvall, T. L. and Harvey, J. T.: Solar Doppler shifts: sources of continuous spectra, Springer,, 1986. a

Eckermann, S. D. and Hocking, W. K.: Effect of superposition on measurements of atmospheric gravity waves: A cautionary note and some reinterpretations, J. Geophys. Res., 94, 6333,, 1989. a, b

Ern, M., Trinh, Q. T., Preusse, P., Gille, J. C., Mlynczak, M. G., Russell III, J. M., and Riese, M.: GRACILE: a comprehensive climatology of atmospheric gravity wave parameters based on satellite limb soundings, Earth Syst. Sci. Data, 10, 857–892,, 2018. a

Eyer, L. and Bartholdi, P.: Variable stars: Which Nyquist frequency?, Astron. Astrophys. Sup., 135, 1–3,, 1999. a

Ferraz-Mello, S.: Estimation of Periods from Unequally Spaced Observations, Astron. J., 86, 619,, 1981. a

Fritts, D. C.: A review of gravity wave saturation processes, effects, and variability in the middle atmosphere, Pure Appl. Geophys., 130, 343–371,, 1989. a

Fritts, D. C. and Alexander, M.: Gravity wave dynamics and effects in the middle atmosphere, Rev. Geophys., 41, 1003,, 2003. a

Fritts, D. C. and VanZandt, T. E.: Spectral Estimates of gravity wave energy and momentum fluxes. Part I: Energy dissipation, acceleration, and constraints, J. Atmos. Sci., 50, 3685–3694,<3685:SEOGWE>2.0.CO;2, 1993. a

Fritts, D. C., Tsuda, T., Sato, T., Fukao, S., and Kato, S.: Observational evidence of a saturated gravity wave spectrum in the troposphere and lower stratosphere, J. Atmos. Sci., 45, 1741–1759,<1741:OEOASG>2.0.CO;2, 1988. a

Gardner, C. S.: Diffusive filtering theory of gravity wave spectra in the atmosphere, J. Geophys. Res., 99, 20601,, 1994. a, b, c, d, e, f

Gardner, C. S., Tao, X., and Papen, G. C.: Simultaneous lidar observations of vertical wind, temperature, and density profiles in the upper mesophere: Evidence for nonseparability of atmospheric perturbation spectra, Geophys. Res. Lett., 22, 2877–2880,, 1995. a, b, c, d, e, f, g, h

Gerding, M., Höffner, J., Lautenbach, J., Rauthe, M., and Lübken, F.-J.: Seasonal variation of nocturnal temperatures between 1 and 105 km altitude at 54 N observed by lidar, Atmos. Chem. Phys., 8, 7465–7482,, 2008. a

Guharay, A. and Sekar, R.: Seasonal characteristics of gravity waves in the middle atmosphere over Gadanki using Rayleigh lidar observations, J. Atmos. Sol.-Terr. Phy., 73, 1762–1770,, 2011. a, b, c

Haar, A.: On the theory of orthogonal function systems, Math. Ann., 69, 331–371, 1910. a

Hall, C. M. and Aso, T.: Mesospheric velocities and buoyancy subrange spectral slopes determined over Svalbard by ESR, Geophys. Res. Lett., 26, 1685–1688,, 1999. a, b, c, d

Hamilton, K.: The role of parameterized drag in a troposphere-stratosphere-mesosphere general circulation model, Springer,, 1997. a

He, Y., Sheng, Z., and He, M.: Spectral Analysis of Gravity Waves from Near Space High-Resolution Balloon Data in Northwest China, Atmosphere, 11, 133,, 2020. a

Hébert, R.: RScaling (1.0.0), Zenodo,, 2021. a

Hébert, R., Rehfeld, K., and Laepple, T.: Comparing estimation techniques for temporal scaling in palaeoclimate time series, Nonlin. Processes Geophys., 28, 311–328,, 2021. a

Hertzog, A. and Vial, F.: A study of the dynamics of the equatorial lower stratosphere by use of ultra-long-duration balloons: 2. Gravity waves, J. Geophys. Res., 106, 22745–22761,, 2001. a, b

Hines, C. O.: Internal atmospheric gravity waves at ionospheric heights, Can. J. Phys., 38, 1441–1481,, 1960. a

Hines, C. O.: The Saturation of Gravity Waves in the Middle Atmosphere. Part II: Development Of Doppler-Spread Theory, J. Atmos. Sci., 48, 1361–1379,<1361:TSOGWI>2.0.CO;2, 1991. a, b, c, d

Holton, J. M.: The influence of gravity wave breaking on the general circulation of the middle atmosphere, J. Atmos. Sci., 40, 2497–2507,<2497:TIOGWB>2.0.CO;2, 1983. a

Horne, J. H. and Baliunas, S. L.: A prescription for period analysis of unevenly sampled time series, Astrophys. J., 302, 757,, 1986. a, b

Houbolt, J. C., Steiner, R., and Pratt, K. G.: Dynamic response of airplanes to atmospheric turbulence including flight data on input and response, vol. 199, National Aeronautics and Space Administration, (last access: 25 January 2024), 1964. a

Houchi, K., Stoffelen, A., Marseille, G.-J., and De Kloe, J.: Comparison of wind and wind shear climatologies derived from high-resolution radiosondes and the ECMWF model, J. Geophys. Res., 115, D22123,, 2010. a

Keisler, S. R. and Rhyne, R. H.: An assessment of prewhitening in estimating power spectra of atmospheric turbulence at long wavelengths, Tech. rep., NASA Langley Research Center, (last access: 25 January 2024), 1976. a, b

Kirchner, J. W.: Aliasing in 1/fα noise spectra: Origins, consequences, and remedies, Phys. Rev. E, 71, 066110,, 2005. a

Klis, M.: Rapid variability in x-ray binaries–towards a unified description, NATO Science Series C, 1995 edn., Springer, Dordrecht, Netherlands, (last access: 25 January 2024), 1994. a, b, c

Lepot, M., Aubin, J.-B., and Clemens, F.: Interpolation in Time Series: An introductive overview of existing methods, their performance criteria and uncertainty assessment, Water, 9, 796,, 2017. a

Lindgren, E., Sheshadri, A., Podglajen, A., and Carver, R. W.: Seasonal and latitudinal variability of the gravity wave spectrum in the lower stratosphere, J. Geophys. Res.-Atmos., 125, e2020JD03285,, 2020. a

Lindzen, R. S.: Turbulence and stress owing to gravity wave and tidal breakdown, J. Geophys. Res., 86, 9707,, 1981. a

Liu, H., McInerney, J., Santos, S. A., Lauritzen, P. H., Taylor, M. J., and Pedatella, N.: Gravity waves simulated by high-resolution Whole Atmosphere Community Climate Model, Geophys. Res. Lett., 41, 9106–9112,, 2014. a

Lomb, N. R.: Least-squares frequency analysis of unequally spaced data, Astrophys. Space Sci., 39, 447–462,, 1976. a

Lovejoy, S.: A voyage through scales, a missing quadrillion and why the climate is not what you expect, Clim. Dynam., 44, 3187–3210,, 2014. a, b

Lovejoy, S. and Schertzer, D.: Haar wavelets, fluctuations and structure functions: convenient choices for geophysics, Nonlin. Processes Geophys., 19, 513–527,, 2012. a, b, c

Maekawa, Y., Fukao, S., Sato, T., Kato, S., and Woodman, R. F.: Internal inertia–gravity waves in the tropical lower stratosphere observed by the Arecibo Radar, J. Atmos. Sci., 41, 2359–2367,<2359:IIWITT>2.0.CO;2, 1984. a, b

Marinna, A. M., Alfredo, L. A., and Christopher, R. S.: Using the Lomb-Scargle method for wave statistics from gappy time series, IEEE Conference Proceedings, 2019, 1–9,, 2019. a

Meisel, D. D.: Fourier transforms of data sampled at unequal observational intervals, Astron. J., 83, 538,, 1978. a

Mossad, M.: GapWaveSpectra (1.0.0), Zenodo [code],, 2023. a, b

Munteanu, C., Negrea, C., Echim, M., and Mursula, K.: Effect of data gaps: comparison of different spectral analysis methods, Ann. Geophys., 34, 437–449,, 2016. a, b

Muraoka, Y., Sugiyama, T., Kawahira, K., Sato, T., Tsuda, T., Fukao, S., and Kato, S.: Cause of a monochromatic inertia-gravity wave breaking observed by the MU radar, Geophys. Res. Lett., 15, 1349–1352,, 1988. a

Narendra Babu, A., Kishore Kumar, K., Kishore Kumar, G., Venkat Ratnam, M., Vijaya Bhaskara Rao, S., and Narayana Rao, D.: Long-term MST radar observations of vertical wave number spectra of gravity waves in the tropical troposphere over Gadanki (13.5 N, 79.2 E): comparison with model spectra, Ann. Geophys., 26, 1671–1680,, 2008. a

Nastrom, G. D., Van Zandt, T. E., and Warnock, J. M.: Vertical wavenumber spectra of wind and temperature from high-resolution balloon soundings over Illinois, J. Geophys. Res., 102, 6685–6701,, 1997. a

Pichon, A. L., Assink, J., Heinrich, P., Blanc, E., Charlton-Perez, A., Lee, C. S., Keckhut, P., Hauchecorne, A., Rüfenacht, R., Kämpfer, N., Drob, D. P., Smets, P., Evers, L., Ceranna, L., Pilger, C., Ross, O. A., and Claud, C.: Comparison of co-located independent ground-based middle atmospheric wind and temperature measurements with numerical weather prediction models, J. Geophys. Res.-Atmos., 120, 8318–8331,, 2015. a

Pinel, J. and Lovejoy, S.: Atmospheric waves as scaling, turbulent phenomena, Atmos. Chem. Phys., 14, 3195–3210,, 2014. a

Podglajen, A., Hertzog, A., Plougonven, R., and Legras, B.: Lagrangian temperature and vertical velocity fluctuations due to gravity waves in the lower stratosphere, Geophys. Res. Lett., 43, 3543–3553,, 2016. a

Qing, H., Zhou, C., Zhao, Z., Chen, G., Ni, B., Gu, X., Yang, G., and Zhang, Y.: A statistical study of inertia gravity waves in the troposphere based on the measurements of Wuhan Atmosphere Radio Exploration (WARE) radar, J. Geophys. Res.-Atmos., 119, 3701–3714,, 2014. a, b

Rao, N. K., Ratnam, M. V., Vedavathi, C., Tsuda, T., Murthy, B. V. K., Sathishkumar, S., Gurubaran, S., Kumar, K. S., Subrahmanyam, K. V., and Rao, S. V.: Seasonal, inter-annual and solar cycle variability of the quasi two day wave in the low-latitude mesosphere and lower thermosphere, J. Atmos. Sol.-Terr. Phy., 152–153, 20–29,, 2017. a

Rice, S. O.: Mathematical analysis of random noise, Bell Syst. Tech. J., 23, 282–332,, 1944. a

Rigling, B. D.: Application of temporal gap filling to natural power law spectrums, IEEE Geosci. Remote S., 9, 624–628,, 2012. a

Roberts, D. A., Lehar, J., and Dreher, J.: Time Series Analysis with Clean – Part One – Derivation of a Spectrum, Astron. J., 93, 968,, 1987. a

Scargle, J. D.: Studies in astronomical time series analysis. II – Statistical aspects of spectral analysis of unevenly spaced data, Astrophys. J., 263, 835,, 1982. a, b, c

Schulz, M. and Mudelsee, M.: REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series, Comput. Geosci., 28, 421–426,, 2002. a

Schulz, M. and Stattegger, K.: Spectrum: spectral analysis of unevenly spaced paleoclimatic time series, Comput. Geosci., 23, 929–945,, 1997. a, b

Schuster, A.: On the investigation of hidden periodicities with application to a supposed 26 d period of meteorological phenomena, J. Geophys. Res., 3, 13,, 1898. a

Shibata, T., Ichimori, S., Narikiyo, T., and Maeda, M.: Spectral analysis of vertical temperature profiles observed by a lidar in the upper stratosphere and the lower mesosphere, J. Meteorol. Soc. Jpn. Ser. II, 66, 1001–1005,, 1988. a

Shinozuka, M.: Simulation of Multivariate and Multidimensional Random Processes, J. Acoust. Soc. Am., 49, 357–368,, 2005. a

Sica, R. J. and Russell, A. G.: How many waves are in the gravity wave spectrum?, Geophys. Res. Lett., 26, 3617–3620,, 1999. a, b

Smith, A. M.: Global Dynamics of the MLT, Surv. Geophys., 33, 1177–1230,, 2012. a

Smith, S. M., Fritts, D. C., and VanZandt, T. E.: Evidence for a saturated spectrum of atmospheric gravity waves, J. Atmos. Sci., 44, 1404–1410,<1404:EFASSO>2.0.CO;2, 1987. a, b, c

Song, I.-S. and Chun, H.-Y.: A Lagrangian Spectral Parameterization of Gravity Wave Drag Induced by Cumulus Convection, J. Atmos. Sci., 65, 1204–1224,, 2008. a

Swenson, G. R., Haque, R., Yang, W., and Gardner, C. S.: Momentum and energy fluxes of monochromatic gravity waves observed by an OH imager at Starfire Optical Range, New Mexico, J. Geophys. Res., 104, 6067–6080,, 1999. a

Vanderplas, J.: Understanding the lomb–scargle periodogram, Astrophys. J. Suppl. S., 236, 16,, 2018. a, b

VanZandt, T.: A universal spectrum of buoyancy waves in the atmosphere, Geophys. Res. Lett., 9, 575–578,, 1982. a, b

Vio, R., Andreani, P., and Biggs, A. D.: Unevenly-sampled signals: a general formalism for the Lomb-Scargle periodogram, Astron. Astrophys., 519, A85,, 2010. a

Weinstock, J.: Nonlinear Theory of Gravity Waves: Momentum Deposition, Generalized Rayleigh Friction, and Diffusion, J. Atmos. Sci., 39, 1698–1710,<1698:NTOGWM>2.0.CO;2, 1982. a

Weinstock, J.: Saturated and unsaturated spectra of gravity waves and scale-dependent diffusion, J. Atmos. Sci., 47, 2211–2226,<2211:SAUSOG>2.0.CO;2, 1990. a, b, c, d

Zechmeister, M. and Kürster, M.: The generalised Lomb-Scargle periodogram. A new formalism for the floating-mean and Keplerian periodograms, Astron. Astrophys., 496, 577–584, 2009. a, b

Zhan, Q., Manson, A. H., and Meek, C.: The impact of gaps and spectral methods on the spectral slope of the middle atmospheric wind, J. Atmos. Terr. Phys., 58, 1329–1336,, 1996. a

Zhang, S., Huang, C., Huang, K., Gong, Y., Chen, G., Gan, Q., and Zhang, Y.: Latitudinal and seasonal variations of vertical wave number spectra of three-dimensional winds revealed by radiosonde observations, J. Geophys. Res.-Atmos., 122, 13174–13190,, 2017. a, b, c, d

Zhang, S. D., Huang, C., and Yi, F.: Radiosonde observations of vertical wave number spectra for gravity waves in the lower atmosphere over Central China, Ann. Geophys., 24, 3257–3265,, 2006.  a, b, c

Zhang, S.-N., Peterson, R. N., Wiens, R. H., and Shepherd, G. G.: Gravity waves from O2 nightglow during the AIDA '89 campaign I: emission rate/temperature observations, J. Atmos. Terr. Phys., 55, 355–375,, 1993. a

Zhu, X.: A new theory of the saturated gravity wave spectrum for the middle atmosphere, J. Atmos. Sci., 51, 3615–3626,<3615:ANTOTS>2.0.CO;2, 1994. a, b

Short summary
This numerical study addresses observational gaps' impact on atmospheric gravity wave spectra. Three methods, fast Fourier transform (FFT), generalized Lomb–Scargle periodogram (GLS), and Haar structure function (HSF), were tested on synthetic data. HSF is best for spectra with negative slopes. GLS excels for flat and positive slopes and identifying dominant frequencies. Accurately estimating these aspects is crucial for understanding gravity wave dynamics and energy transfer in the atmosphere.