Comparison of Spectrum Estimation Methods for the Accurate Evaluation of Sea State Parameters †

The monitoring of sea state conditions, either for weather forecasting or ship seakeeping analysis, requires the reliable assessment of the sea spectra encountered by the ship, either as a final result or intermediate step for the measurement of the relevant wave-motion parameters. In current analyses, different spectrum estimation methods, namely the Welch, Thomson and ARMA models, have been applied and compared based on a set of random wave signals, with different durations, representative of several sea state conditions. Subsequently, two sea spectrum reconstruction techniques were described and applied in order to detect the main sea state parameters, namely the significant wave height, the mean wave period and the spectrum peak enhancement factor. The performances of both spectral analysis and sea state reconstruction methods are discussed in order to provide some preliminary guidelines for practical application purposes. In this respect, based on current results, the Welch and Thomson methods seem to be the most promising techniques, combined with the nonlinear least-square reconstruction technique.


Introduction
Nowadays ships are equipped with a variety of sensors useful to obtain information about relevant motions and accelerations in a seaway and to increase, by means of proper active weather routing techniques, both the level of onboard comfort experienced by passengers [1] and crew safety during routine working operations. In this respect, the correct estimate of ship motion spectra is a basic issue to provide reliable information to the active weather routing decision system and increase the onboard comfort level and the safety of ship and navigation.
Before assessing the main parameters of the ship motion spectra in a seaway, the basic problems of sea spectrum reconstruction from wave time history needs to be investigated, concerning: (i) the proper selection of the most reliable technique that allows for efficiently reconstructing the wave spectrum; (ii) the minimum duration of the wave time record, required to obtain a robust estimate of the wave spectrum parameters; (iii) the amplitude of the window embodied in the Short-Time Fourier transform (STFT) when evolving sea states need to be investigated. These basic issues are essential for the optimal fitting of wave spectra to measured sea states, as widely discussed, e.g., by Mansard and Funke [2]. Therefore, a study has been undertaken in order to compare the performances of some key spectrum estimation methods, whose preliminary results were anticipated in the 2019 International having denoted by: f the wave frequency, A γ = 1 − 0.287lnγ the normalizing factor, σ the spectral width parameter, equal to 0.07 if f ≤ f p and 0.09 otherwise, f p the peak frequency, γ the peak enhancement factor and S PM the Pierson-Moskowitz spectrum: S PM ( f ) = 5 16 Hence, the PM spectrum for fully developed seas (γ = 1) is a special case of the JONSWAP one.

Outline of Some of the Main Spectrum Estimation Methods
Spectrum measurement is a key tool for the monitoring of dynamic phenomena. To perform such a measurement, signals have firstly to be acquired, through a measurement system with adequate dynamic characteristics (frequency response), and processed in order to obtain an estimate of the spectrum. Several estimation methods have been proposed since the appearance of the Fast Fourier Transform (FFT), in the late 1960s, which made their implementation as digital procedures possible. They may be parsed in two main groups-namely non parametrical and parametrical [6,7]. The seminal idea under the first group was the "periodogram", that is the square of the Fourier Transform (FT) of the series of observations, normalized in respect of the observation duration, T 0 . The periodogram, first proposed by Schuster [8] as a way for identifying hidden periodicities in a highly noisy signal, Sensors 2020, 20, 1416 3 of 15 constitutes a rough estimate of the power spectral density (PSD), since it is both biased, for a limited observation time, and it shows a high variance, which does not decrease by increasing the observation time. In fact, a longer observation time allows better spectral resolution, since the frequency spacing in the measured spectrum is ∆ f = T −1 0 , but the variance remains the same. The basic periodogram can, and should, be improved to reduce bias, which is accomplished typically by tapering or pre-whitening procedures, and to reduce variance, by averaging or smoothing. An important method that moves along these lines, is the one proposed by Welch in 1967 [9], which has been considered in the present work.
Another important method was proposed by Thomson [10], still in the area of the non-parametrical approach. The basic idea was to taper the data with different tapers to highlight different features of the signal, hence the denomination of multi-taper method (MTM).
An approach alternative to the non-parametrical is called parametrical and has been developed since the 1970s. It basically consists of considering parametrical models of the observed time series, such as the auto-regressive moving-average (ARMA) model, and obtaining estimates for the parameters involved. Once such estimates are available, the PSD of the signal may be approximated by that of the output of the model when driven by a white noise input. Methods differ in the type of model (AR, MA, or ARMA) and in the way the parameters are estimated.

Welch Method
The method consists in parsing the data record, corresponding to an (overall) observation duration T, in smaller segments of duration T 0 , with partial overlap, typically from 20% to 50%. Each segment is pre-treated by tapering with a smooth window, to reduce the bias due to spectral leakage, then the (modified) periodogram is calculated for each of them. The spectrum is obtained by averaging over such periodograms. In this way, bias is reduced by tapering and variance is reduced by averaging.
Let us then denote the series of measurements by x i = x(i∆t), where ∆t is the sampling interval, and i = 1, . . . N, with T = N∆t, and T 0 = N 0 ∆t. Let w 1 , . . . , w N 0 be a data taper, then the modified periodogram for the l-th segment is: where j is the imaginary unit. The spectral estimator is then: where n is the number of segments and m is an integer-valued shift factor, satisfying 0 < m ≤ N 0 and m(n − 1) = N − N 0 . Let us briefly discuss the application of the method and the criteria for choosing the involved parameters. Basically, two main "metrological" characteristics of the method must be considered, namely the effective bandwidth and the variance (or, equivalently, the standard deviation). It should be noted that the bandwidth of a spectral estimator is a measure of the minimum separation in frequency between approximately uncorrelated spectral estimates [7]. Therefore, the wider such bandwidth is, the worse spectral resolution is. Let us then briefly discuss the choice of the main parameters of the method. The degree of overlap is strictly related the kind of taper adopted. Basically, the smoother the taper is, the higher the degree of overlap can be adopted, which allows us to increase the number of segments and to reduce the variance. On the other hand, the smoother the window is, the larger its bandwidth is and, consequently, the worse its spectral resolution results. Welch suggested adopting a 50% overlap, and to apply a cosine (Hanning) window: we consider this a very good compromise and thus we have adopted it in the present study. Concerning the effective bandwidth, for Welch's method it can be expressed as ∆ f e = α w T −1 0 where α w is a factor that depends upon the kind of the selected taper and on the way bandwidth is defined. In the case of the Hanning window and considering a half-power bandwidth, we obtain α w = 1.44 [6]. Concerning the variance of the estimator, with a 50% overlap, a relative standard uncertainty (standard deviation) u S ( f )/S( f ) = (11/18)N 0 N −1 can be assumed, where S( f ) is the PSD and u S ( f ) is the absolute standard uncertainty [10]. These results are of great metrological import, since they allow us to keep the quality of the result under control. Once the record duration T is fixed and the kind of taper and the degree of overlap have been decided, the duration of the observation window, T 0 , remains the only design parameter to be optimised. Such optimization can be done using a trial-and-error approach, with a trade-off between the need to have a good spectral resolution, which demands for a large T 0 , and a small variance, which requires a small T 0 . It is important to note that in each trial the associated effective bandwidth and variance can be computed thanks to the above formulae, which are valid only for the Hanning window with a 50% overlap. Details of the choice made in this study will be given in Section 5.

Thomson Method
This method generalizes the tapering issue by adopting multiple orthogonal tapers. The aim is to recover the information that may be lost when using a single taper. The estimator is the average of K direct spectral estimators, each acting on the whole data record (rather than on a signal segment, as happens in Welch method) and applying a different taper. Each (partial) estimator is defined by: where h i,k is the kth data taper, usually chosen as the kth discrete prolate spheroidal sequence with parameter W, where 2W is the normalized bandwidth of the tapers, i.e., the bandwidth for ∆t = 1 s. The final estimator is thus:Ŝ where K is typically chosen to be equal to 2NW − 1. The metrological characteristics of the procedure can be kept under control by assuming an effective bandwidth ∆ f e = 2W/∆t (Hz) and considering that the estimator is approximately equal in distribution to S( f )χ 2 2K /2K, which yields a relative standard uncertainty equal to K − 1 2 . Therefore, in respect to Welch method, there is here much less arbitrariness, since, for a fixed observation time, T, the only parameter to be chosen is the half-bandwidth W, which influences both spectral resolution and relative standard uncertainty. Detail on the choices made in this study will be given in Section 5.2

Parametric Estimation Methods (ARMA)
An alternative class of methods, called parametrical, consider discrete parametrical time models of the series of observations and look for estimators of the involved parameters. Once that estimates are obtained for them, the spectrum can be calculated as a function of such parameters. For example, in the special case of (zero-mean) auto-regressive moving-average (ARMA) models, the series of measurements x i are modelled as: where ε i is a white noise process with zero mean and variance σ 2 , and a k and b k are the parameters of the model. Their estimation involves solving a non-linear optimization problem, for which some methods are available, none of which, to the authors' knowledge, prevails over the others. Once estimated,â k andb k , are available and σ 2 has also been estimated as the mean squared error (MSE) of the optimization process, the spectrum can be calculated as follows:  (8) Historically, parametrical methods were developed mainly to overcome limitations in achievable spectral resolution that are inherent in the non-parametrical approach, especially for short records. Here, the main design parameter is the choice of the model order n, which is a somewhat critical issue. In fact, the ARMA model basically refers to the response of a linear system to a broad-band, flat-spectrum excitation, which may be an acceptable assumption, e.g., in vibration measurement, in many cases. If the system under investigation may be modelled as linear, one could consider the order of the assumed model. Yet even in this favorable case, there may be problems due, for example, to minor modes in the real system not predicted by the model, of by a not flat spectrum excitation-both circumstances requiring additional degrees of freedom for being properly described, which implies a higher model order. Furthermore, when the above assumption is not justified, as in the case of the JONSWAP spectrum, the ARMA model can only be regarded as a flexible mathematical device capable of approximating a wide range of situations. Therefore, criteria have been proposed for helping model order selection, based on purely statistical/informational criteria. One the most popular in Akaike's final prediction error (FPE) criterion [6]. In this study, we have searched, for each test case, the best order both by comparing the estimated and the theoretical spectrum through visual inspection, and by considering the FPE also. Results will be discussed in Section 5.2.

Iterative Procedure
The iterative procedure proposed by Mansard and Funke [2,4] is mainly based on three steps, devoted to the assessment of the spectrum peak frequency, peak enhancement factor and significant wave height, respectively, according to the flow chart shown in Figure 1. In order to assess the peak frequency at the first step of the iterative procedure, a tentative value equal to 1.385 is assumed for the peak enhancement factor, as detailed in the following. Historically, parametrical methods were developed mainly to overcome limitations in achievable spectral resolution that are inherent in the non-parametrical approach, especially for short records. Here, the main design parameter is the choice of the model order , which is a somewhat critical issue. In fact, the ARMA model basically refers to the response of a linear system to a broadband, flat-spectrum excitation, which may be an acceptable assumption, e.g., in vibration measurement, in many cases. If the system under investigation may be modelled as linear, one could consider the order of the assumed model. Yet even in this favorable case, there may be problems due, for example, to minor modes in the real system not predicted by the model, of by a not flat spectrum excitation-both circumstances requiring additional degrees of freedom for being properly described, which implies a higher model order. Furthermore, when the above assumption is not justified, as in the case of the JONSWAP spectrum, the ARMA model can only be regarded as a flexible mathematical device capable of approximating a wide range of situations. Therefore, criteria have been proposed for helping model order selection, based on purely statistical/informational criteria. One the most popular in Akaike's final prediction error (FPE) criterion [6]. In this study, we have searched, for each test case, the best order both by comparing the estimated and the theoretical spectrum through visual inspection, and by considering the FPE also. Results will be discussed in Section 5.2.

Iterative Procedure
The iterative procedure proposed by Mansard and Funke [2,4] is mainly based on three steps, devoted to the assessment of the spectrum peak frequency, peak enhancement factor and significant wave height, respectively, according to the flow chart shown in Figure 1. In order to assess the peak frequency at the first step of the iterative procedure, a tentative value equal to 1.385 is assumed for the peak enhancement factor, as detailed in the following.

Estimation of the Peak Frequency
The assessment of the spectrum peak frequency can be performed based on different algorithms, as discussed by Mansard and Funke [2,4]. The detection of the peak period based on the largest value of the reconstructed wave spectrum could lead to a large variability of the peak frequency and should be avoided for practical purposes. In this respect, Mansard and Funke concluded that the best estimate of the peak frequency can be obtained by applying the method proposed by Read [11], and mainly based on the 5th order spectral moment of the reconstructed wave spectrum. Furthermore, they introduced in the formula provided by Read a bias corrective factor depending on the spectrum peak enhancement factor , according to the following equation: where is defined as follows: The corrective factor provided by Equation (10) can be determined once the peak enhancement factor is known, which implies that an iterative procedure is required. Hence, Mansard and Funke [2,4]

Estimation of the Peak Frequency
The assessment of the spectrum peak frequency f p can be performed based on different algorithms, as discussed by Mansard and Funke [2,4]. The detection of the peak period based on the largest value of the reconstructed wave spectrum could lead to a large variability of the peak frequency and should be avoided for practical purposes. In this respect, Mansard and Funke concluded that the best estimate of the peak frequency can be obtained by applying the method proposed by Read [11], and mainly based on the 5th order spectral moment of the reconstructed wave spectrum. Furthermore, they introduced in the formula provided by Read a bias corrective factor C f depending on the spectrum peak enhancement factor γ, according to the following equation: Sensors 2020, 20, 1416 where C f is defined as follows: The corrective factor provided by Equation (10) can be determined once the peak enhancement factor is known, which implies that an iterative procedure is required. Hence, Mansard and Funke [2,4] suggest assuming it equal to 1.02, namely γ = 1.385, at the first step of the iterative procedure to compute the peak enhancement factor γ at the subsequent step.

Estimation of the Peak Enhancement Factor
The peak enhancement factor can be determined after assessing the shape parameter k f of the bivariate Rayleigh distribution provided by Battjes and van Vledder [12]: where τ = √ m 0 /m 2 depends on the zeroth and second order spectral moments of the reconstructed wave spectrum, determined based on a lower and upper bounds of the wave frequency equal to 0.5 f p and 2.5 f p , respectively. After assessing the shape parameter k f , the tentative peak enhancement factor γ 0 is assessed by the following equation: After determining γ 0 , the peak enhancement factor is determined by the following formula provided by Mansard and Funke [2,4]: This corrective term was provided to reduce the bias in the assessment of γ, based on random generation of sea state elevations starting from theoretical JONSWAP spectra. Figure 2 provides the plot of the peak enhancement factor γ versus the shape parameter k f of the bivariate Rayleigh distribution that generally ranges from 0.4 up to 0.7 for practical purposes. Hence, after assuming in Equation (10) a tentative value of the peak enhancement factor equal to 1.385, the peak frequency f p is determined by Equation (9) and the γ value is updated by Equation (13). The procedure is iteratively performed until the relative variation of the peak enhancement factor between two subsequent steps is lower than a given threshold that can be assumed equal to 1% for practical purposes.

Estimation of the Significant Wave Height
The significant wave height H s is determined by the following formula: where m 0,c is the corrected zeroth spectral moment [2,4]: The term provided in Equation (15) allows removing the bias in the significant wave height assessment.
Sensors 2020, 20, 1416 7 of 15 generation of sea state elevations starting from theoretical JONSWAP spectra. Figure 2 provides the plot of the peak enhancement factor versus the shape parameter of the bivariate Rayleigh distribution that generally ranges from 0.4 up to 0.7 for practical purposes. Hence, after assuming in Equation (10) a tentative value of the peak enhancement factor equal to 1.385, the peak frequency is determined by Equation (9) and the value is updated by Equation (13). The procedure is iteratively performed until the relative variation of the peak enhancement factor between two subsequent steps is lower than a given threshold that can be assumed equal to 1% for practical purposes. The significant wave height is determined by the following formula:

Nonlinear Least-Square Method
The nonlinear least-square method (NLSM), proposed by Rossi et al. [3], is mainly based on the preliminary assessment of the significant wave height based on zeroth order spectral moment of the reconstructed wave spectrum, according to the following equation: The spectrum peak frequency and enhancement factor are subsequently detected by curve fitting of the theoretical JONSWAP spectrum based on the least-square method, which is based on the iterative trust-region-reflective algorithm, according to the interior-reflective Newton method [13]. Each iteration involves the approximate solution of a large linear system using the method of preconditioned conjugate gradients.

Selection of Test Cases and Random Wave Generation
In the numerical study the test cases reported in Table 1 have been considered, as they cover a wide range of sea state conditions, including both fully (γ = 1) and partly (γ > 1) developed sea states and corresponding to grades 3, 4, 5 and 6 of the Douglas (DG) Scale that allows connecting the significant wave height with roughness of sea for navigation. In Table 1 H s denotes the significant wave height, while T m (T p ) is the mean (peak) wave period. For each sea state condition, two random wave time histories have been generated based on 600 (short) and 3600 s (long) time durations, starting from the theoretical JONSWAP spectra. Hence, Figure 3 reports the input JONSWAP spectrum for the DG5 sea state condition, while Figure 4a report the short and long random wave time histories. In this respect, it must be pointed out that the wave time history was determined based on the following equation: After partitioning the theoretical input spectrum into a discrete set of components. In Equation (17)

Spectral Analysis
Spectral analysis of the time series generated as specified in the previous section was performed according to the three methods described in Section 3, namely Welch, Thomson and ARMA. Results for the DG5 sea state are provided, in graphical form, in Figures 5-7, whilst results from all test cases with reference to the recovered sea state parameters are provided, in tabular form, in Tables 2-5, in Section 5.3, where they are also discussed.
As anticipated in Section 3.2, the Welch method was applied with the classical Hanning (cosine) data window and with a 50% superposition of adjacent segments. For the long ( = 3600 s) time series, a segment duration 0 = 120 s was adopted, which corresponds to an effective bandwidth ∆ = 0.012 Hz and to a relative standard uncertainty (standard deviation) ( ) ( ) = 0.14 ⁄ . In the case of short records, a higher (worse) estimator's bandwidth was accepted, to limit the uncertainty of the estimate. Therefore 0 = 80 s was adopted, which implies an effective bandwidth ∆ = 0.018 Hz and yields a relative standard uncertainty (standard deviation) ( ) ( ) = 0.29 ⁄ . These choices were based on the previous experience gained in the preliminary study reported in Ref [3] and its validity was checked and confirmed in the present investigation. Results from this method

Spectral Analysis
Spectral analysis of the time series generated as specified in the previous section was performed according to the three methods described in Section 3, namely Welch, Thomson and ARMA. Results for the DG5 sea state are provided, in graphical form, in Figures 5-7, whilst results from all test cases with reference to the recovered sea state parameters are provided, in tabular form, in Tables 2-5, in Section 5.3, where they are also discussed.
As anticipated in Section 3.2, the Welch method was applied with the classical Hanning (cosine) data window and with a 50% superposition of adjacent segments. For the long ( = 3600 s) time series, a segment duration 0 = 120 s was adopted, which corresponds to an effective bandwidth ∆ = 0.012 Hz and to a relative standard uncertainty (standard deviation) ( ) ( ) = 0.14 ⁄ . In the case of short records, a higher (worse) estimator's bandwidth was accepted, to limit the uncertainty of the estimate. Therefore 0 = 80 s was adopted, which implies an effective bandwidth ∆ = 0.018 Hz and yields a relative standard uncertainty (standard deviation) ( ) ( ) = 0.29 ⁄ . These choices were based on the previous experience gained in the preliminary study reported in Ref [3]

Spectral Analysis
Spectral analysis of the time series generated as specified in the previous section was performed according to the three methods described in Section 3, namely Welch, Thomson and ARMA. Results for the DG5 sea state are provided, in graphical form, in Figures 5-7, whilst results from all test cases with reference to the recovered sea state parameters are provided, in tabular form, in Tables 2-5, in Section 5.3, where they are also discussed.   appearance, which is somewhat typical of this method, originally conceived for dealing with short records, they are reliable and comparable, if not even superior, to the ones provided by the Welch method, as concerns their capability of restoring the shape of the main peak of the spectrum, which we consider the most important goal, in this context. A deeper comparison will be possible in reference to the recovery of sea state parameters, to be discussed in Section 5.3. records, they are reliable and comparable, if not even superior, to the ones provided by the Welch method, as concerns their capability of restoring the shape of the main peak of the spectrum, which we consider the most important goal, in this context. A deeper comparison will be possible in reference to the recovery of sea state parameters, to be discussed in Section 5.3. As anticipated in Section 3.2, the Welch method was applied with the classical Hanning (cosine) data window and with a 50% superposition of adjacent segments. For the long (T = 3600 s) time series, a segment duration T 0 = 120 s was adopted, which corresponds to an effective bandwidth ∆ f e = 0.012 Hz and to a relative standard uncertainty (standard deviation) u S ( f )/S( f ) = 0.14. In the case of short records, a higher (worse) estimator's bandwidth was accepted, to limit the uncertainty of the estimate. Therefore T 0 = 80 s was adopted, which implies an effective bandwidth ∆ f e = 0.018 Hz and yields a relative standard uncertainty (standard deviation) u S ( f )/S( f ) = 0.29. These choices were based on the previous experience gained in the preliminary study reported in Ref [3] and its validity was checked and confirmed in the present investigation. Results from this method for sea state DG5, are presented in Figure 4. The results are quite good, especially for the long record, whilst for the shorter one the trade-off between resolution and statistical variability is more critical, as the quoted figures clearly show. In fact, this analysis is quite demanding, since in addition to a limited variance a good spectral resolution is required, to properly represent the narrow peak of the spectrum. In fact, differences between the theoretical and the estimated spectrum are predicted even in a conservative way by the relative uncertainty figures calculated for the method. Furthermore, in Figure 7c,d, results for the short signal for = 5 and for = 7, respectively, are shown. In this case, model order 6 was selected. It can be noted that for the lower order, there is an evident mismatch in the localization of the main peak, whilst with the higher order, a spurious second peak appears. In this case, for example, the FPE criterion rated order = 7 better than order = 6, which is clearly not the case.  In the application of Thomson method, a similar trade-off was made for the choice of the leading parameter for this analysis, that is the normalized half bandwidth of the taper W. In the case of the long duration signal, W = 0.0042 was assumed, which corresponds to an effective bandwidth ∆ f e = 0.017 Hz [14] and to relative standard uncertainty u S ( f )/S( f ) = 0.13 [15]. In the case of the short duration signal, W = 0.0063 was taken instead, yielding ∆ f e = 0.025 Hz and u S ( f )/S( f ) = 0.27 respectively. Results from this analysis are shown in Figure 6. In spite of their quite "noisy" appearance, which is somewhat typical of this method, originally conceived for dealing with short records, they are reliable and comparable, if not even superior, to the ones provided by the Welch method, as concerns their capability of restoring the shape of the main peak of the spectrum, which we consider the most important goal, in this context. A deeper comparison will be possible in reference to the recovery of sea state parameters, to be discussed in Section 5.3.
Concerning parametrical spectrum estimation, a preliminary investigation was reported in [3] where Autoregressive (AR) models were considered [16], with two well established estimation methods, namely Burg's algorithm and the covariance approach [6]. Results were not fully satisfactory since for low model orders (2-6) some bias in the localization of the spectrum peak was experienced, whilst when increasing the order spurious artefact peaks emerged. Therefore, in the present study ARMA models were considered instead, since they can usually provide comparable or superior approximation in respect to AR models, with lower model orders [17]. In fact, the results were better, since peak localization, with the proper model order, can now be achieved, although the resulting peak shape is still not fully satisfactory. As mentioned in Section 3.4, the search of the best model order was done both by comparing the estimated spectrum with the reference one, for various, progressively increasing orders, and by calculating the corresponding Akaike's FPE. A general pattern appeared: for low orders a poor reconstruction of the main peak resulted; this improved to some extent but proceeding further spurious peaks appeared. This behavior is similar to what happened with the AR estimator [3], but the best results here were much better than in the AR approach. In most cases order 6 led to the best results, yet in some case orders up to 10 had to be assumed. The FPE was also computed, which provided results in agreement with the visual inspection in 40% of cases. Results for this approach as applied to signal DG5 are shown in Figure 7a for a short time duration, and Figure 7b, for a long time one. In both cases model order n = 6 yielded the best results.
Furthermore, in Figure 7c,d, results for the short signal for n = 5 and for n = 7, respectively, are shown. In this case, model order 6 was selected. It can be noted that for the lower order, there is an evident mismatch in the localization of the main peak, whilst with the higher order, a spurious second peak appears. In this case, for example, the FPE criterion rated order n = 7 better than order n = 6, which is clearly not the case.

Sea State Reconstruction
Sea state reconstruction was performed based on the sea spectrum analysis carried out in Section 5.2. Particularly, the significant wave height H s , the mean wave period T m , and the peak enhancement factor γ were determined by the iterative and NLSM procedures, outlined in Section 4.1 and Section 4.2. The reconstructed parameters, corresponding to the DG3 sea state, are reported in Table 2 for both the short and long time durations, corresponding to 600 and 3600 s, respectively. Based on current results, it can be seen that: (i) The Welch and Thomson methods allow reconstructing the sea state parameters with relative errors lower than 3% for the significant wave height and mean period and 8% for the peak enhancement factor. (ii) Higher errors arise if the ARMA method is applied, especially for the assessment of the spectrum peak enhancement factor, with absolute percentage errors up to 50%. (iii) The accuracy of the methods generally increases with the time duration, even if the reliability of reconstructed parameters is sufficiently accurate for practical purposes, even if the short time history, corresponding to 600 s, is embodied in the spectral analysis. Nevertheless, this trend is not always clear, provided that a certain dependence on both the spectrum estimation method and reconstruction technique is also recognized.
Almost similar outcomes can be stressed by the analysis of the results reported in Tables 3-5, corresponding to the DG4, DG5 and DG6 sea state conditions. In fact, in all cases it is confirmed that the Welch and Thomson methods allow the efficient reconstruction of the sea state parameters, while some issues arise when using the ARMA method. Concerning the selection of the most suitable sea state reconstruction technique, no substantial differences arise between the iterative and NLSM methods. In fact, the percentage errors are comparable-even if it seems that the latter technique is slightly more accurate than the former and thus preferable, according to the results of the simulations. Starting from the results listed in Tables 2-5, Figures 8-10 report a comparative analysis between the theoretical and reconstructed JOSNWAP spectra for the DG5 sea state condition, with reference to both the short and long time durations. The reference spectrum is plotted in a continuous line, while the square and circle lines refer to the sea spectra reconstructed by the iterative and NLSM methods, respectively.
corresponding to the DG4, DG5 and DG6 sea state conditions. In fact, in all cases it is confirmed that the Welch and Thomson methods allow the efficient reconstruction of the sea state parameters, while some issues arise when using the ARMA method. Concerning the selection of the most suitable sea state reconstruction technique, no substantial differences arise between the iterative and NLSM methods. In fact, the percentage errors are comparable-even if it seems that the latter technique is slightly more accurate than the former and thus preferable, according to the results of the simulations. Starting from the results listed in Tables 2-5, Figures 8-10 report a comparative analysis between the theoretical and reconstructed JOSNWAP spectra for the DG5 sea state condition, with reference to both the short and long time durations. The reference spectrum is plotted in a continuous line, while the square and circle lines refer to the sea spectra reconstructed by the iterative and NLSM methods, respectively.
The current results confirm the main outcomes already stressed by the analysis of the results reported in Tables 2-5. In fact, it is not only confirmed that the Welch and Thomson methods allow the efficient reconstruction of the theoretical sea spectrum-it is also verified that the ARMA method fails to adequately predict the spectrum peak enhancement factor. Obviously, the reconstructed sea spectrum, corresponding to the long time duration, matches the theoretical JONSWAP spectrum slightly better, even if the results obtained based on the short time duration are not fully reliable for practical purposes.

Conclusions
The paper focused on the application of three spectral analysis techniques, combined with two sea state reconstruction methods, with the main aim of deriving the sea state parameters, namely the significant wave height, the mean wave period and the spectrum peak enhancement factor, starting from a random wave time history. In particular, two random sea surface elevations, corresponding to 600 and 3600 s, were generated, based on four theoretical JONSWAP spectra, representative of the DG3, DG4, DG5 and DG6 sea state conditions. The Welch, Thomson and ARMA methods were applied to derive the sea spectra from the wave time histories, representative of the selected sea state conditions. In this experiment, it was possible to identify optimal values for the key analysis parameters, that may constitute a useful indication for The current results confirm the main outcomes already stressed by the analysis of the results reported in Tables 2-5. In fact, it is not only confirmed that the Welch and Thomson methods allow the efficient reconstruction of the theoretical sea spectrum-it is also verified that the ARMA method fails to adequately predict the spectrum peak enhancement factor. Obviously, the reconstructed sea spectrum, corresponding to the long time duration, matches the theoretical JONSWAP spectrum slightly better, even if the results obtained based on the short time duration are not fully reliable for practical purposes.

Conclusions
The paper focused on the application of three spectral analysis techniques, combined with two sea state reconstruction methods, with the main aim of deriving the sea state parameters, namely the significant wave height, the mean wave period and the spectrum peak enhancement factor, starting from a random wave time history. In particular, two random sea surface elevations, corresponding to 600 and 3600 s, were generated, based on four theoretical JONSWAP spectra, representative of the DG3, DG4, DG5 and DG6 sea state conditions. The Welch, Thomson and ARMA methods were applied to derive the sea spectra from the wave time histories, representative of the selected sea state conditions. In this experiment, it was possible to identify optimal values for the key analysis parameters, that may constitute a useful indication for the practical application of these methods. Subsequently, the sea state parameters of the reconstructed JONSWAP spectra were determined by the iterative and NLSM methods, in order to detect the most suitable combination of spectral analysis and sea state reconstruction techniques. Based on current results, the Welch and Thomson methods seem to be the most promising techniques, combined with the NLSM method, as they allow us to efficiently reconstruct the sea state parameters, even in the short-time duration case. Contrastingly, it was verified that the ARMA method fails to adequately predict the spectrum peak enhancement factor. Current outcomes are promising for further research activities, devoted to investigating the effectiveness of the proposed techniques with reference to real sea state conditions, often characterized by double-peaked wave spectra, obtained by combining wind sea and swell waves coming from different directions. This topic will be the subject of future work.