Tuning of transfer functions for analysis of wave-ship interactions

Abstract Transfer functions are often used together with a wave spectrum for analysis of wave–ship interactions, where one application addresses the prediction of wave-induced motions or other types of global responses. This paper presents a simple and practical method which can be used to tune the transfer function of such responses to facilitate improved prediction capability. The input to the method consists of a measured response, i.e. time series sequences from a given sensor, the 2D wave spectrum characterising the seaway in which the measurements are taken, and an initial estimate of the transfer function for the response in study. The paper presents results obtained using data from an in-service container ship. The 2D wave spectra are taken from the ERA5 database, while the transfer function is computed by a simple closed-form expression. The paper shows that the application of the tuned transfer function leads to predictions which are significantly improved compared to using the transfer function without tuning.


Introduction
Transfer functions, equivalently RAOs, describing the relationship between waves and wave-induced ship responses [1] play an essential role in many applications focused on safe and energy efficient ship operation. Examples include onboard decision support systems [2][3][4], where the transfer function of a particular response can be used together with a known wave spectrum for the prediction of, say, the expected significant value ahead of time in a horizon of 30-60 min. In a similar field, transfer functions are often used in the creation of a digital twin used to simulate the behaviours of a physical counterpart [5][6][7]. Transfer functions are also applied in methods addressing real-time sea state estimation, where measurements of wave-induced responses are processed together with corresponding transfer functions to facilitate wave spectrum estimation [8][9][10].
Accurate calculation of a ship's transfer functions requires detailed information about the hull lines. In early stages of ship design projects, the hull lines are not available and, in case of existing, in-service ships, the hull lines are not necessarily available to the ship operator or third-party performance optimisation companies. Another concern of transfer functions relates to the fact that they assume small-amplitude waves and, in the strict sense, apply to one specific wave amplitude only, implying that nonlinear effects in the wave-ship interaction are not caught. This concern means that the more severe the sea state becomes, the more unreliable the transfer functions are expected to be. In turn, this conflicts with the very purpose of using transfer functions in the given case; namely for computing response levels to guide reliable operational decisions in critical situations with high waves. In simple terms, when the transfer functions are needed the most, they are the most unreliable.
Mathematically, a transfer function (TRF) is a filter that transfers waves into responses, often considered for frequency domain calculations. Symbolically, this can be written as the relation between three quantities: Response = TRF * Waves, leaving it unnoticed that a spectrum is considered on the left and right-hand sides of the equality. From a theoretical point of view, it is appreciated that, when two of the quantities are known, the third can be theoretically determined, leaving aside any practical difficulties. Thus, it is realised that if measurements of a response exist, and information about the corresponding transfer function and waves is available at the same time, a comparison of the measurement with the theoretical estimate is readily established.
Today, global wave reanalyses are made freely accessible [11,12], and through the ERA5 database it is possible to download 2D wave spectra for all parts of the oceans.

Content and novelty of study
This study presents an approach from which a vessel's transfer functions can be tuned to yield a better prediction capability when they are applied together with a given wave spectrum. In the study, it is decided to focus on a set of semi-analytically established transfer functions [13], given in closed-form and taking as input only the main particulars of the vessel. The study contains an analysis made with in-service data that has been collected from a larger container ship on a voyage crossing the Pacific Ocean.
The reason to study transfer functions that, per se, are not necessarily believed to be highly accurate is deliberate. On the other hand, it could have been a choice to study transfer functions computed by more advanced methods, e.g. strip theory or panel codes based on potential theory. In either case, the presented approach can be used to tune the available transfer functions.
In the analysis, the 2D wave spectra obtained from the ERA5 database [12] are used, but, in principle, 2D wave spectra obtained from any available source could be considered; obviously with the requirement that the source must be reliable both in terms of accuracy and in terms of availability and accessibility at the actual geographical position, in time and space, of the studied ship.
Independent on the choices related to transfer functions and 2D wave spectra, the central point of the outlined approach is to improve, i.e. tune, the available transfer functions. In turn, their use leads to a more accurate wave-response analysis, also in cases where the waves are not necessarily small-amplitude waves. It is noteworthy that the present study has some resemblance to [14,15].

Composition of paper
Section 2 presents the methodology, including the definitions and assumptions, that leads to the mathematical formulations from which the tuning coefficient can be determined. The mathematical formulations themselves are simple, but a few technical and practical concerns are noteworthy, and Section 3 discusses the most important ones. The case ship and the studied data are presented in Section 4, whereas the analysis and all associated results are given in Section 5. Finally, Section 6 draws conclusions and lists future work.

Definitions and assumptions
• The experienced seaway is the result of a short-crested wave system. Time series sequences of measured wave-induced responses are assumed to be stationary and ergodic processes. Although changing environmental and operational conditions compromise this assumption, a limitation to, say, 30 min time windows means that the assumption can be met for practical purposes. • The directional wave spectrum ( 0 , ) fully characterises the energy density of the wave system; 0 is the circular wave frequency [rad/s] and [deg] is the direction where the energy comes from, relative to the North-East-Down (NED) frame of reference. The index 0 on 0 is used to emphasise that the wave spectrum is a function of the intrinsic, and absolute, frequency. • Any global wave-induced ship response can be considered; for instance, one of the 6 DOF motion components.
• The corresponding transfer function is is the angle between the direction of the incident wave and the ship's longitudinal axis, where = 180 deg corresponds to head sea. Onwards, will be termed the relative wave heading. The transfer function depends implicitly on operational parameters such as the ship's forward speed and draught . • In general, numerical methods (closed-form, strip theory, 3D panel code, . . . ) can be used to compute an estimate ,0 ( 0 , ) of the transfer function. Linearity is assumed between waves and the induced responses. As already mentioned in the Introduction, this study calculates the transfer function by semi-analytical closed-form expressions [13].

Realistic conditions in a seaway
Because of uncertainties in operational parameters and simplifications and assumptions connected to the numerical method used for transfer function estimation, including violation of the linear assumption, an improved estimate of the transfer function is, where the tuning coefficient ≡ ( 0 , |"actual conditions"); emphasising that the tuning coefficient applies specifically to the response being considered, and that the coefficient depends explicitly on both frequency and relative wave heading. The tuning coefficient may also depend on other given parameters, such as vessel forward speed, loading condition, and possibly also significant wave height, all described by the actual operational and environmental conditions. From measurements of the response , the response spectrum can be estimated as̃( ); for instance, it can be computed by FFT or similar methods. Note, the spectrum comes as the encountered spectrum and is the encounter frequency [rad/s]. A theoretical estimatê( ) of the response spectrum can be obtained by the combined use of the wave spectrum and the transfer function of the given response in question. In this case, where is the compass heading of the ship in the NED frame of reference, and̂( , ) is given by Eq. (1). The acceleration of gravity is . The frequency conversion introduced in Eqs. (5) and (6) via Eq. (4) is necessary since the wave spectrum per se is given as a function of the intrinsic, and thus absolute, frequency 0 . For the same reason, the straight-line brackets ⟨⋯⟩ with as index is used to emphasise that evaluation happens for a given frequency of encounter ; this is discussed in more details in Section 3. The error ( ) between the theoretical estimatê( ) and the measured response spectrum̃( ) is assumed to be normally distributed. Note also that, in general, the ship's compass course and heading are different but herein they are assumed identical.
When the directional wave spectrum is known, the tuning coefficient can be estimated by minimising ( ) for any given encounter frequency . Thus, consideration of the full range of relevant encounter frequencies leads to the following nonlinear problem with least squares objective function, noting that the interval of encounter frequencies is discretised as { ,1 , ,2 , … , , } and that enters viâ( 0 , ), cf. (1). It is stressed that the integrand must be computed for given , as introduced above and further outlined in the next section.

Mapping between absolute and encounter domain
When a ship sails with non-zero forward speed relative to propagating waves, the mapping between the two frequency domains, absolute domain vs. encounter domain, is controlled by the Doppler shift as given by Eq. (4). In head sea conditions with 90 deg ≤ ≤ 270 deg, this is trivial. On the other hand, in following seas with < 90 deg and > 270 deg, and for fixed forward speed > 0 knots, a given can be the result of up to three absolute frequencies. This is shown in Fig. 1, which is a graphical illustration of the Doppler shift for the particular condition of following sea and non-zero forward speed. It is recognised that if, and only if, > 1 4 there is a 1-to −1 relationship between the encountered and the absolute frequency; otherwise there is a 1-to −3 relationship. Mathematically, the triple-valued relationship is an elementary fact that must hold in the response spectrum estimation in Eq. (2). In practical computations, this is rather delicate [16][17][18][19], but in a situation where the wave spectrum is known/given, like considered herein, the solution to the problem at hand is unique after all, in contrast to the case when the wave spectrum is unknown and is sought through an inverse problem [8,9,20,21].
The physical inference from Fig. 1 is that three distinct energy densities from the wave spectrum at 0,1 , 0,2 , and 0,3 contribute to the encounter-response spectrum at , when a ship sails in following sea if < 1 4 . In general, for any given and following sea, the theoretical response spectrum should therefore be computed by, underlining again that the straight-line brackets with as index indicate that the integrand must be computed for given . The result shows how to map frequency components of a theoretical response spectrum in the encounter domain. The central point about the mapping in Eq. (9) is to allow a comparison with a measured response spectrum, cf. Eqs. (2) and (8). On the other hand, from a theoretical point of view, it could be a choice to map the other way; that is, rather than mapping the theoretical response spectrum to the encounter domain, the measured spectrum could, in principle, be mapped to the absolute domain. In this case, however, the measured spectrum becomes heavily distorted and, from a practical point of view, the computational complications increase [16,18]. The actual implementation of the mapping from absolute domain to encounter domain of the theoretical response spectrum is sketched by the pseudo code shown in Algorithm 1, noting that the mathematical operations themselves, in relation to the spectrum estimation, are just indicated without details, such as interpolation in pre-specified discrete ranges of absolute frequencies and in relative wave headings.

Closed-form transfer functions
Semi-analytical transfer functions have been presented by [13], and closed-form expressions of wave-induced responses are given both for motion components and derived responses. The present study serves as a conceptual illustration on how to tune transfer functions, and the results herein focus on just the pitching motion. The following briefly presents how the transfer function for pitch is established by [13]. Note that, in a complete form, any transfer function is complex-valued or equivalently with two parts being the modulus and the argument. While [13] considers only the modulus of the transfer function, [22] derives also the argument for some of the responses but this is not studied herein.
The ship hull is approximated by a rectangular box with dimensions × × for length, breadth, and draught, respectively, and there is assumed to be no coupling between the motion components. The equation of motion for the pitching angle is written as, Here is the wave amplitude, = 2 0 is the wave number, assuming deep water, and is time, while differentiation with respect to time is denoted by a dot. The Doppler shift is introduced through the parameter defined by, and thus = 0 , where the Froude number = ∕ √ has been introduced. The shape effect of the hull geometry, given alone by the block coefficient , can be included most appropriately by simply taking the breadth as: where 0 is the maximum waterline breadth. Thus, the ship is represented by a homogeneously loaded box-shaped barge with the beam modified so that the total mass of the ship equals the buoyancy of the actual ship. = ∕ cos ; 10: if 90 deg ≤ ≤ 270 deg then % Head sea condition 11: else % Following sea condition 15: Fig. 1 16: ; 17: Fig. 1 20: Fig. 1 23: In the following, note that all parameters depend explicitly on the absolute frequency 0 for given vessel speed and relative wave heading , although this is not necessarily specified.
The sectional hydrodynamic damping is modelled by the dimensionless ratio between the incoming and the diffracted wave amplitudes through the following approximation, [23]: The forcing function is given by, where The Smith correction factor is approximated by Finally, the solution of Eq. (10), in terms of the modulus of the pitch transfer function, is It is realised that Eq. (18) yields zero pitching when the encounter angle between ship and waves is 90 deg, equivalently 270 deg, since a box-shaped vessel is considered. In reality, normal ships do not have fore-aft symmetry and, although little, pitching will to some degree occur while sailing in (long-crested) beam waves. In contrast to the original work [13], this has been accounted for in the present study by assuming Obviously, the many simplifications and approximations, as introduced in reaching Eq. (18), potentially lead to inaccurate results for the pitch transfer function. Nonetheless, the closed-form expression of pitch and also other motion components can be relevant for early stages of the ship design process and also other practical applications [13,[24][25][26]. As such, it makes sense to develop a method from which the result of the closed-form expression can be tuned to be in better agreement with measured data.

Optimising for the tuning coefficients
The measured response spectrum̃( ) corresponding to a time series sequence of length, say, 30 min is considered, and this constitutes a sample. In this case, the sample-specific tuning coefficient ( 0 , ) can be found by solving Eq. (8) as an inverse problem [27][28][29], considering the individual frequency-directional components as unknowns and realising that the equation system is highly underdetermined. In principle, this problem is similar to the inverse problem solved when a wave spectrum is estimated by the wave buoy analogy [8,10,20]. In the present work, however, it has already been explained that the sample-specific tuning coefficient is found from the nonlinear minimisation problem given by Eq. (8). It can be noticed that a reasonable starting guess in the optimisation is ( 0 , ) = 0, cf. Eq. (1), and note also that the size of the optimisation problem depends on the disctretisation of 0 and . That is, the number of variables to optimise is equal to the product between the number of discrete relative wave headings and the number of discrete absolute frequencies. The optimised sample-specific tuning coefficient applies specifically to the considered time series sequence; or, more precisely, to the measured response spectrum̃( ). If a set of sequential response spectrã( ){ }, = 1, 2, … , is obtained, for instance during the progression of a voyage, so that the voyage in total consists of 30-minutes periods, the mean tuning coefficient over the voyage is simply determined as, Obviously, there is little need to mention that a set of mean tuning coefficients can be determined if the complete data stream from the entire voyage, or several voyages, is decomposed into segments depending on, say, the significant wave height.
The practical implementation of the optimisation is made in MATLAB ® using the function fminunc. It is noteworthy that bounds are not imposed on the single elements of the tuning coefficient (matrix). It could be sensible to work with bounds if a more advanced method were used for computing the transfer function to be tuned; assuming such method, per se, to be in reasonable agreement with the real physics. Working with the closed-form expression, on the other hand, it is imagined that its values will not necessarily be highly accurate, in absolute terms. It is therefore decided to leave out bounds or other constraints on the tuning coefficient. In a future study, it could be of relevance to impose bounds/constraints in case unrealistic values of the tuning coefficient would be a concern. As an alternative further work, it could also be of interest, partly indicated above, to solve for the frequency-directional components directly; simply by imposing smoothness priors in a similar way as made with the Bayesian method used in the wave buoy analogy [9,30] to account for the fact that the governing equation system is highly underdetermined.

ERA5 wave spectra
The ERA5 dataset [12] provides hourly estimates of a large number of atmospheric, ocean-wave and land-surface quantities, on a discrete spatial grid spaced by 0.5 degrees in latitude and longitude (at least for ocean waves). The reanalysis combines model data with satellite observations -especially from the Copernicus EU-project [11] -into a globally complete and consistent dataset.
In this study, the global-scale hourly ERA5 directional wave spectra, together with the associated sea state parameters, have been downloaded from [12] using a routine in Python to cover the entire duration of the studied voyage, see Section 4. Each wave spectrum comprises components distributed in a frequency-direction grid, with 30 wave frequencies in the range [0.0345-0.5478] Hz and 24 circular directions spread evenly with a step of 15 • . The available spectra which were located closest to the vessel's exact GPS position were selected for each time stamps in Coordinated Universal Time (UTC). A maximum distance of approximately 25 km is observed between the ERA5 measurement point and the ship's position at the time stamp, as seen in Fig. 2. The set of spectra obtained by this nearest-neighbour interpolation was considered as ground truth for the directional wave spectra near the ship's geographic position.

Ship and route
A 7200 TEU container vessel is considered, see Table 1.
The container ship is installed with a motion sensor (XSENS, MTi-30-6A5G4), and the measurements of this study have been recorded in-service. Seven days of consecutive data were obtained while the ship followed an east-bound route across the Northern Pacific Ocean, with measurements from the Sea of Japan to the Graham Island in Canada, see map in Fig. 3, in the period of April 1st 2016 to April 7th 2016. The particular motion sensor was installed in a position off the centreline in the engine control room [31,32]. In this study, it is decided to give focus to the measured pitch motion exclusively. Hull flexibility is neglected, and a linear relationship between linear motion and angular motion is assumed.  speed-over-ground (SOG). The two measurements show rather good agreement with STW generally being the larger, indicating a projected sea current in the opposite direction of advance for the most part of the voyage. The logged mean draught amidships is also plotted in Fig. 4. The particular voyage is interesting in the sense that the STW is (nearly) constant and the variation in draught is insignificant: = 21 knots and = 14.1 m, respectively.

Encountered sea states
The integral wave parameters were computed throughout the voyage, using their mathematical definition based on the spectral moments, e.g. [33]. The results are shown in Fig. 5 for the significant wave height, the mean wave energy period, and the mean relative wave heading (RWH), emphasising that results apply to the nearest ERA5 grid point. The consequence of being off the native grid point(s) has been investigated by [34] in a study focused on integral parameters, and it is concluded that it can be influential but generally the effect will be small. Anyhow, this study considers the actual 2D wave spectra in the analysis, and results are thus constrained to rely on conditions at the nearest grid point. As such, it is left to note that small uncertainties will be caused as a result of a mismatch in ship position and the nearest native grid point, but such uncertainties are not considered herein.
The next section contains a more detailed analysis of the pitch measurements collected during the voyage. At this point, however, it is interesting to note that the sea state is relatively low in the beginning of the voyage, the first day or so, with a small significant wave height and a low mean wave period. The low period means that waves with a short wave length, equivalently high frequency, characterise the wave systems. In theory, this means in turn that no substantial motions are induced, since the ship acts as a low-pass filter. This is easily confirmed by inspection of the transfer function for pitch, as calculated by the closed-form expression, cf. Eq. (18), without tuning, see Fig. 6.

Results and discussions
The following section contains the analysis of the study, and thereby the section presents and discusses the results. As a practical comment it should be noted that the study considers pitch only. It is therefore unnecessary to specify that it is the pitch response spectrum which is under investigation. Instead, pitch spectrum and response spectrum are used interchangeably. Besides, for convenience, any time series sequence and its associated data is described by a sample index rather than a date-time stamp. Thus, it follows that the data of the seven days, equivalent to 168 h, is characterised by 168 samples.

About discretisation
As indicated already, calculations follow in a discretised format, where the transfer function of pitch is computed for discrete ranges of absolute frequencies and relative wave headings. Specifically, the results apply to discretisations: * 0 = (0.01 ∶ 0.01 ∶ 0.30) Hz and * = (0 ∶ 10 ∶ 350) deg, where the frequency is given as cycles per second 0 = 0 2 . In this case, = 36 relative wave headings and = 30 frequencies are considered. It is realised that this particular discretisation is inconsistent with that used for the ERA5 2D spectra, cf. Section 3.4, and calls for interpolation of the latter.
For a given time series sequence, the measured response spectrum is obtained from the estimated autocovariance function with frequency smoothing using a Parzen window function [35]. The resulting spectrum is interpolated to values at * = (0.01 ∶ 0.002 ∶ 0.40) Hz.    It is noteworthy that for the particular discretisation of absolute frequencies and relative wave headings, the tuning coefficient corresponding to one single sample, i.e. one time series sequence, will come as a matrix with × elements. Note also that when referring to a tuning coefficient corresponding to a single sample, singular form of coefficient is used, despite the coefficient in fact contains multiple (matrix) elements.

Sample-specific tuning
Figs. 7 and 8 show the effect of tuning, when theoretical estimates of response spectra are compared with corresponding measured spectra for some arbitrarily selected cases taken from samples 1-120. Two preliminary remarks: (1) The tuned response spectra U.D. Nielsen et al. have been computed using the sample-specific tuning coefficients, directly resulting from the minimisation in Eq. (8). The results corresponding to using a mean tuning coefficient, cf. Eq. (21), for response prediction is studied later in Section 5.3. (2) Throughout all samples of response spectra, including those not shown, the measured spectrum shows a low-frequency peak located at a frequency of encounter in the range 0.02-0.025 Hz, and the value of the peak itself is about 0.25-0.3 deg 2 s. This observation is also discussed in more details later.
The 0-th order spectral moment 0 corresponding to a spectrum is calculated by The 0-th order spectral moment, equivalent to the variance of the response time series, has been computed for the measured spectra and for the theoretical estimates with and without tuning, respectively. The outcome is shown in Fig. 9, but emphasising that sample-specific tuning is made only together with samples 1-120. In the plot, each point corresponds to a given 30-minutes time series sequence. As mentioned in Section 3.4, the ERA5 2D wave spectra are updated on an hourly basis. Therefore, it is just every second response spectrum, equivalently time series sequence, which is considered; noticing that the 30-minutes sequence is selected from a window with 15 min prior to the hourly ERA5 update, and 15 min succeeding the ERA5 update. Fig. 9 shows that the use of a sample-specific tuned transfer function leads to nearly exact matches between the measured and theoretically estimated results, confirming the good agreement observed from the arbitrarily selected samples of response spectra shown in Figs. 7 and 8. A more direct comparison between the measured and theoretically estimated response spectra is illustrated in Fig. 10 presenting the normalised summed-squared error (SSE): The SSE should be considered as a relative measure, since the single number in itself has little meaning. On the other hand, it is evident from Fig. 10 that tuning of the transfer function leads to improved response spectrum estimates. The observed trend of decreasing SSE, especially for the tuned results, is a coincidence since, at this stage, it is not attempted to optimise (''train'') depending on the amount of available data.

Response prediction based on mean tuning coefficient
It was introduced by Eq. (21) that a mean tuning coefficient can be computed from a number of sample-specific tuning coefficients. In fact, the equation suggests to use samples. For the particular data stream one mean tuning coefficient can for example be computed from all of the first 120 samples, realising that, in this case, the mean tuning coefficient is obtained from about 70% of all the available data and with no special account to sea state whatsoever. Obviously, this kind of division of the data is not the most sensible, since the tuning coefficient, as mentioned, is expected to depend on, notably, significant wave height. The reason to not work with a more meaningful division is simply due to the limited amount of data included in the analysis. Notwithstanding, conceptually this situation can be used to illustrate the potential for making response predictions with a tuned transfer function, but it is clear that a future study, based on a (much) larger dataset, must carefully examine the mean tuning coefficient's dependency on ''external'' parameters ( , , , …). Fig. 11(a) shows the frequency-variation of the mean tuning coefficient for different relative wave headings, emphasising that absolute frequencies are considered on the horizontal axis. It can be inferred that for certain combinations of frequency and relative U.D. Nielsen et al.   wave heading the transfer function undergoes considerable tuning, and it is also noted that the mean tuning coefficient takes negative values only for a few combinations of frequency and relative wave heading. It is observed that for frequencies around 0.05 Hz, the transfer function will be the most amplified in stern-quartering seas ( ≃ 40-60 deg), while at a higher frequency band around 0.10 Hz the situation is somewhat changed, since the most amplification of the transfer function happens for relative wave headings corresponding to beam seas to head seas ( ≃ 90-180 deg). Fig. 11(b) presents the frequency-variation of the standard deviation of the sample-specific tuning coefficients obtained over the 120 samples. The specific plot reveals that the sample-specific tuning coefficients are not constant over the considered samples; that is, the sample-specific tuning coefficients do not take the same value for given combinations of frequency and wave headings. This is likely a result of the changing sea state during the progression of the voyage, see Fig. 5. For the actual application of the mean tuning coefficient towards response prediction it should be noted that the smaller the standard deviation, the better the mean tuning coefficient is expected to be. This will be discussed below. For illustrative purposes, polar versions of the plots in Fig. 11 are shown in Fig. 12, where it is noted that the (absolute) frequency along any radial line goes from 0 Hz at the circle's centre to 0.3 Hz at the perimeter.
The capability of the tuned transfer function can be assessed by studying the predictions obtained on the samples 121-168, emphasising that these samples have not themselves been used in the calculation of the mean tuning coefficient. Fig. 13 shows an arbitrary selection of response spectra from samples 121-168, and each plot presents three spectra: The theoretical estimate based on a transfer function that has not been tuned, the measured spectrum, and the theoretical estimate based on a tuned transfer function, where it is the mean tuning coefficient which is applied. It can be seen that the theoretical estimates generally match better with the measured spectra if the transfer function is tuned. However, it is also evident that the agreement generally is not as good as when the sample-specific tuning coefficients are applied, cf. Figs. 7 and 8. A quantification of the comparison is again established by comparing the 0-th order spectral moment, and the outcome is shown in Fig. 14. Overall, good agreement exist.
Returning to the details reflected by the response spectra, cf. Fig. 13, it is clear that a distinct and intensive low-frequency peak at 0.03-0.04 Hz (encounter frequency) almost consistently occurs in the tuned theoretical estimates. This peak is also appearing in all of the measured spectra, albeit less pronounced. In fact, in the measured spectra, two low-frequency peaks are appearing; the one at 0.03-0.04 Hz and the other peak at 0.02-0.025 Hz of which the latter was introduced already in Section 5.2. If focus initially is on the peak at a frequency of encounter = 0.03-0.04 Hz, this peak is simply a direct result of the present wave spectrum. This can be seen from the plots in Fig. 15 presenting the wave spectra that match the response spectra in Fig. 13. As can be seen, the wave system has its main peak located at a frequency ∼ 0.1 Hz and with waves coming from a southwesterly direction. Converted to measurements on the ship, by considering the ship heading and forward speed through the Doppler shift in Eq. (4), this means that the ship encounters the (most significant) waves with a frequency = 0.03 − 0.04 Hz, noting that in the given cases the relative wave heading corresponds to following to stern-quartering waves ( = 0 − 60 deg). Obviously, the given wave systems will also contribute to response energy densities at lower encounter frequencies than 0.03 Hz, cf. Eq. (4), but from the wave spectra there are no indications that an actual second peak in the pitch spectra should be appearing at = 0.02 − 0.025 Hz. It is therefore speculated that the pitch sensor could have a malfunction resulting in low-frequency noise around frequencies 0.02-0.025 Hz and with a somewhat constant density about 0.25-0.3 deg 2 s. Unfortunately, it is difficult to draw any definitive conclusions in this direction. This is substantiated by returning to Figs. 7 and 8, showing the examples of response spectra from samples 1-120. The matching wave spectra have been included in the appendix in Fig. A.18, and it can be seen that, by and large, almost all of the wave spectra show wave energy of both smaller and larger density in the spectral region located in the southwesterly corner, obviously noting that the magnitude of the energy density clearly depends on the case in study. The observed low-frequency peak in the pitch spectra at a frequency of encounter = 0.02 − 0.025 Hz could therefore also be the result of a somewhat persistent swell system in the first part of the voyage, represented by the wave spectra in Fig. A.18(a). In the later part, as the voyage progresses, the swell system fades away, but, by coincidence, the present wave system at the ship position changes its direction, in a counter-clockwise direction, cf. Fig. A.18(b).
As already partly indicated, only further detailed investigations will be able to conclude about a possible malfunction in the pitch sensor, but this is left to future work. Instead, it is here important to mention the consequence of all of the discussions above. Thus, the implication in relation to the tuning coefficient is that, when computing the mean tuning coefficient from samples 1-120, the optimisation leads to a need for amplifying (i.e. tuning) especially the transfer functions for relative wave headings corresponding to following to stern-quartering seas. This is a finding naturally resulting from the optimisation of the tuning coefficient, since the optimisation simply ensures the best possible fit to the measured spectrum, in all cases of samples 1-120, no matter the presence  of any noise, if any. In turn, when the mean tuning coefficient is applied for pitch motion prediction in samples 121-168, the transfer function (matrix) elements, specifically corresponding to following to stern-quartering seas, become too significant; simply because during the samples 121-168 the ship is in fact sailing in following to stern-quartering waves, as reflected by the mean relative wave heading, cf. Figs. 5 and A.18. Altogether, this results in theoretical tuned estimates of the response spectra where the lower-frequency part tends to be overestimated. It will be beyond the scope of this study but it could be of relevance to apply a high-pass filter to the measurements in order to study the effect of the low-frequency ''noise''. At this stage, suffice it to say that the important overall finding from the previous analysis is that the introduction of a tuning coefficient, in terms of Eq. (1), leads to improved response predictions.
Studying Fig. 14, it is remarkable how the drop in the agreement between the tuned predictions and the measurements of 0 in samples 158-168 follows the quite sudden change in sea state, cf. Fig. 5, starting just before sample 158, equivalent to noon on April  7th. It is believed that this observation is not necessarily directly connected to a possible malfunctioning sensor. Instead it could be the result of overfitting, understood in a similar way like discussed above, i.e. that the sample-specific tuning coefficients are overfitted to some particular conditions, with the consequence that the tuned theoretical prediction gives poor results for different conditions. Alternatively, it could be speculated that the observation results because of inconsistencies in the ERA5 wave spectra around samples 156-168, realising that the wave period drops steeply while the wave height, in contrast, increases. On the other hand, these changes are also followed by a sudden change, albeit with a slight delay, in the (relative) wave direction, which could indicate a change from, say, a swell-dominated sea to a distinct wind-dominated sea, as partly confirmed by the actual 2D wave spectra in Fig. 15. On top of this, to add a level of complexity, it is noticed that the vessel forward speed suddenly is increased, followed by an immediate reduction, around the same period of time (noon, 7th April), cf. Fig. 4. Altogether, it is beyond doubt that only very detailed investigations of the exact conditions can possibly explain exactly what is correct. While this is not part of the present paper's scope, an additional investigation has been made involving the transfer function calculated by strip theory.

Validation with strip theory and general comments
For validation purposes, the pitch transfer function has also been computed by strip theory using an in-house code based on [36]. The result is included in Fig. 16 together with the closed-form transfer function without tuning. It is appreciated that the agreement is somewhat reasonable, although the closed-form expression tends to start filtering away waves at a lower frequency than the strip theory method, especially for wave encounter angles from 0 deg (following sea) to 90-100 deg (beam sea). The immediate consequence of this observation is twofold: (1) If a tuning coefficient were to be calculated for the strip theory results, the tuning coefficient would presumably be somewhat similar to what has been found in the preceding for the closed-form expression. (2) The fairly good agreement could be an(other) indication that the measurements are corrupted due to a possible sensor malfunction. As a supplementary comment to (1), it must be realised that, despite the ''overall'' good agreement, the deviations between the transfer functions (strip theory vs. closed-form) at the higher frequencies could lead to changes in the resulting tuning coefficient, if calculated based on strip theory. This is so, since in the calculation of the predicted response spectra, there are contributions from all wave encounter angles (0-360 deg), cf. Eq. (2). In fact, the importance of the higher-frequency contributions can be seen from Fig. 17 that shows the computed response variance, i.e. 0 , predicted using both sets of transfer functions, in either case without tuning. It is noticeable that the strip theory results are different to the closed-form results only for parts of the data. Thus, the largest deviations are observed for samples (app.) 84-168, where the strip theory results consistently match better with the measurements. As indicated previously and with reference to Fig. 5, those samples reflect operating conditions where the ship encounters, in sequential order, beam sea over following sea to beam sea, and, as realised from Fig. 16, this set of wave encounter angles is also where the largest deviations between the transfer functions occur; emphasising that the deviations are not significant in the low-frequency region. Returning to Fig. 17, it is appreciated that samples 1-84, on the other hand, represent cases with good agreement between the two sets of predicted results, but noticing that the predictions themselves are off compared to the measurements. Overall, the main point to note from the past discussion is that, unfortunately, the inclusion of -or validation with -strip theory based transfer functions has not shed so much additional light on the issue about the possibility of a malfunctioning sensor. At present, it is therefore repeated to note that further future analyses are required to make any solid conclusions in this direction. In the event of a future study, it is important to also consider the uncertainty associated with the ERA5 2D spectra; noticing that the ERA5 spectra do not represent the actual ground truth, although assumed herein. Besides, the analysis relies on the use of the wave spectrum at the nearest grid point, meaning that the mere distance from this position to the exact ship position, cf. Fig. 2, will be influential and adds uncertainty to the analysis.
As a final point of discussion, related to the (dis)agreement between the (tuned) predictions and the measurements of 0 , it is of relevance to point out the main assumption of the current analysis, which is that of linearity between wave and response, expressed by Eq. (2). Despite the inclusion of a tuning coefficient, to partly correct the linear result, still the analysis assumes that the total response can be obtained simply by adding together an infinite number of component regular responses [1], each with fixed amplitude, frequency, direction, and phase. However, for truly multidirectional seas, as observed with some of the 2D polar plots in Figs. 15 and A.18, it could be that the neglect of nonlinear effects resulting because of both directional and frequency wave-wave interactions is simply inappropriate/too approximative. In addition, the analysis uses a sum of just 30 wave components, frequency-wise, on the interval (0.01-0.30 Hz), cf. Section 5.1. Above all, the varying level of discrepancies along the voyage between the measured pitch variance, on one hand, and the pitch variance computed from not-tuned as well as tuned transfer functions, on the other hand, reveals that there can be multiple independent possible causes for a need to correct/tune a transfer function. In this study, the causes include: the varying accuracy of the interpolation of ERA5 wave spectra at the ship's exact position; the inaccuracy of linear theory in multidirectional seas; an early drop of the closed-form transfer function magnitude (compared to strip theory) at high frequencies in following to beam sea.

Conclusions and future work
The results and discussions in the preceding have shown that it is possible to tune the transfer function of pitch, so that theoretical predictions based on the tuned transfer function agree with a reasonable accuracy to corresponding measurements. In particular, it has been observed that the theoretical predictions will be significantly improved compared to predictions based on the original transfer function without inclusion of a tuning coefficient.
The presented study has many limitations, among which the most important ones are: (1) The study has been focused on one response only, and, for illustrative purposes, it was decided to work with a semi-analytical transfer function. (2) The studied inservice data is from one voyage only, and optimisation of the mean tuning coefficient was not made with account to sea state. (3) The sample-specific tuning coefficient was obtained through optimisation considering a nonlinear minimisation problem, and it was not investigated if a true global minimum was reached. (4) The ERA5 2D wave spectra are taken as the ground truth. While being at positions exactly coinciding with the native grid points this can be a reasonable assumption, but clearly, for an arbitrary position off the grid, it will be associated with uncertainty to use the wave spectrum ''just'' at the nearest grid point.
The limitations inherently calls for future work, and it is therefore suggested to perform additional studies addressing the following questions and comments: • If additional responses (heave, roll, accelerations, . . . ) are considered, will it then be advantageous to compute one overall (mean) tuning coefficient? • What results are obtained if the tuning is associated with transfer function(s) computed by more advanced numerical software?
• When semi-analytical closed-form expressions are considered for transfer functions, it could be of interest to directly optimise some of the parameters entering the mathematical formulations. • Rather than computing the mean tuning coefficient from a somewhat arbitrary set of samples, like in this study due to scarcity in data, it should be attempted to compute the mean tuning coefficient as a function of sea state, notably as function of significant wave height. This could be a way to better take account of any nonlinear behaviour not considered when transfer functions are used in spectral calculations. In addition, results from the study also indicated that it could be sensible to consider the mean tuning coefficient's sensitivity to relative wave heading. • If the sample-specific tuning coefficient is obtained from an optimisation, it will be important to ensure that the global optimum is found. As an alternative to optimisation, it can be a possible solution to directly solve for each of the matrix elements in the tuning coefficient. This approach can follow ideas from the wave buoy analogy [20]. • In line with the previous bullet, it could be of relevance to impose bounds on the tuning coefficient, if solved from a minimisation problem like in this study. This could counteract unrealistic values of a resulting tuned transfer function, obviously with the potential cost of a reduced agreement between measurements and tuned predictions. • All uncertainty associated with any inconsistency between the measured response spectrum and the corresponding theoretical estimate has been ascribed to the transfer function. Part of the uncertainty will be associated with the wave spectrum [10,34]. This must be further investigated.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.