Uncertainty of Wave Spectral Shape and Parameters Associated with the Spectral Estimation

Abstract

The uncertainty in estimating the wave spectrum from the records of wave elevation by heave–pitch–roll buoys is studied, considering the effects of the estimation method and the spectral resolution adopted in the process. This investigation utilizes measurements from a wave buoy moored in deep water in the South Atlantic Ocean. First, the spectra are computed using the autocorrelation function and the direct Fourier method. Second, the spectral resolution is tested in terms of degrees of freedom. The degrees of freedom are varied, and the resulting spectra and integrated parameters are computed, showing significant variability. A simple and robust methodology for determining the wave spectrum is suggested, which involves calculating the average energy density in each frequency band. The results of this methodology reduce the variability of the estimated parameters, improving overall accuracy while preserving frequency resolution, which is crucial in complex sea states. Additionally, to demonstrate the feasibility of the implemented approach, the final spectrum is fitted using an empirical model ideal for that type of spectrum. Finally, the performance and the goodness of the fit process for the final averaged curve are checked by widely used statistical metrics, such as R² = 0.97 and root mean square error = 0.49.

1. Introduction

Waves in the ocean are one of the most popular and prominent sources of energy in the marine environment. The studies of this energy and its propagation are important for forecasts and alerts and play an essential role in hydrodynamics due to environmental loads. Independent of the application, their behavior is highly complex and random. Past theoretical advances and recent computational progress in applying statistics to ocean and coastal engineering allow the assessment of the responses of various ocean structures to the loads generated by the waves. In this sense, the wave spectrum is an essential input to different evaluations concerning the design or the operation of marine structures and ship navigation, being indispensable in the calculations of stochastic loading and the response of floating bodies to wave effects [1,2].

The design stage of any marine structure requires the distributions of significant wave heights (H_s) and characteristic periods (mean or peak) to effectively model the safety domain and operational criteria. Significant wave height (H_s) and, e.g., associated peak period (T_p) are essential wave parameters that are used for assessing the response spectrum of the structures, and different approaches are available to model the probabilistic relation between these variables [3,4,5]. Nevertheless, the waves, like many other ocean-meteorological processes, usually present patterns of variation at different time scales (e.g., daily, annual, and seasonal). Consequently, from a practical point of view, knowing the variability and average properties of the wave conditions is of great interest. In particular, knowledge of the climatic average and the variability of the wave spectrum represents an aspect of great practical interest. For example, robust ways to estimate and represent climatic wave spectra have been presented by Rodriguez et al. [6] and Lucas et al. [7].

It is known that the wave parameters obtained from the spectrum are statistical estimates with different levels of associated uncertainty. The confidence intervals can also be calculated, quantifying thus the robustness of the results. Several authors addressed these questions such as Donelan and Pierson [8], who studied the sampling variability of estimates of spectra, Medina et al. [9] who discussed distortions associated with random sea simulators, Young [10] investigated the probability distribution of spectral integrals, while Elgar [11] studied the bias of effective degrees of freedom of a spectrum. The uncertainty of the shape of the tail of the spectrum was also studied by several authors such as Donelan et al. [12], Prevosto et al. [13], as well as Rodriguez and Guedes Soares [14].

Rodriguez et al. [15] quantified the uncertainty of the sea-state parameters resulting from the methods of spectral estimation and demonstrated that these intervals are functions of the spectral density curve integration process and are dependent on the shape of the spectrum. As a consequence, the stability and accuracy of the wave parameters are functions of the number of degrees of freedom used to estimate the spectrum. In short, the expected confidence intervals in the wave parameter computation performance are functions of the effective degrees of freedom.

The accuracy of the wave parameters obtained from the spectral density function is influenced by some external factors that cause uncertainty. These uncertainties must be identified, studied, and reduced whenever possible, as they are from different sources and natures, and they can induce uncertainties in the spectral estimations in various ways. Usually, the uncertainties related to sea wave description can be divided into two groups: (i) aleatory uncertainty (random nature) and (ii) epistemic uncertainty (knowledge-based). Uncertainties type (i) represents the natural randomness of a quantity, also known as intrinsic or inherent uncertainty, e.g., the wave height changeability over time. This particular aleatory uncertainty cannot be reduced or eliminated. Epistemic uncertainty represents modeling errors that can be reduced by collecting more information about the quantity or refining the methods of measurement as well as the statistics used.

Bitner-Gregersen and Hagen [16] have discussed these uncertainties, which they classified as data uncertainty, statistical uncertainty (related to sampling variability), model uncertainty, and climatic uncertainty. Data uncertainty is derived by instrumental imperfection employed to measure the quantity or due to incorrect models used to estimate the data. In other words, if a quantity considered is not obtained directly from a measurement but via a numerical empirical estimation process. Statistical uncertainty, also known as sampling variability, is a consequence of a limited number of observations or when relevant periods in the time series are missing, generating sampling biases. This may refer to short-term time series or long-term histories of sea wave characteristics. Climatic uncertainty or climatic variability, often regarded as model uncertainty, addresses the representativeness of a measured or simulated variable. The climate uncertainty is due to the natural variability of met-ocean climate and anthropogenic climate change. These general definitions of uncertainties are included in the DNV rules on reliability analysis of marine structures and detailed in DNV [17] and by Bitner-Gregersen et al. [18,19] and Bitner-Gregersen and Magnusson [20].

The above-mentioned uncertainties, which result from the environment’s random variations, cannot be reduced; however, the uncertainties due to the use of an incorrect methodology must be validated and reduced whenever possible. Consequently, the wave parameters, as the final product of some methodology to estimate the spectra, contain model uncertainty. Rodriguez et al. [21] who examined different methods of spectral estimation (Blackman–Tukey, direct Fourier transform, and maximum entropy) concluded that there was a significant model uncertainty for all sea-state parameters considered. The sea wave description will be affected to varying degrees by all types of epistemic uncertainties. In this vein, the identification of these uncertainties and their quantification represents important information for the planning of any marine activities.

The most common uncertainties in the wave parameters computed from the spectral density function can be considered the mentioned model uncertainties, which add statistical variability to the approximations. The implementation and performance of the model affect the stability and accuracy of the final parameters. The model uncertainty is induced by (a) the procedure for numerical integration, (b) frequency cut-off selection, (c) the wave record length, and (d) frequency resolution, or similarly, degrees of freedom of the spectral estimates. Realizing that these effects are extremely important in determining the stability of the estimated spectral parameters, several studies in the scope of the characterization and quantification of these factors were carried out.

Several authors have addressed this type of uncertainty. Rye [22] investigated the stability of the main wave parameters. Chakrabarti and Cooley [23] studied the statistical distribution of wave periods and heights, and Arhan [24] analyzed the variability in spectral estimation using long continuous wave records. In addition, Cavanié [25] focused on the standard error in the estimation of mean and significant wave heights of finite length. Mansard and Funke [26] addressed only the statistical variability. Young [27] studied the confidence limits associated with estimates of the peak frequency, and Peña [28] studied the stability of the moments of the maximum entropy wind wave spectrum. Complementary, Rodriguez et al. [29] investigated the uncertainty of the sea-state parameters resulting from the methods of spectral estimation and Gomes and Guedes Soares [30] the estimation of directional spectrum by maximum entropy. Moreover, according to Gao et al. [31,32], irregular waves can influence oscillations and reflections due to the wave interaction with periodic structures or with the bottom topography, inducing resonance in harbors.

The previously mentioned references are highly relevant in ocean and coastal engineering projects, where safety and operational efficiency depend on the control of incident waves. However, some pertinent examples of applications in which the suggested approach could help as a guide to well-practice are cited next. So, correctly modeling the effects of uncertainties for defining the wave climate is vital, e.g., designing maritime structures, impact studies on port operations, assessment of coastal risks and erosion, fatigue analysis and useful life of offshore structures, sediment transport, and coastal dynamics simulations. In general, all the topics mentioned above depend on a correct and robust estimation of wave parameters, and most of the time using spectral evaluation.

The present work deals with the uncertainties that influence the variability of the wave parameters, due to the frequency resolution, or correspondingly, the degrees of freedom of the spectral estimation. The objective is to support the decision of the number of degrees of freedom applied in estimating the spectrum. Knowing that this is an important source of uncertainty, a simple approach is suggested for determining the representative spectrum robustly and accurately. Additionally, the new suggested approach aims to define a final spectrum to be used as a characteristic, resulting from a robust statistical process that reduces uncertainty.

2. Ocean Data

The wave data analyzed are measurements from a heave–pitch–roll buoy (Axys 3M) moored at the position: latitude 31.57° S and longitude 49.87° W. The dataset consists of a heave signal containing 20 min of measurements of the buoy movements (heave, pitch, and roll), per hour, with a sampling interval of 0.78 s (D_f = 1/0.78). The local wave characteristics and the uncertainties of the Axys measurement system have been addressed by Campos et al. [33]. The wave process is assumed stationary over the 20 min intervals of recordings.

The period that registered an extreme event occurred on the 20th and 21st of August 2009. In this sense, the short-term raw data considered to consist of sea surface elevation occurred exactly on 21 August 2009. Figure 1 shows the map of the location as well as the moored buoy managed by the Brazilian Navy. Figure 2 illustrates the storm evolution recorded by the buoy, indicating the maximum significant height on day 21. Figure 3 illustrates the displacement of the sea surface as well as the number of waves recorded by the buoy elevation movement. The method is exemplified by considering the real oceanic data records of these conditions. The methodology used to process the raw data and extract the wave parameters is based on the distribution of the energy over the frequencies, and only a heave signal is used for the power spectrum estimation. The approaches used are described in the next section.

3. Wave Spectrum from Water Surface Movement

The spectrum of an ocean wave can be obtained by the decomposition of the random sea-surface displacement motion into a series of sinusoidal components associated with a specific frequency and elevation. After the process of decomposition of the raw signal, the associated unique frequency and amplitude of each component are computed. This analysis is very powerful since it condenses the energy component present in the raw signal per frequency (period) in the wave decomposition. The final spectrum can be defined as the summation of the amplitude squared over a range of angular frequency. The units are square meters per second (m²s), representing the amount of energy in a given sea state per frequency. This analysis is one of the most important and commonly used in ocean engineering as well as in many other scientific fields.

Nowadays, the literature defines the quantity S(ω) with different names, such as energy spectrum, spectral density, and power spectrum, among others, such as the variance spectrum, since the variance spectrum is defined by autocovariance-based. The most well-known applied approaches to estimate the spectrum from an ocean wave observation are the Autocorrelation Function (ACF) and the Fast Fourier Transform (FFT). This work applies both methods to obtain the wave spectrum [34,35].

3.1. Autocorrelation Function

To compute the spectral density S(ω) from the Autocorrelation Function, R(τ), the Wiener–Khintchine theorem is applied. The Wiener–Khintchine theorem represents R(τ) and S(ω) as a pair of Fourier transforms:

(1)$S\left(\omega \right)={\displaystyle \frac{1}{2\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}R\left(\tau \right)e^{-i\omega \tau }d\tau ={\displaystyle \frac{1}{2\pi }}{\int }_0^{\mathrm{\infty }}R\left(\tau \right)cos\omega \tau d\tau$

(2)$R\left(\tau \right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}S\left(\omega \right)e^{+i\omega \tau }d\omega =2{\int }_0^{\mathrm{\infty }}S\left(\omega \right)cos\omega \tau d\omega$

To specify the ACF numerically by using the sea surface elevation data (ղ), the following relation is used:

(3)$\widehat{R}\left(r∆t\right)={\displaystyle \frac{1}{N_t-r}}{\sum }_{n=1}^{N_t-r}{ղ}_n{ղ}_{n+r}\hspace{1em}\hspace{1em}\hspace{1em}r=0,1,2\dots m$

where r represents the lag number, m is the maximum lag number (m << N_t), Δt is timestep, and N_t is the total number of samples. The magnitude of m is chosen according to the first global minimum of the ACF envelope.

The Hanning Window is used to modify the autocorrelation function and suppress the spectral energy leakage through the side lobes. The windowing process aims to taper the ACF by eliminating the discontinuity at the end of R(τ). The mathematical expression of the Hanning Window is as follows:

(4)$u_h\left(r∆t\right)={\displaystyle \frac{1}{2}}\left(1+cos{\displaystyle \frac{\pi r}{m}}\right)\hspace{1em}\hspace{1em}\hspace{1em}r=0,1,2\dots m$

Eventually, the modified (tapered) ACF can be evaluated by:

(5)$\tilde{R}\left(r∆t\right)=\widehat{R}\left(r∆t\right)·u_h\left(r∆t\right)$

The normalized ACF is equal to:

(6)$\rho \left(\tau \right)={\displaystyle \frac{R\left(\tau \right)}{R\left(0\right)}}$

Equation (1) is rewritten allowing the wave spectrum estimation from tapered ACF. Moreover, the trapezoidal rule has been applied for the numerical integration of Equation (1) in the following equation.

(7)$\widehat{S}\left({\omega }_k\right)={\displaystyle \frac{∆t}{\pi }}\left(\tilde{R}\left(0\right)+2{\sum }_{r=1}^{m-1}\tilde{R}\left(r∆t\right)cos\left({\displaystyle \frac{\pi rk}{m}}\right)+\tilde{R}\left(m∆t\right)cos(k\pi )\right)$

where ${\omega }_k=k∆\omega ={\displaystyle \frac{k\pi }{m∆t}}$.

To reduce the statistical variability and smooth the spectrum obtained from Equation (7), the Moving Average and Hanning Window methods are applied separately.

Evaluation of ACF Due to Lag Number Variation

The ACF evaluates the relation between the spatial lag and the water surface elevation. It shows how strongly the surface elevation at one location is correlated to the elevation at a distance of a lag number. The evaluation of the lag number in the performance of the parameters for spectral estimation is crucial. After that, they can be adopted in the probabilistic models of wave variability such as Boccotti [36], Naess [37], and Vinje [38].

In this subsection, the statistical expectation of the ACF parameters is studied as a function of the lag variation. The number of lags is increased to evaluate this relation. Firstly, various maximum lag numbers are tested and the parameters of normalized ACF are computed, such as time lag when the first minimum of normalized ACF occurs [R(τ*)], the value of the first minimum of normalized ACF [ρ(τ*)], the value of the second derivative of the normalized ACF at the first minimum of the normalized ACF [ρ″(τ*)], as well as the value of the envelope of the normalized ACF at half of the mean wave period (r_m). The estimated values are presented in Table 1 for several lag numbers tested.

The minimum lag parameter [a = ρ(τ*)] was introduced by Boccotti [36] and represents the narrow-band parameter (Boccotti [39]) and the value of the second derivative of the normalized ACF at first (also the global) minimum, [b = ρ″(τ*)] is a parameter to be applied in the large wave height distribution model of Boccotti [36,39]. The parameter (r_m) is used for large wave heights and long-crested waves and is required for the Tayfun distribution model [40].

Various maximum lag numbers are initially tested to compute the spectra by ACF based on the Hanning Window approach, which can represent non-stationary data. Figure 4 exemplifies the spectrum obtained with the mentioned approach. The smoothed modified spectrum with a different maximum lag number is shown. Then m = 20,000 is selected based on a trade-off between calculation time and noise accuracy to derive the curve and compute the parameters. This approach represents the amplitude variability in time, known as the Hilbert Spectrum, H(ω,t). The same units of the Fourier-based variance spectra are obtained using the amplitude squared. The main difference is that the frequency is now a function of time. In this context, the Hilbert spectrum represents the frequency content, which varies over time. It is obtained through the Hilbert Transform, generally in conjunction with the Empirical Mode Decomposition (EMD) (Veltcheva and Guedes Soares [41]). It combines the instantaneous amplitudes and frequencies over time to form the Hilbert spectrum. This spectrum is a three-dimensional representation (time, frequency, amplitude) that shows how the signal energy is distributed as a function of time and frequency. However, there are some limitations, such as computational complexity, where the EMD process and the subsequent application of the Hilbert transform can be computationally intensive, and noise sensitivities, where the method may be sensitive to the noise present in the signal, affecting the accuracy of the process.

After increasing the maximum lag number applied to estimate the spectra, the resulting spectrum represents more variance in the wave amplitudes per frequency range when the lags are increased to 30,000. As the lag number increases, the signal’s noise also increases, making it impossible to fit the curve. The computed ACF, the normalized ACF, and the second derivative of normalized ACF for a lag number of 20,000 are shown in Figure 5, in which the envelope of the normalized ACF is estimated by using the mentioned Hilbert transform and shown jointly. Figure 6 presents the original spectrum obtained by ACF considering 20,000 lags jointly with the respective spectrum smoothed by the Moving Average and Hanning Window approaches. They are computed separately to compare the visual performance; however, no considerable variation is observed between the methods in the comparative visual inspection.

3.2. Fast Fourier Transform

It is possible to estimate the spectral density function by a direct Fourier transformation of the observed time series. The power spectral function is the Fourier transform of the observed series, which can be efficiently estimated by employing the Fast Fourier Transform (FFT) algorithm. It is considered a natural estimator of the spectral density function, being a raw density distribution. The spectral estimates obtained with this procedure characterize a large variance and do not present a consistent estimate of the wave parameters from the density function. Consequently, it becomes necessary to apply some smoothing techniques to reduce the variance of these raw estimates.

In this sense, the same approach of the previous subsection is used to improve the variance properties of the raw estimations, that is, the Moving Average and Hanning Window methods. Figure 7 shows the estimated spectrum obtained by FFT and the respective modified smoothed spectra. According to Borgman [42] and Priestley [43], the spectral estimates from measurements of sea surface elevation follow a chi-square distribution in which the distribution spectra of random coordinates depend on the degrees of freedom and confidence intervals assumed in the procedure.

The characteristics of the density frequencies and the associated periods, the ordinary spectral moments (asymmetry (m₃) and kurtosis (m₄) not expressed), and spectral bandwidth computed using both approaches (estimated by ACF and FFT separately) are shown in Table 2. This information consists of the mean spectral frequency ω_m, the mean period T_m, the average zero-up-crossing frequency ω_z, the zero-up-crossing wave period T_z, the peak frequency ω_p, and the peak period T_p for the analyzed series as well as information about the spectral shape by $\epsilon$ and $\vartheta$ parameters.

The last two parameters represent the width of a spectrum. They can be evaluated using two parameters: the first one ($\epsilon$) introduced by Cartwright and Longuet-Higgins in [44] and the second one ($\vartheta$) introduced by Longuet-Higgins in 1975 [45]. For a very narrow spectrum, it should be $\epsilon$ << 1. Since it involves the fourth-order moment, which may depend rather critically on the behavior of the spectrum at high frequencies, it is not convenient to use this parameter for some practical purposes (Longuet-Higgins [46]). A better estimation of the spectral width is provided by the parameter $\vartheta$, which depends only on the lower-order moments and, for a very narrow spectrum, should be $\vartheta$ << 1. In the case study, the values of the parameters are $\vartheta =0.806$ for the FFT-based spectrum and $\vartheta =0.334$ for the ACF-based, respectively. The computed bandwidth parameter corresponds to a spectrum considering 8 degrees of freedom. Figure 8 shows the final smoothed spectra estimated by both methods (ACF and FFT) for a comparative evaluation. The resulting spectra of the smoothing processes, when compared, show similar results. The differences in estimating the raw (unmodified spectrum) are evident, but both the Moving Average method and the Hanning Window present similar results when smoothing the raw spectrum.

4. Spectral Resolution and Accuracy

The spectral parameters and the spectral functions are all dependent on the approach adopted to obtain the spectrum. Moreover, the spectral shape (form of the curve) obtained in the calculations varies considerably depending on the degree of resolution used when carrying out the calculations. In the present section, the resolution of the spectral density due to the number of degrees of freedom is studied. The spectral density function of a stationary and ergodic random process η(t) is defined as follows:

(8)$S\left(f\right)={\displaystyle \frac{2}{T}}{\left|X\left(f\right)\right|}^2$

where X(f) is the Fast Fourier Transform of η(t), and S(f) represents the quantitative contribution to the process of each frequency component. The complex Fourier coefficients, X(f), can be estimated directly from the time series η(t), through the Fast Fourier Transform of this one, which for a discrete time series can be expressed as follows:

(9)$X\left(f\right)={\sum }_{-\mathrm{\infty }}^{\mathrm{\infty }}x(n∆t)e^{-i2\pi fn∆t}∆t$

Hence, considering a discrete time series constituted by N-points sampled at equidistant time intervals Δt.

$\eta \left(t_n\right)=\eta (n\Delta t)\hspace{1em}\hspace{1em}\hspace{1em}n=0,1,2,\dots ,N-1$

For a times series of duration T = NΔt, the relation between the frequency resolution, Δf and the sampling interval, Δt is given by:

(10)$∆f={\displaystyle \frac{1}{N\Delta t}}$

So, the Fourier Transform for η(t) can be written as follows:

(11)$X\left(f_m\right)={\sum }_{n=0}^{N-1}{x}_n{e}^{-i2\pi fmn/N}$

Each value of f_m is given by:

$f_m={\displaystyle \frac{m}{N\mathsf{\Delta }t}}\hspace{1em}\hspace{1em}\hspace{1em}m=0,1,2,\dots ,N-1$

However, for a sampling period, Δt, the maximum frequency that can be solved is given by the frequency of Nyquist, f_N:

$fmax={\displaystyle \frac{1}{2\mathsf{\Delta }t}}=f_N$

In other words, the maximum value of m must be N/2, so that the values of X(f) are associated with the frequencies belonging to the range.

$0\le f_m\le f_N$

In this way obtaining all the information included in the sampling interval Δt. However, this process requires a very high computing time, mainly for high values of N. Therefore, in general, these expressions are evaluated using the FFT algorithm previously mentioned, which allows for a substantial reduction in the number of multiplications involved in the resolution of the previous expressions, thus achieving a drastic reduction in the computing time required.

4.1. Smoothing the Raw Spectrum

The analysis of the statistical properties of the Fourier coefficients reveals that they present a high variance. In this way, the “raw” spectrum does not provide statistically significant results. The simplest way to reduce the variability of spectral estimates is to use the “average” operator to smooth spectral estimates. To obtain several spectral estimates for the same frequency from a discrete time series, the sequence of observations must be subdivided into smaller segments. This automatically results in a decrease in frequency resolution—which occasionally can compromise the information under complex multimodal sea states. When averaging estimates around a given frequency, they must be statistically independent of each other for the average. In practice, ensuring statistical independence may involve dividing the original time series into non-overlapping segments or using methods that generate overlapping segments but adjust for loss of independence. Therefore, there is a need for statistical independence when averaging estimates to obtain productive and meaningful results.

If the spectrum varies rapidly over these frequencies, the average will generate a significant bias. The methods used for this purpose can be classified into two groups:

1. Average of different estimates for the same frequency.

2. Average of the estimates centered around a given frequency.

Processing as in Ref. (1) corresponds to the Bartlett method [47], later improved by Welch [48]. The averages from the procedure (2) constitute the smoothing method proposed by Daniell [49].

The key advantage of Welch’s method is its ability to provide a more accurate and reliable estimation of the power spectral density by reducing variance, improving frequency resolution, and accommodating non-stationary signals through segmentation and averaging techniques. The main difference between the two methods lies in how they perform their spectral estimation. However, Welch’s method provides a trade-off between frequency resolution and variance reduction. It is more robust against variance fluctuation compared to conventional periodogram methods. In this section, only Welch’s method is presented briefly due to the mentioned advantages, as it is used in the present study to compute the smoothed spectrum from several degrees of freedom.

Welch Method

Welch (1967) [48] proposed a modification of Bartlett’s method, which consists of subdividing the original series η(t) into K-segments of M-data each, then applying a window, w(t), permitting such segments to overlap. The objective of using the data setup windows is to reduce the bias of the estimates, minimizing the “leakage” effect, although this causes a slight decrease in the frequency resolution. Overlapping segments allow an increase in the number of segments to be averaged. With the application of these improvements, together with the computational efficiency of the FFT algorithm, the Welch method became the most used spectral estimation procedure for a given degree of freedom today. The Welch method can be described as follows. Given a record of N-data points, η(t):

$\{\left(0\right),\left(1\right),\dots ,\left(N-1\right)\}$

By dividing η(t) into K-segments of M-data, thus overlapping by an amount (Seg = M − D), where D represents the displacement between adjacent segments, the beginning of the second sequence is located in the ordinate D, the third sequence in 2D, and so on, until reaching the position N−M₋₁. In this way, K-segments can be expressed as follows:

(12)${\left(I\right)}_j=(I-JD)$

where

(13)$I=\mathrm{0,1},2,\dots M-1;\hspace{1em}J=\mathrm{0,1},2,\dots K-1$

In this case, JD is the origin of the J-sequence. That is the difference between the amount of data per segment and that of overlapping data.

(14)$JD=J(M-Seg)$

For an M and a given data, the number of subseries that is obtained from a record of N-data is given by:

(15)$K=INT\left[{\displaystyle \frac{N-Seg}{M-Seg}}\right]$

and the estimation of the spectrum corresponding to each of them is expressed as follows:

(16)$\widehat{S}j\left(f\right)={{\displaystyle \frac{2}{\mu M∆t}}\left|Xj\left(f\right)\right|}^2\hspace{1em}\hspace{1em}\hspace{1em}0\le f\le {\displaystyle \frac{1}{2∆t}}$

where X_j (f) is the discrete Fourier transform of the J^th segment,

(17)$X_j\left(f\right)=∆t{\sum }_{I=0}^{M-1}\omega \left(I\right)J\left(I\right)e^{(-i2\pi fI∆t)}$

The $u$ is a correction factor introduced to consider the energy (variance) reduction in the subseries when applying a data setup window. It can be estimated by:

(18)$u={\displaystyle \frac{1}{M}}{\sum }_{I=0}^{M-1}({\omega \left(I\right))}^2$

The used data setup window in the Welch procedure is the truncated cosine window expressed as follows:

(19)$\begin{array}{ll}{\omega }_j={\displaystyle \frac{1}{2}}\left[1-cos\left({\displaystyle \frac{\pi j}{k}}\right)\right]& j=0,\dots ,k-1\\ {\omega }_j=1& j=k,\dots ,M-k-1\\ {\omega }_j={\displaystyle \frac{1}{2}}\left[1-cos\left({\displaystyle \frac{\pi (M-j)}{k}}\right)\right]& j=M-k,\dots ,M-1\\ k=0.1M& \end{array}$

After calculating the periodograms corresponding to each of the K-subseries, the average of these is determined to obtain the spectral estimator of Welch.

(20)$\tilde{S}\left(f\right)={\displaystyle \frac{1}{K}}{\sum }_{j=0}^{k-1}\widehat{S}j\left(f\right)$

For fixed values of N and M, the total reduction in variance for Seg = 50% is approximated by:

(21)$Var\left[\tilde{S}\left(f\right)\right]\approx {\displaystyle \frac{11}{18K}}{\sigma }^2\left(f\right)$

The procedure of overlapping the data segments is used to cause a decrease in the variance of the estimates. In this sense, the amount of overlap (Seg) used in the analysis is usually expressed as a percentage of M.

(22)$Seg={\displaystyle \frac{X}{100}}·M$

This represents one of the most conflicting aspects of this method. Welch [48] recommends that this should not exceed 50%. However, some authors consider an overlap between 50% and 75% (e.g., Kay [50]), while others reject such a high percentage (e.g., Yuen and Fraser [51]). The present work opted for an overlap percentage of 50%. Typically, when a given spectrum results from smoothing functions such as Welch, some of these functions affect the degree of resolution induced by the degree of freedom considered. The variation in the degrees of freedom modifies the shape of the smoothed spectra, including over-smoothed or under-smoothed peaks found in the spectra. Chakrabarti and Cooley [23] found that the spectral parameters and average periods show large variations with different values of frequency and degrees of freedom used when calculating the spectra. They recommended that any statistical parameters obtained from wave spectral density distribution should appear with their associated degrees of freedom. The number of degrees of freedom of the spectral obtained by this method can be expressed as a function of the percentage of overlap, such as:

(23)$\upsilon =2K=2\left[{\displaystyle \frac{N-Seg}{M-Seg}}+1\right]$

From a Welch function over the raw spectrum, the resolution of the calculated function is tested due to the change in degrees of freedom. The degrees of freedom to test the variance spectra are 8, 16, 32, 64, and 128. The spectra corresponding to each degree of freedom are calculated and reported (see Table 3 and Figure 9). In Table 3, it is important to differentiate the concept when calculating significant wave height. Basically, H_m (Significant Wave Height) is a generic term that can refer to various definitions of significant wave height, depending on the calculation method used. The most common is the significant wave height (H_1/3), which is calculated as the average of the upper third of wave heights in a record, usually denoted as H_s. Still, here, it is just denoted as H_m. Otherwise, H_m0 (Significant Spectral Wave Height) refers to the significant wave height derived from the wave energy spectrum calculated from the area under the wave energy density spectrum S(f) curve using the formula: H_m0 = 4√m₀, where m₀ is the zeroth moment of the spectrum. Additionally, in the context of use, H_m is often used in practical and operational contexts where one deals directly with observed wave heights. On the other hand, H_m0 is used in contexts where spectral analysis is relevant, such as in wave modeling and wave energy studies. Both measurements are important to characterize the state of the sea, but H_m0 offers a view based on energy distribution, which may be more informative in specific scientific and engineering contexts.

The number of degrees of freedom used induces the process of smoothing the raw signal. For curves established with several degrees of freedom, such as 128, 64, or even 32 degrees of freedom, the shape of the resulting curves seems too smooth, and then information regarding density variation per frequency is lost, as illustrated in Figure 9(top). Moreover, it is possible to appreciate that the energy density contained is spread over two peaks, one with the highest intensity located in the low-frequency zone (occurring at a frequency around 0.09 H_z), and the other secondary peak located in the highest frequency zone, being computed over the frequency of 0.115 H_z. The figure presents a bimodal spectrum that generally represents two different wave systems. When the degrees of freedom are limited and reduced to 16 or 8, the curves seem to better represent the raw spectrum by following the bimodal trend, as shown in Figure 9(bottom). Attention is paid to the small number of degrees of freedom presenting a bimodal behavior.

4.2. Fitting Smoothed Spectrum

The general procedure used to characterize the variance spectral density function of the wind wave is based on existing knowledge about the physical mechanisms processes applied to fitting experimental observations. It is required to define certain analytical functions capable of describing the distribution of energy content along the frequency. Over the last few years, some formulations have generally been accepted to characterize the wave spectrum in terms of a few parameters to represent wave conditions and to deal with the representation of the design conditions of marine structures. The accepted models are the fully developed wave spectral model of Pierson and Moskowitz [52] and the developing wave spectrum JONSWAP (Hasselmann et al. [53]), the latter being widely used during the last three decades. The JONSWAP is the generic single peak wave spectrum as it reduces to the Pierson Moskowitz when the peak enhancement parameter is equal to one. In this section, the model used to fit the smoothed spectrum is the JONSWAP spectral model, particularly the wind speed parameterization by Goda [54].

The spectra measured during the JONSWAP are suggested to characterize developing sea states while Pierson Moskowitz described fully developed waves. It is a spectral formulation similar to the one proposed by Pierson and Moskowitz [52], modified by introducing an additional term that allows the intensification of the peak on the spectra. This term is called the spectral peak enhancement factor. The JONSWAP model can be represented by Hasselmann et al. [53]:

(24)$E\left(f\right)=\left(\alpha g^2{\left(2\pi \right)}^{-4}{f}^{-5}\right){\gamma }^{exp\left({\displaystyle \frac{-{(f-f_p)}^2}{{2\sigma }^2{f}^2}}\right)}$

where $\sigma =\left\{{\sigma }_a for f\le f_p\right\};\left\{{\sigma }_b for f>f_p\right\}$.

E(f) is the variance energy spectrum and f is frequency. The previous equation contains five spectral parameters that define the spectral curve. The parameters f_p and α are scale parameters, representing the frequency at the maximum of the spectrum and the Phillips constant (Phillips [55]), respectively. The remaining three parameters are crucial to define the shape of the curve of the spectrum: the mentioned peak enhancement factor γ and σa and σb define the left–right side widths of the spectral peak region, respectively. The described models are parameterized based on wind speed, in which Goda [54] suggests a modification in the expression of the JONSWAP spectrum to re-define the expression in terms only of the significant wave height, H_m0 and the corresponding spectral peak frequency, T_p. The functional form can be expressed as follows:

(25)$S\left(f\right)=\left(\beta {H_{m0}}^2{T_p}^{-4}{f}^{-5}\right){\gamma }^{exp{[-\left(T_{p}f-1\right)}^2}/{2\sigma }^2$

where $\beta =\left({\displaystyle \frac{0.06238}{0.230+0.0336\gamma -0.185{(1.9+\gamma )}^{-1}}}\right)\left[1.094-0.01915ln\gamma \right]$.

The parameter $\beta$ is a dimensionless constant being a function of enhancement factor (γ), and T_p denotes the period corresponding to the spectral peak frequency. Remember that the peak enhancement factor (γ) depends on the stage of wave development and can vary between 1 and 7.

4.2.1. Fitting 8 and 16 Degrees of Freedom

In Welch’s method for spectral analysis, the term “degrees of freedom” refers to the number of data segments used in computing the power spectral density estimate. As mentioned, the method involves dividing a signal into overlapping segments, applying a window function to each segment, computing the Fourier Transform of each windowed segment, and averaging the squared magnitude of these transforms to obtain the spectral estimate. So, the number of degrees of freedom in Welch’s method is determined by the number of segments used in the computation. In practice, choosing the appropriate number of degrees of freedom (or segments) involves a trade-off between frequency resolution and variance reduction.

Increasing the number of segments decreases frequency resolution but improves variance reduction. Therefore, understanding the relationship between the number of segments, frequency resolution, and variance reduction in spectral analysis is essential for the correct formulation. Conversely, decreasing the number of segments reduces frequency resolution but increases variance reduction. Thus, many degrees of freedom appear to lead to a white noise signal in which it has no physical meaning. In this vein, as the idea is to study uncertainties and reduce the variance of estimates, the degrees of freedom evaluated are limited to 8 and 16. Some crucial aspects to consider, such as decreasing frequency resolution to increase the number of segments, typically involve using shorter segments, and the ability to resolve closely spaced frequencies diminishes. Importantly, to improve variance reduction, more segments mean more independent estimates contributing to the average, thus reducing the overall variance of the spectral estimate. So, by considering these aspects, analysts can make informed decisions and optimize their spectral analysis procedures to suit their data’s specific characteristics and requirements. Moreover, there is no universally “best” or fixed number of degrees of freedom for estimating the Welch spectrum, as it depends on various factors such as the characteristics of the signal, the desired frequency resolution, the length of the signal, and the trade-offs between resolution and variance. Engineers and researchers commonly experiment with different numbers of segments (degrees of freedom) to strike a balance between achieving sufficient frequency resolution and reducing variance in the spectral estimation. The choice often involves empirical considerations based on the specific characteristics of the signal being analyzed and the requirements of the analysis.

The present subsection is dedicated to fitting the Welch spectra smoothed by using 8 and 16 degrees of freedom (DoF), respectively. The model used is the JONSWAP modified by Goda, as described in the previous section. Table 4 shows the results of the parametric fit for spectra computed with 8 DoF and the confidence intervals for 99%. The goodness of fit is also calculated and based on robust statistical metrics. Figure 10 (top panel) illustrates the data fitting process by the JONSWAP model dependent on the H_s and T_p and a slight visual improvement in confidence when intervals based on 99% are included.

Table 5 indicates the parameters resulting from the fitting process for 16 DoF. The fit is slightly better than the previous 8 DoF spectrum as noted in the goodness of fit statistic. Visually, the fitted curve can reach the main peak; however, the small secondary peak (but relevant in energy terms) is underestimated. This bias in fitting the peak improves when the confidence interval is considered for future predictions. The fitted curve is computed considering 99% of boundaries, thus achieving robust confidence intervals. Figure 11 shows the fitted curve to spectral density corresponding to 16 DoF and the computed confidence intervals.

The previously studied spectra are composed of 8 and 16 DoF and appear slightly bimodal in which two peaks are apparent when the spectral density data are plotted. However, the JONSWAP model used to fit the data has only one peak of the spectrum, even though increasing the confidence intervals up to 99% the secondary peak is excluded from the boundaries. In this way, it is evident that a theoretical formulation capable of characterizing the energetic content in the complete range of frequencies is needed. Naturally, the observed wave structures for 8 and 16 DoF are satisfactorily fitted but cannot be described by a simple peak spectral model.

The progressive development of wave measurement techniques allowed many authors to propose different empirical formulations that are functions of the degree of wave development to a given sea system, such as fully or not fully wave development conditions. The formulations differ in aspects such as the number of free parameters to fit, a consequence of different mathematical assumptions when the variable is parametrized, normally specified by sea-state dependence or function of external variables such as wind speed or fetch length.

Occasionally, wave spectra exhibit three defined peaks, one associated with the local wave field and another corresponding to swell waves with different energy contents, as Boukhanovsky and Guedes Soares [56] presented. Thus, to better describe the energy content defined by the spectrum density, the model requires flexibility enough to represent the density distribution acceptably, in particular when two or three sea systems are computed. The double peak recorded in the spectrum characterizes a mixed-sea situation and requires extended models to fit such a double peak, such as the model proposed by Guedes Soares [57] that builds upon two JONSWAP spectra. This model could be considered the generic one for a specific value of one parameter; it is reduced to the single peak JONSWAP model, which in turn includes the Pierson-Moskowitz.

4.2.2. Fitting Double-Peaked Spectra

Regularly, the wave spectra in the ocean are multipeaked peaked due to remarkable dependence on meteorological and oceanographic factors. In particular, local wind waves often develop in the presence of some background low-frequency swells from distant storms [58]. Titov [59] showed that wind waves and swells for the North Atlantic during winter and summer prove that zones with heavy swells extended over large open oceans and coastal areas are typical features of all oceans. Thompson [60] analyzed wave records in different locations, including the North Atlantic Ocean, and observed that multi-peaked spectra are common at all locations. Cumming et al. [61], using hindcast data, determined in the North Atlantic Ocean that 25% of spectra are double-peak. Guedes Soares [57] used measured data from the North Atlantic to assess the probability of occurrence of double peak spectra and observed in global situations the percentage of occurrence of bimodal spectra is about 20–25%, and Guedes Soares and Nolasco [62] conducted a similar assessment in offshore zones in Portugal and obtained 23–26%. Boukhanovsky and Guedes Soares [56] and Lucas et al. [7] have analyzed even more complex spectral shapes and obtained comparable percentages in offshore and coastal waters.

All previously mentioned justifies why it is not possible to establish a universal spectral model that is validated for all random conditions. In other words, while the wave spectrum is a unique description of the wave system generated by a defined wind condition, the wave spectra that exit at different points in the ocean have a considerable variation as a function of time for a given place, resulting from varying wind conditions and combinations of different wave systems. Furthermore, an aspect less mentioned is the existence of ocean currents that will change the shape of the wave spectrum by significant amounts when they reach a considerable speed (Ponce de Leon and Guedes Soares [63,64]). In short, it is not possible to use a single formulation to characterize the different energetic structures of the waves, even at the same point in the ocean.

Entirely the spectral formulations proposed to represent bimodal spectra are based on the idea of decomposing the spectrum into two parts: one associated with the spectrum in the low frequencies area and another one for the higher frequencies, or what is equivalent, characterize the bimodal spectra by the linear superposition of two unimodal spectra (Guedes Soares [57]):

(26)$S\left(f\right)={S\left(f\right)}_{sw}+{S\left(f\right)}_{ws}$

where ${S\left(f\right)}_{sw}$ indicates the spectrum corresponding to the low-frequency band (swell) and ${S\left(f\right)}_{ws}$ represents the spectrum associated with the high-frequency wave system (wind-sea). Once the spectra are computed, a spectral formula needs to be attributed to each component and combined to reproduce the total wave spectrum. In this work, the spectral components are fitted by a JONSWAP model to each wave system component, which allows the individual spectral parameters to be compared. Aiming to obtain a good spectral description of the bimodal wave system, the peaks are then fitted with a modified JONSWAP spectral function, and the summation of these is used to describe the whole frequency spectrum.

Rodriguez and Guedes Soares [14] have suggested an approach for the automatic identification of multimodal sea wave spectra. Having identified the peak frequencies of the spectral components, a JONSWAP model can be fitted around each of the component’s peak frequencies, as Guedes Soares and Carvalho [65] described. This scheme can be applied to any number of peaks defining different sea-state components. Figure 12 and Figure 13 show the observed spectrum estimated using 8 degrees of freedom and 16 degrees of freedom, respectively. Table 6 describes the parameters fitted and the goodness of fit parameters.

Considering the metrics used in the goodness of fit, the spectrum with 16 degrees of freedom performs a little better. Moreover, the approximation of the empirical model with the observed data is considered very good, indicating robustness and confidence in the computed output parameters. Great flexibility is also appreciated in the bimodal model used, which can satisfactorily fit the two representative peaks in the spectrum with 8 degrees of freedom, where the highest values for the peak enhancement factors are computed.

Torsethaugen [66] used the model of two JONSWAP formulations proposed by Guedes Soares [57], and for each sea state characterized by H_m0 and T_p, he provided the values of the spectral parameters that resulted from fitting the JONSWAP model to actual ocean data. In addition, Torsethaugen proposed a criterion for the classification of the swell-dominated sea and wind-dominated sea for both spectra:

$\begin{array}{l}\mathrm{Swell}:T_p>T_f\\ \text{Wind-sea}:T_p\le T_f\\ \mathrm{and} T_f=a_f{\displaystyle \frac{H_{m0}}{3}}\end{array}$

where T_p is the peak period, H_m0 is the significant wave height calculated from the zeroth moment of the spectrum, and a_f = 6.6 is adapted from the JONSWAP experiment (Hasselmann et al. [53]). Logically, if T_p ≤ T_f, the local wind-sea dominates the spectral peak, if T_p > T_f the swell dominates the spectral peak. In this way, considering the values of H_m0 exposed in Table 3 estimated to degrees of freedom 8 and 16 that are computed by H_m0 = 4√m₀.

Tf_8DoF = 6.6 (3.54)^1/3 = 7.78 where the T_p = 10.94, swell dominates the spectral peak.

Tf_16DoF = 6.6 (3.63)^1/3 = 7.98 where the T_p = 10.79, swell dominates the spectral peak.

This behavior is also verified visually when the main observed peak is located at lower frequencies compared to the secondary peak’s frequency.

5. Robust Estimation of Representative Spectra

After evaluating the bimodal behavior present in spectra with small numbers of degrees of freedom, a robust resulting spectrum is computed. The robust spectrum is the spectrum that results from the average of all 5 degrees of freedom considered. The spectral densities present in every frequency band are averaged, resulting in a unique representative spectrum. Figure 14 shows all the degrees of freedom applied in the average process and the resulting averaged robust spectrum overlaid with the raw spectrum.

The energy content, the frequency distribution, and the shape of the robust spectrum resulting from averaging several spectra over the frequency band reveal and reinforce the existence of the bimodal condition since the resulting spectrum presents two peaks. It can be related to atmospheric variability, which generates different types of wave systems in the study place. The main idea of the process is to define a robust way to estimate the representative spectrum for a location. It can be used whenever there are uncertainties in estimating the climatic spectrum concerning the number of degrees of freedom considered. Considerable variability is observed between spectra derived with varying degrees of freedom, as shown in the previous figures.

This variability plays an extremely important role in practical cases in which only one spectrum is used. The average of the densities measurements for each frequency band is suggested to reduce this variability and provide robustness in estimating the representative climatic spectrum. So, for example, at a frequency of 0.2 H_z, five density measurements are considered, not just one. Additionally, the ensuing average of these five densities results in a value more robust for this specific frequency than one. This understanding is extended over the entire frequency range and then the resulting spectrum is obtained.

In this way, when obtaining five measurements of spectral density for each frequency band, the robustness is increased since the final spectrum is defined by five measurements and not by a single degree of freedom, as used in the practical cases where only one spectrum (generally derived using 16, 32, or 64 DoF) is considered to defining the sea states.

From the set of five spectra corresponding to wave spectra computed with five different numbers of degrees of freedom, a maximum and minimum value of the spectra for each frequency band is obtained. Naturally, maximum and minimum spectra are heavily affected by outliers, and their representativeness is very low. It happens because the envelope curve is defined by just one spectral estimation for every frequency band. The spectrum resulting from averaging five spectra from various DoF is fitted to determine the robust estimator.

The Nonlinear Least Squares method is applied to fit the spectral curves. This method seeks to minimize the sum of squares of the differences between the observed values and those predicted by the fitted model. The standard Levenberg–Marquardt approach is used for nonlinear least squares problems. It combines the Gauss–Newton method and gradient descent approaches, providing robust and efficient convergence even when the model is highly nonlinear. The mentioned methods are used to compute the parameters from Table 7. Table 7 presents the computed fitted parameters for the averaged final spectrum and the confidence intervals for 99%. Figure 15 shows the final spectrum fitted using the JONSWAP model, which is a function of sea-state parameters such as H_s and T_p. It embodies a summation of two typical JONSWAP (three parameters), each of them representing the peak distributions over the frequencies and consequently increasing the complexity of the fit process since the representation of two peaked curves requires six parameters. Nevertheless, despite the complexity, the model used manages to represent the bimodal curve excellently.

Figure 15 shows the final spectrum and the fit curve computed by the model. On visual inspection, it is possible to see that the model has enough flexibility to fit the two peaks observed in the averaged final spectrum. Furthermore, the good performance of the spectral model in representing the computed spectrum is confirmed by the reasonable values estimated by the metrics, such as sum squared error performance function (SSE), coefficient of determination (R²), adjusted coefficient of determination (Adj. R²), and root mean square level (RMSE) that are briefly explained following.

The Sum Squared Error Performance Function (SSE) measures the total deviation in the predicted values from the empirical values. It sums the squares of the differences between predicted and empirical observations. It indicates the overall error of a model, so lower SSE values indicate a better fit. The Coefficient of Determination (R²) is a statistical measure that explains the proportion of variance in the dependent variable that is predictable from the independent variable(s). It assesses the goodness of fit of a model. Their values range from 0 to 1, with higher values indicating a better fit. Adjusted Coefficient of Determination (Adj. R²) modifies the R² value to account for the number of predictors in the model. Unlike R², it can decrease if unnecessary predictors are added. So, it provides a more accurate measure of model fit by adjusting for the number of predictors, preventing overfitting. Finally, the Root Mean Square Error (RMSE) calculates the square root of the average of the squared differences between predicted and actual values. It measures the standard deviation in prediction errors, providing insights into the model’s prediction accuracy. In this vein, lower RMSE values indicate better predictive performance. Such metrics were used by Guedes Soares [57,67] to represent double-peaked sea wave spectra and evaluate the occurrence of double-peaked wave spectra. They were also used by Rodriguez and Guedes Soares [14,15] to determine the uncertainty in estimating the slope of the high-frequency tail of wave spectra and in the criteria for the automatic identification of multimodal sea wave spectra. The references cited above offer a detailed mathematical formulation of the metrics used.

However, the structure of the spectrum that is not located at the peak, more precisely in high-frequency regions, the model curve does not fit perfectly due to the small fluctuations in the energy distributions in this frequency band. Nevertheless, this apparent lack of flexibility in these cases can be disregarded when the confidence intervals are defined and displayed jointly. Figure 16 illustrates the previously commented, where the small fluctuations in energy density observed in the high-frequency bands are covered by the 99% confidence intervals, further increasing the robustness of the computed parameters.

Moreover, a graphical display of the residuals for a JONSWAP model averaged spectrum fit is shown in Figure 17. Assuming the model is applied to fit the data, the residuals approximate the random errors. Therefore, if the residuals appear to behave randomly, it suggests that the model fits the data well. However, if the residuals display a systematic pattern, it is a clear sign that the model fits the data poorly. The plot shows that the residuals are calculated as the vertical distance from the data point to the fitted curve and displays the residuals relative to the fit, which is the zero line. The residuals appear well scattered around zero in the frequencies where the peaks are observed, indicating that the model describes the data reasonably well. Finally, it is important to highlight that the exposed method has limitations in applications that evaluate the directional spectrum. However, it could also be limited to normal non-extreme regime applications, as it was evaluated and implemented for a specific extreme event considering a single storm.

6. Discussion

Considering 20 min records of real oceanic wave data, focusing on the sea-surface elevation on 21 August 2009 (when a storm occurred), the autocorrelation function and the wave spectrum were calculated to estimate the parameters of the sea states. Then, assuming the well-known Welch method for estimating the wave spectrum, various degrees of freedom (DoF) applied to estimate the spectrum are compared. The DoF fundamentally accepts the same understanding as the degree of smoothing of the raw signal. The main idea of the technique is to segment the signal into many subsections, reducing the bias due to the leakage effect. The original time series is segmented into K-blocks of equal length M, which can be overlapped in a given fraction (Seg), usually close to or lower than 50% of the segment length, and for this work, the overlapped windowed segment is averaged in 50% of the length. The frequency resolution tested using different numbers of DoF applied to smooth the raw spectrum is performed successfully (DoF = 8, 16, 32, 64, and 128). As expected, each estimated spectrum presented a different shape from the other and also presented important variability when compared.

The degrees of freedom induce the stability of the model in the smoothing process. It is observed that as the DoF is increased, practically the signal is getting over-smoothed. The smallest DoF applied is 8, and despite having few degrees of freedom, the model can coherently smooth the raw signal. Then DoF = 16, which is also considered a small number of degrees of freedom, results in a coherent smoothing process. The two smoothed spectra by 8 and 16 were selected to perform theoretical density curve fit by a suitable spectral model. Through this evaluation, it is possible to determine the ordinary spectral moments (m_i) used to compute wave height parameters. It also helps to assess the frequencies and characterize the associated periods. The evaluation of the spectral bandwidth is verified by considering the Cartwright and Longuet-Higgins [44] approach and the results of Longuet-Higgins [45]. For most practical cases, simple checking of spectral width results can be performed considering that it is a narrow band.

In the spectral evaluation, a crucial subsection is added in which it considers the variability in the spectral resolution. Various degrees of smoothing are tested to assess the accuracy. Empirical spectral models are applied to fit the spectra curves and estimate the free parameters. When the resulting spectrum has only one peak, the traditional JONSWAP model is considered in the curve fitting. However, in the spectra with few degrees of freedom, the curve has two peaks (although the second peak appears faint, two wave systems are noticeable), so the simplified JONSWAP model did not perform well when fitting these two peaks. Thus, a modified model is required to test the curve fits. The model considered is the JONSWAP modified by Goda, which presented an excellent performance when fitting the two peaks observed in the spectrum smoothed with 8 and 16 DoF. Either of these two spectra can be considered representative of determining wave parameters. However, this type of uncertainty dealt with in spectral resolution when smoothing the spectrum is nowadays considered a challenge. Due to the variability computed between the different smoothing degrees, it is difficult to decide which one to use as the representative sea-state spectrum for a given location. With this in mind, a methodology considered robust contemplates all spectra previously calculated in the raw spectrum smoothing process (five spectra in total). In summary, for each and all frequency bands, the contained energy densities were computed using the classic central tendency statistical estimator: the average operator. As a consequence, a final spectrum was generated as a result of average estimations.

Finally, the final spectrum can be considered valid, alternative, and robust, showing a bimodal behavior with two peaks. Thus, the spectral density curve fitting process is performed, and the modified JONSWAP two-peaked model is considered in this process. The final spectrum curve that averaged the densities of five spectra (one curve per DoF) for each frequency was excellently fitted by the JONSWAP empirical model. The metrics used in the goodness-of-fit test demonstrate that the process was completed, and the 99% confidence in the parameters is considered truthful. So, the robust alternative is applying the average and assuming the averaged spectrum as the most accurate to consider the climatic wave for a given location. It is a good indicator of central tendency because it is fairly robust and stable under various external factors.

7. Conclusions

Evaluating actual recordings of the sea surface enables both theoretical and practical applications. For example, the ACF approach, parameters such as time lag, and the first (global) minimum of normalized function are parameters used in Boccotti’s [36] models and adopted by Naess [37] and Vinje [38] to parameterize their wave height distribution models. Furthermore, the values of the normalized ACF envelopes provide important information that makes it possible to compute characteristic parameters of the sea states. This parameter is used in the distribution of Tayfun [40] for large wave heights. The value of the second derivative of the normalized ACF is important and required whenever the Boccotti distribution model for wave heights [36,39] is used in the modeling. In summary, the Autocorrelation Function provides insights into the correlation between a signal and its delayed versions, allowing the identification of repeating patterns and periodicities within the signal. It can estimate the dominant frequencies by analyzing the correlation at different lag times. However, its computational complexity can be relatively high for longer data sets, making it less efficient for large-scale analysis. Nevertheless, for short-term investigations, it is robust and comparable with the Fast Fourier Transform method to model extreme mixed sea conditions generated by the passage of a storm.

On the other hand, the Fast Fourier Transform efficiently converts a signal from its time-domain representation into the frequency domain, clearly visualizing the signal’s frequency components and their respective amplitudes. FFT is highly efficient computationally, making it particularly advantageous for processing large datasets quickly. It simplifies the identification of specific frequency components, which can be directly comparable (see Table 3). Often, a combination of these techniques might be employed to obtain a comprehensive understanding of the wave spectrum, leveraging the strengths of each method.

Finally, the choice between ACF and FFT in determining wave spectra depends on the specific characteristics of the data, the level of temporal information required, the computational resources available, and the depth of frequency analysis needed for the study or application at hand. Integrating both methods or choosing one over the other should be based on carefully assessing these factors to derive the most accurate and informative wave spectrum analysis. Thus, with such factors considered, the study addressed the uncertainties generated by the choice in the number of degrees of freedom when smoothing the spectrum using the classical Welch method. So, the most common uncertainties in the estimation of wave parameters from the Welch spectral method are treated here. The variability in the resulting spectrum due to the degrees of freedom of resolution is clear, generating uncertainties in the final results.

This uncertainty affects the estimation of wave parameters directly from the spectrum, which can be reduced considering the simple and robust methodology presented in the current study. Finally, a simple approach that averages the densities of five curves for each frequency band reducing the uncertainty is recommended. This generates a final spectrum, which has been fitted with an ideal empirical model. In this way, the questions on which spectral resolution to use in obtaining the representative spectrum have been discussed and uncertainties have been reduced, which allows establishing a robust final wave spectrum.

Author Contributions

Conceptualization, R.M.C.; Methodology, G.C. and R.M.C.; Validation, G.C.; Formal analysis, G.C. and R.M.C.; Resources, R.M.C.; Data curation, G.C.; Writing—original draft, G.C.; Writing—review & editing, R.M.C. and C.G.S.; Visualization, G.C. and R.M.C.; Supervision, G.C., R.M.C. and C.G.S.; Project administration, C.G.S.; Funding acquisition, R.M.C. and C.G.S. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The buoy data can be obtained in the links (accessed on 4 October 2022): https://www.marinha.mil.br/chm/dados-do-goos-brasil/pnboia-mapa and http://www.goosbrasil.org/pnboia/dados/.

Acknowledgments

The authors would like to acknowledge the Brazilian Navy and the program PNBOIA for providing the buoy data, as well as the Oceanographic Observation and Modelling Network (REMO). This work contributes to the Strategic Research Plan of the Centre for Marine Technology and Ocean Engineering (CENTEC), which is financed by FCT under contract UIDB/UIDP/00134/2020. The authors would like to express their gratitude to Professor Carlos Eduardo Parente Ribeiro for his valuable support with spectral analysis and signal processing, including both theoretical statistical guidance and assistance with coding.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Location of the moored buoy at the South Coast of Brazil and the Axys 3M data series buoy.

Figure 2. Significant wave height recorded by the buoy during the passage of a storm on 20 and 21 August 2009.

Figure 3. Illustration of sea surface displacement on 21 August 2009 as well the number of waves recorded by the heave movement of the buoy.

Figure 4. The wave spectra estimated by ACF exemplified by various lag numbers.

Figure 5. Functions and parameters related to the autocorrelations.

Figure 6. Wave spectrum obtained from ACF (original spectrum) and the smoothed spectrum by Hanning Window and Moving Average approaches.

Figure 7. Wave spectrum obtained from FFT (original spectrum) and the smoothed spectrum computed by Hanning Window and Moving Average.

Figure 8. Comparison of wave spectra obtained from ACF and FFT.

Figure 9. (a) Smoothed spectra computed by the Welch method for several degrees of freedom. (b) Attention to the small number of degrees of freedom to assess the bimodal system.

Figure 10. Top panel: 8 degrees of freedom spectra fitted by the JONSWAP model. Bottom panel: 99% of confidence interval (predictions boundaries).

Figure 11. Top panel: 16 degrees of freedom spectra fitted by the JONSWAP model. Bottom panel: 99% of confidence interval.

Figure 12. Welch spectra for 8 degrees of freedom as well as the JONSWAP bimodal fit.

Figure 13. Welch spectra for 16 degrees of freedom as well as the JONSWAP bimodal fit.

Figure 14. The wave spectrum for various degrees of freedom and the robust spectra resulting from average density over frequencies, respectively.

Figure 15. The fitted final spectrum resulting from averaging several DoF.

Figure 16. The fit process considering the confidence intervals in order to cover all frequency bands.

Figure 17. Residual plot of fitting the JONSWAP spectral model data.

Table 1

Parameters from the autocorrelation function computed for several lag numbers.

Lag Number	R(τ*)	ρ(τ*)	ρ″(τ*)	rm
300	3.53	−0.343	0.246	0.307
500	5.01	−0.177	−0.021	0.052
1000	3.75	−0.425	0.236	0.409
2000	4.11	−0.579	0.239	0.587
3000	4.22	−0.635	0.247	0.644
10,000	4.24	−0.643	0.248	0.652
20,000	4.24	−0.645	0.249	0.655
30,000	4.24	−0.646	0.249	0.655

Table 2

Spectral frequencies and associated periods, moments, and bandwidth computed based on ACF and FFT.

Parameters	ACF	FFT
Zero spectral Moment $m_0$	0.830	0.848
First spectral Moment $m_1$	0.680	0.705
Second spectral Moment $m_2$	0.620	0.966
${\omega }_m={\displaystyle \frac{m_1}{m_0}} \left({\displaystyle \frac{rad}{s}}\right)$	0.820	0.830
$T_m={\displaystyle \frac{2\pi }{{\omega }_m}} (s)$	7.662	7.563
${\omega }_z=\sqrt{{\displaystyle \frac{m_2}{m_0}}} \left({\displaystyle \frac{rad}{s}}\right)$	0.864	1.067
$T_z={\displaystyle \frac{2\pi }{{\omega }_z}} (s)$	7.266	5.888
${\omega }_p\left({\displaystyle \frac{rad}{s}}\right)$	0.722	0.718
$T_p={\displaystyle \frac{2\pi }{{\omega }_p}} (s)$	8.695	8.746
$\epsilon =\sqrt{1-{\displaystyle \frac{{m_2}^2}{m_0{m}_4}}}$	1.00	1.00
$\vartheta =\sqrt{{\displaystyle \frac{m_0{m}_2}{{m_1}^2}}}-1$	0.806	0.5334

Table 3

Wave statistical parameters computed for different degrees of freedom.

DoF	m 0	m 1	m 2	m 3	Tm	Tp	Hm	Hm 0	$\mathit{\vartheta }$
8	0.781	0.104	0.016	0.0035	7.492	10.941	2.219	3.545	0.459
16	0.823	0.108	0.017	0.0036	7.595	10.793	2.277	3.639	0.464
32	0.842	0.110	0.017	0.0037	7.618	9.984	2.304	3.681	0.472
64	0.882	0.115	0.018	0.0039	7.627	9.623	2.358	3.767	0.490
128	0.936	0.117	0.020	0.0045	8.002	9.180	2.428	3.879	0.614

Table 4

The 8 degrees of freedom parametric fit for the JONSWAP model and the goodness of fit.

Coefficients (99% Confidence Bounds)			Goodness of Fit
Parameters	Values	Intervals	R2	0.8196
β	1.596	(1.4605, 1.9366)	Adj R2	0.8161
γ	1.944	(1.6855, 2.2033)	SSE	737.2
σ	0.053	(0.0376, 0.0696)	RMSE	1.892

Table 5

The 16 degrees of freedom parametric fit for the JONSWAP model and the goodness of fit.

Coefficients (99% Confidence Bounds)			Goodness of Fit
Parameters	Values	Intervals	R2	0.9499
β	1.121	(0.4223, 3.9357)	Adj R2	0.9490
γ	1.940	(1.6890, 2.1901)	SSE	199.5
σ	0.169	(0.1696, 0.2053)	RMSE	0.9396

Table 6

Statistical fitting parameters for 8 and 16 degrees of freedom spectra. f_p = 1/T_p.

8 Degrees of Freedom				Goodness of Fit
Parameters	Values	Parameters	Values
β 1	0.690	β 2	0.547
γ 1	9.722	γ 2	3.953	R2	0.9633
σ 1	0.091	σ 2	0.035	Adj R2	0.9617
Hs 1	1.626	Hs 2	1.232	SSE	150
Tp 1	10.730	Tp 2	8.701	RMSE	0.8638
fp 1	0.093	fp 2	0.114
16 Degrees of Freedom				Goodness of Fit
Parameters	Values	Parameters	Values
β 1	1.306	β 2	0.011
γ 1	2.332	γ 2	3.292	R2	0.9872
σ 1	0.062	σ 2	0.131	Adj R2	0.9817
Hs 1	1.217	Hs 2	1.664	SSE	50.77
Tp 1	8.778	Tp 2	10.780	RMSE	0.4793
fp 1	0.113	fp 2	0.092

Table 7

Fitted parameters for robust averaged DoF estimator approach using JONSWAP two peaked spectral models. The performance and the goodness of the fit process are checked by widely used statistical metrics, respectively.

Averaged DoF Spectrum				Goodness of Fit
Parameters	Values (99% Confidence Bounds)	Parameters	Values (99% Confidence Bounds)
β 1	0.760 (0.631, 0.906)	β 2	1.007 (0.831, 1.628)
γ 1	1.707 (1.430, 1.984)	γ 2	0.869 (0.684, 1.054)	R2	0.9743
σ 1	0.116 (0.077, 0.155)	σ 2	0.040 (0.028, 0.052)	Adj R2	0.9738
Hs 1	1.389 (1.060, 2.789)	Hs 2	1.280 (1.048, 1.681)	SSE	70.8
Tp 1	8.730 (8.673, 8.786)	Tp 2	10.750 (10.671, 10.840)	RMSE	0.4975

References

1. Guedes Soares, C. Stochastic Models of Load Effects for the Primary Ship Structure. Struct. Saf.; 1990; 8, pp. 353-368. [DOI: https://dx.doi.org/10.1016/0167-4730(90)90052-Q]

2. Guedes Soares, C. On the Uncertainty in Long Term Predictions of Wave Induced Loads in Ships. Mar. Struct.; 1999; 12, pp. 171-182. [DOI: https://dx.doi.org/10.1016/S0951-8339(99)00025-8]

3. Ferreira, J.A.; Guedes Soares, C. Modelling Bivariate Distributions of Significant Wave Height and Mean Wave Period. Appl. Ocean Res.; 2002; 24, pp. 31-45. [DOI: https://dx.doi.org/10.1016/S0141-1187(02)00006-8]

4. Ferreira, J.A.; Guedes Soares, C. Modelling the long-term distribution of significant wave height with the Beta and Gamma models. Ocean. Eng.; 1999; 26, pp. 713-725. [DOI: https://dx.doi.org/10.1016/S0029-8018(98)00022-5]

5. Muraleedharan, G.; Lucas, C.; Guedes Soares, C. Spectral wave energy period and peak period statistics concomitant with maximum significant wave heights. Coast. Eng.; 2023; 183, 104260. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2022.104260]

6. Rodriguez, G.R.; Clarindo, G.; Guedes Soares, C. Robust estimation and representation of climatic wave spectrum. Progress in Maritime Technology and Engineering—Martech 2018; Guedes, S.; Santos,. Taylor & Francis Group: Abingdon, UK, 2018; pp. 705-713.

7. Lucas, C.; Silva, D.; Guedes Soares, C. Climatic directional wave spectra in coastal sites. Coast. Eng.; 2023; 180, 104255. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2022.104255]

8. Donelan, M.; Pierson, W.J. The sampling variability of estimates of spectra of wind-generated gravity waves. J. Geophys. Res.; 1983; 88, pp. 4381-4392. [DOI: https://dx.doi.org/10.1029/JC088iC07p04381]

9. Medina, J.R.; Aguilar, J.; Diez, J.J. Distortions associated with random sea simulators. J. Waterw. Port Coast. Ocean Eng.; 1985; 111, pp. 603-628. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-950X(1985)111:4(603)]

10. Young, I.R. Probability distribution of spectral integrals. J. Waterw. Port Coast. Ocean Eng.; 1986; 112, pp. 338-341. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-950X(1986)112:2(338)]

11. Elgar, S. Bias of effective degrees of freedom of a spectrum. J. Waterw. Port Coast. Ocean Eng.; 1987; 113, pp. 77-82. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-950X(1987)113:1(77)]

12. Donelan, M.A.; Hamilton, J.; Hui, W.H. Directional spectra of wind-generated waves. Phil. Trans. Roy. Soc.; 1985; 315, pp. 509-562.

13. Prevosto, M.; Krogstad, H.E.; Barstow, S.; Guedes Soares, C. Observations of the High Frequency Range of the Wave Spectrum. J. Offshore Mech. Arct. Eng.; 1996; 118, pp. 89-95. [DOI: https://dx.doi.org/10.1115/1.2828828]

14. Rodriguez, G.R.; Guedes Soares, C. A Criterion for the Automatic Identification of Multimodal Sea Wave Spectra. Appl. Ocean Res.; 1999; 21, pp. 329-333. [DOI: https://dx.doi.org/10.1016/S0141-1187(99)00007-3]

15. Rodriguez, G.R.; Guedes Soares, C. Uncertainty in the Estimation of the Slope of the High Frequency Tail of Wave Spectra. Appl. Ocean Res.; 1999; 21, pp. 207-213. [DOI: https://dx.doi.org/10.1016/S0141-1187(99)00014-0]

16. Bitner-Gregersen, E.; Hagen, Ø. Uncertainties in data for the offshore environment. Struct. Saf.; 1990; 7, pp. 11-34. [DOI: https://dx.doi.org/10.1016/0167-4730(90)90010-M]

17. Veritas, D.N. DNVGL-RPC205 recommended practice: Environmental conditions and environmental loads. Tech. Rep.; 2017; 259.

18. Bitner-Gregersen, E.; Ewans, K.C.; Johnson, M.C. Some uncertainties associated with wind and waves description and their importance for engineering applications. Ocean Eng.; 2014; 86, pp. 11-25. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2014.05.002]

19. Bitner-Gregersen, E.; Waseda, T.; Parunov, J.; Yim, S.; Hirdaris, S.; Ma, N.; Guedes Soares, C. Uncertainties in long-term wave modelling. Mar. Struct.; 2022; 84, 103217. [DOI: https://dx.doi.org/10.1016/j.marstruc.2022.103217]

20. Bitner-Gregersen, E.; Magnusson, A.K. Effect of intrinsic and sampling variability on wave parameters and wave statistics. Ocean Dyn.; 2014; 64, pp. 1643-1655. [DOI: https://dx.doi.org/10.1007/s10236-014-0768-8]

21. Rodríguez, G.R.; Guedes Soares, C.; Machado, U. Uncertainty of the sea state parameters resulting from the methods of spectral estimation. Ocean Eng.; 1999; 26, 10. [DOI: https://dx.doi.org/10.1016/S0029-8018(98)00027-4]

22. Rye, H. The stability of some currently used wave parameters. Coast. Eng.; 1977; 1, pp. 17-30. [DOI: https://dx.doi.org/10.1016/0378-3839(77)90004-7]

23. Chakrabarti, S.K.; Cooley, R.P. Statistical distribution of periods and heights of ocean waves. J. Geophys. Res.; 1977; 82, pp. 1363-1368. [DOI: https://dx.doi.org/10.1029/JC082i009p01363]

24. Arhan, M.F. Analysis of the variability of spectral estimation using long continuous wave records. Proceedings of the International Conference on Sea Climatology; Paris, France, 12–23 February 1979; pp. 27-47.

25. Cavanié, A.G. Evaluation of the standard error in the estimation of mean and significant wave heights as well as mean period from records of finite length. Proceedings of the International Conference on Sea Climatology; Paris, France, 12–23 February 1979; pp. 73-85.

26. Mansard, E.P.D.; Funke, E.R. On the statistical variability of wave parameters. Technical Report, TR-HY-015, NRCC; Hydraulics Laboratory: Ottawa, ON, Canada, 1986.

27. Young, I.R. The determination of confidence limits associated with estimates of the spectral peak frequency. Ocean Eng.; 1995; 22, pp. 669-686. [DOI: https://dx.doi.org/10.1016/0029-8018(95)00002-3]

28. Peña, H.G. On the stability of the moments of the maximum entropy wind wave spectrum. Ocean Eng.; 1983; 10, pp. 79-96. [DOI: https://dx.doi.org/10.1016/0029-8018(83)90015-X]

29. Rodríguez, G.R.; Grisolía, D.; Díaz, M. On the statistical variability of some spectral bandwidth and non-linearity parameters of wind generated gravity waves. Proceedings of the 6th International Symposium on Stochastic Hydraulics; Taipei, Taiwan, 18–20 May 1992; pp. 353-360.

30. Gomes, D.; Guedes Soares, C. Estimation of directional spectrum by maximum entropy method. Proceedings of the Offshore Mechanics and Arctic Engineering; Yokohama, Japan, 13–17 April 1997; Volume II, pp. 9-17.

31. Gao, J.; Ma, X.; Zang, J.; Dong, G.; Ma, X.; Zhu, Y.; Zhou, L. Numerical investigation of harbor oscillations induced by focused transient wave groups. Coast. Eng.; 2020; 158, 103670. [DOI: https://dx.doi.org/10.1016/j.coastaleng.2020.103670]

32. Gao, J.; Hou, L.; Liu, Y.; Shi, H. Influences of Bragg reflection on harbor resonance triggered by irregular wave groups. Ocean Eng.; 2024; 305, 117941. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.117941]

33. Campos, R.M.; Islam, H.; Ferreira, T.R.S.; Guedes Soares, C. Impact of heavy biofouling on a nearshore heave-pitch-roll wave buoy performance. Appl. Ocean Res.; 2021; 107, 102500. [DOI: https://dx.doi.org/10.1016/j.apor.2020.102500]

34. Ochi, M.K. Ocean Waves the Stochastic Approach; Cambridge University Press: Cambridge, UK, 1998.

35. Goda, Y. Random Seas and Design of Maritime Structures. Advanced Series on Ocean Engineering; World Scientific Publishing Co.: London, UK, 2000; Volume 15.

36. Boccotti, P. Sea waves and quasi-determinism of rare events in random processes. Proc. Acad. Ligure Di Sci. E Lett. Ser. VIII LXXVI; 1984; 76, pp. 119-127.

37. Naess, A. On the distribution of crest-to-trough wave heights. Ocean Eng.; 1985; 12, pp. 221-234. [DOI: https://dx.doi.org/10.1016/0029-8018(85)90014-9]

38. Vinje, T. The statistical distribution of wave heights in a random seaway. Appl. Ocean Res.; 1989; 11, pp. 143-152. [DOI: https://dx.doi.org/10.1016/0141-1187(89)90024-2]

39. Boccotti, P. Wave Mechanics for Ocean Engineering; Elsevier Science: Amsterdam, The Netherlands, 2000.

40. Tayfun, M.A. Distribution of large wave heights. J. Waterw. Port Coast. Ocean Eng.; 1990; 116, pp. 686-707. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-950X(1990)116:6(686)]

41. Veltcheva, A.; Guedes Soares, C. Identification of the Components of Wave Spectra by Hilbert Huang Transform Method. Appl. Ocean Res.; 2004; 26, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.apor.2004.08.004]

42. Borgman, L.E. Confidence intervals for ocean wave spectra. Proc. 13th Conf. Coast. Eng.; 1972; 1972, pp. 237-250.

43. Priestley, M.B. Spectral Analysis and Time Series; Academic Press: New York, NY, USA, 1981.

44. Cartwright, D.E.; Longuet-Higgins, M.S. The Statistical Distribution of the Maxima of a Random Function. Proc. R. Soc. Lond. Ser. A; 1956; 237, pp. 212-232.

45. Longuet-Higgins, M.S. On the Joint Distribution of Wave Periods and Amplitudes of Sea Waves. J. Geophys. Res.; 1975; 80, pp. 2688-2694. [DOI: https://dx.doi.org/10.1029/JC080i018p02688]

46. Longuet-Higgins, M.S. On the joint distribution of wave periods and amplitudes in a random wave field. Proc. R. Soc.; 1983; 389, pp. 241-258.

47. Bartlett, M.S. Smoothing Periodograms from Time Series with Continuous Spectra. Nature; 1948; 161, pp. 686-687. [DOI: https://dx.doi.org/10.1038/161686a0]

48. Welch, P.D. The use of fast Fourier transforms for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electroacoust.; 1967; 15, pp. 70-73. [DOI: https://dx.doi.org/10.1109/TAU.1967.1161901]

49. Daniell, P.J. Discussion on the symposium on autocorrelation in time series. J. Roy. Statist. Soc.; 1946; 8, pp. 88-90.

50. Kay, S.M. Modern Spectral Estimation; Prentice Hall: Englewood Cliffs, NJ, USA, 1988.

51. Yuen, C.K.; Fraser, D.A. Digital Spectral Analysis; Pitman: London, UK, 1979; 156.

52. Pierson, W.J.; Moskowitz, L. A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. J. Geophys. Res.; 1964; 69, pp. 5185-5190. [DOI: https://dx.doi.org/10.1029/JZ069i024p05181]

53. Hasselmann, K.T.P.; Barnett, E.; Bouws, H.; Carlson, D.E.; Cartwright, K.; Enke, J.A.; Ewing, H.; Gienapp, D.E.; Hasselmann, P.; Kruseman, A. et al. Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP). Dtsch. Hydro. Z.; 1973; 12, 95.

54. Goda, Y. On the methodology of selecting design wave height. Coast. Eng. Proc.; 1981; 21, pp. 899-913.

55. Phillips, O.M. Spectral and physical properties of the equilibrium range in wind-generated waves. J. Fluid Mech.; 1985; 156, pp. 505-531. [DOI: https://dx.doi.org/10.1017/S0022112085002221]

56. Boukhanovsky, A.; Guedes Soares, C. Modelling of Multipeaked Directional Wave Spectra. Appl. Ocean Res.; 2009; 31, pp. 132-141. [DOI: https://dx.doi.org/10.1016/j.apor.2009.06.001]

57. Guedes Soares, C. Representation of Double-Peaked Sea Wave Spectra. Ocean Eng.; 1984; 11, pp. 185-207. [DOI: https://dx.doi.org/10.1016/0029-8018(84)90019-2]

58. Lucas, C.; Guedes Soares, C. On the modelling of swell spectra. Ocean Eng.; 2015; 108, pp. 749-759. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2015.08.017]

59. Titov, L.F. Wind-Driven Waves; Israel Program for scientific Translations: Jerusalem, Israel, 1969; Volume 53.

60. Thompson, E.F. Energy spectra in shallow US coastal waters. Tech. Pap.; 1980; pp. 80-82.

61. Cumming, W.E.; Bales, S.L.; Gentile, D.M. Hindcasting waves for engineering applications. Proceedings of the International Symposium on Hydrodynamics in Ocean Engineering; Trondheim, Norway, 24–28 August 1981; pp. 70-89.

62. Guedes Soares, C.; Nolasco, M.C. Spectral Modelling of Sea-States with Multiple Wave Systems. J. Offshore Mech. Arct. Eng.; 1992; 114, pp. 278-284. [DOI: https://dx.doi.org/10.1115/1.2919981]

63. Ponce de Leon, S.; Guedes Soares, C. Numerical modelling of the effects of the Gulf Stream on the wave characteristics. J. Mar. Sci. Eng.; 2021; 9, 42. [DOI: https://dx.doi.org/10.3390/jmse9010042]

64. Ponce de Leon, S.; Guedes Soares, C. Numerical study of the effect of current on waves in the Agulhas Current Retroflection. Ocean Eng.; 2022; 264, 112333. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.112333]

65. Guedes Soares, C.; Carvalho, A.N. Probability distributions of wave heigths and periods in measured sea-states from the Portuguese coast. J. Offshore Mech. Arct. Eng.; 2003; 125, pp. 198-204. [DOI: https://dx.doi.org/10.1115/1.1576816]

66. Torsethaugen, K. Model for Double Peaked Wave Spectrum; Rep. No. STF22 A96204 SINTEF Civil and Environmental Engineering: Trondheim, Norway, 1996.

67. Guedes Soares, C. On the Occurrence of Double Peaked Wave Spectra. Ocean Eng.; 1991; 18, pp. 167-171. [DOI: https://dx.doi.org/10.1016/0029-8018(91)90040-W]