Precise analysis of gravitational waves from binary neutron star coalescence using Hilbert–Huang transform based on Akima spline interpolation

The equation of state (EOS) information of neutron stars can be obtained by analyzing the post-merger phases of gravitational waves resulting from the coalescence of neutron star binaries. In a previous study, we proposed a method to discriminate the EOS using the Hilbert–Huang transform (HHT). The HHT comprises empirical mode decomposition (EMD) and Hilbert spectrum analysis. An essential aspect of the EMD involves the generation of envelopes through interpolating extrema values. The original EMD (used in the previous study) utilizes cubic spline interpolation. However, the cubic spline occasionally produces pseudo oscillations and overshoots that may decrease the performance of the EMD. In this study, we propose an extended version of the HHT by substituting Akima spline interpolation for the cubic spline. We compared the ability of the original HHT (based on the cubic spline) and the proposed HHT(based on the Akima spline) to discriminate the EOS. The results reveal that the proposed HHT yields a more precise analysis than the original HHT. With the proposed HHT, the number of events for discriminating the EOS is enhanced by 11.4. ..


Introduction
The Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) succeeded in recording the first direct observation of gravitational waves (GWs) in 2015 [1]. The observed GW was emitted from a binary-black-hole coalescence. In 2017, GW from binary neutron ORIGINAL UNEDITED MANUSCRIPT star (BNS) coalescence was observed by Advanced LIGO and Advanced Virgo [2]. One of the aims behind observing GWs resulting from BNS coalescence is to obtain restrictions on the equation of state (EOS) of a neutron star (NS). Several hypotheses regarding the EOS have been considered and a different model of the EOS derives from each of them. Under the assumption that a NS is cold, spherically symmetric, relativistic, and in hydro-statistical equilibrium, the EOS, mass, and radius of the NS exhibit a one-to-one correspondence among themselves [3]. Therefore, if the mass and radius of the NS can be estimated using GW data from BNS coalescence, the EOS of the NS can be constrained.
After coalescence of a BNS, a massive neutron star (MNS) may be formed (e.g., [4][5][6][7][8][9]). Certain studies have revealed that the peak frequency in the Fourier spectrum of the GW from the MNS is dependent on the radius of the fiducial NS of the BNS [6][7][8]. Accordingly, estimation of the peak frequency is essential to constrain the EOS of the NS.
Kaneyama et al. [10] attempted to obtain information on the EOS from the time evolution of the instantaneous frequency of GWs from MNSs through the Hilbert-Huang transform (HHT) [11]. In the research, two types of EOS, namely, Hyp-EOS and Shen-EOS, were considered. For the waveform of BNS based on Hyp-EOS, the instantaneous frequency of the post-merger signal owing to an MNS increases with time. However, the instantaneous frequency of the post-merger signal does not increase as for the waveform based on Shen-EOS. When the source is located at 10 Mpc, the frequency evolution of the post-merger signals on the time-frequency map can be clearly identified, and the difference in the frequency evolution between the two EOS models validated.
HHT [11] comprises two processes: a mode decomposition step called empirical mode decomposition (EMD) and a spectral analysis step called Hilbert spectrum analysis (HSA). The EMD approach decomposes the original input time-series data into certain intrinsic mode functions (IMFs); this is the preconditioning step in the subsequent step. The second part comprises the HSA: the Hilbert transform is applied to each IMF, and an analytic complex representation of each IMF is obtained. Thus, the instantaneous frequency (IF) and instantaneous amplitude (IA) of each IMF are derived. IF and IA offer a time-frequency representation of data that is well-suited for resolving non-linear and transient features of the original data. Therefore, this approach has advantages over the Fourier and wavelet transforms in analyzing non-stationary and non-linear data. Additionally, HHT is not limited by the time-frequency uncertainty principle enabling high-resolution time-frequency analysis. Considering these advantages, HHT has been utilized in various aspects of GW data analysis (e.g., [12][13][14][15][16][17][18][19][20]).
However, because HHT is a relatively new analysis method, certain improvements are required. For example, in addition to the lack of a solid mathematical description/foundation in EMD, this method is limited by several issues in terms of envelope estimation, choice of stopping criterion, mode-mixing problem, and handling of the data edges. Among these problems, envelope estimation is the most prominent challenge that requires special attention. The upper and lower envelopes play an important role in EMD. The envelope estimation has been implemented using cubic spline (CS) interpolation, and envelope estimation through CS has yielded successful results in several situations (e.g., [11,[21][22][23]). However, the CS algorithm is usually "extremely smooth" because the CS curve is a second derivative, and 2/18 ORIGINAL UNEDITED MANUSCRIPT this often results in "overshoot problems". Overshooting not only shifts the mean values of the upper and lower envelopes but also adversely affects IMF decomposition.
In this study, we extended the HHT based on Akima spline interpolation (AS) [29], which is a piecewise cubic Hermite interpolation, instead of CS. AS is based on a cubic or piecewise function that locally determines the slope of a junction point under geometric conditions. Therefore, the AS generates natural curves that appear as hand-drawn sketches. When AS is implemented in HHT, the accuracy of HHT is expected to be improved over that yielded through CS. As the primary aim, GW data from BNS coalescence was analyzed in this study to compare the results of HHT based on AS (AS-HHT) and HHT based on CS (CS-HHT) and to evaluate the performance of AS-HHT.
The remainder of this paper is organized as follows. The (CS-)HHT is briefly reviewed in Section 2. The implementation of AS on HHT is explained in Section 3.1. Section 4 presents a brief description of the numerical waveforms and simulation method employed in this study. Subsequently, the accuracy of the derived instantaneous frequency of GWs from an MNS through CS-HHT and AS-HHT are compared and discussed. In addition, we compared and discussed the differences in frequency evolution between the two EOS models. Section 5 reports a summary of the findings.

Hilbert-Huang Transform
This section briefly describes the overview of the Hilbert-Huang transform (HHT) [11].
Notably, we assume that the input s(t) is given by sampling a continuous signal in a discrete-time series: t j = j∆t for j = 0, 1, · · · , N − 1, where N is the number of data points and ∆t is the sampling interval. In this paper, subscript j is omitted for simplicity unless it is specifically needed.

Empirical Mode Decomposition
Notably, the Hilbert transform does not always yield physically meaningful results. Huang et al. [11] demonstrated that the instantaneous frequency (IF) and instantaneous amplitude (IA) yield physically meaningful results given that the time-series data satisfy the following conditions (hereinafter referred to as the IMF condition): (A) The numbers of extrema and zero crossings are either identical or differ by a count of one. (B) The mean values of the upper and lower envelopes defined using the local maxima and minima are zero at all data points.
These empirical conditions ensure that the data are symmetrical with respect to the axis origin. Because observational data s(t) do not always satisfy these conditions, decomposition via EMD is required. When using EMD, the data are implicitly assumed to be composed of many oscillatory modes of significantly different frequencies. Each mode satisfying the IMF condition is defined as an intrinsic mode function (IMF). Thus, we can decompose any data through EMD, which, in a sense, is a sifting process using a series of high-pass filters. Steps (1)-(7) mentioned hereafter and the left side of Fig. 1 present an overview of the EMD algorithm.
(1) Let the target time-series data s(t) be the i-th input of the EMD s i (t) [ Fig. 1 where i is set to 1, and s i (t) is set as the input of the following sifting process g(t). (5) Calculate the difference between g(t) and m(t) and update g(t) as the difference. (6) Check if g(t) meets the IMF conditions. In practice, absolutely satisfying the IMF conditions is a challenging task. In this regard, various approximation criteria have been proposed. In this study, we employed the stoppage criterion based on a Cauchy convergence test with a predetermined value ϵ [11]: Although this stoppage criterion appears to be mathematically rigorous, the strategy to determine the proper ϵ is not established. If this stoppage criterion is not met, then repeat steps (2) to (6) for g(t). Once the conditions are satisfied, the obtained difference g(t) is set as the i-th IMF c i (t) [ Fig. 1 (d)]. (7) Subtract c i (t) from s i (t), and feed the difference as the next input of the EMD s i+1 (t).
Increment i (i = i + 1) and repeat steps (1) to (7) to determine the next IMF until the input data become monotonic.

4/18
The remnant r(t) = s NIMF+1 (t) of subtraction of the IMFs from the input is monotonic or has only one extremum, where N IMF is the number of extracted IMFs. Consequently, the target data s(t) can be expressed as a composite of the IMFs and the remnant: The right side of Fig. 1 presents an example result of the EMD.
2.1.1. Ensemble Empirical Mode Decomposition. The decomposition result yielded via EMD is not always desirable. For example, two or more frequency modes are sometimes decomposed into the same IMF (called mode-mixing), and in other instances, a single frequency mode is split into two or more IMFs (named mode-splitting). To address this challenge, Wu and Huang [22] proposed an ensemble EMD (EEMD) method to mitigate these phenomena. EEMD is a method that utilizes white Gaussian noise to control IMF frequency bands. The EEMD algorithm is outlined in Steps (i)-(iv).
(i) Prepare N e patterns of white Gaussian noise sequences n (k) (t) (k = 1, . . . , N e ) whose mean is 0 and standard deviation is σ e . (ii) Construct an ensemble of input data by adding the noise sequences to the target data (iv) Obtain the conclusive IMFs and residue by calculating the ensemble mean as follows: In the EMD/EEMD algorithm, pre-determinant parameters such as ϵ, σ e , and N e are present. The method for setting the parameters has been discussed in the literature [23].

Hilbert Spectral Analysis
HHT can yield a time-frequency-energy expression for the given data. With this approach, non-linearity and non-stationarity are not exclusive, unlike some conventional analysis methods, including Fourier transform. For instance, non-stationarity can be treated through instantaneous amplitude (IA) and instantaneous frequency (IF) by means of Hilbert spectral analysis (HSA).
Assuming v(t) as the Hilbert transform of function u(t), it is defined as where P and * denote the Cauchy principal value of the singular integral and the convolution, respectively. If a function u(t) corresponds to the Lebesgue space L p for 1 < p < ∞, then the Hilbert transform is well defined and g(t) = u(t) + iv(t) denotes the boundary value of the holomorphic function g(z) = a(z)e iθ(z) in the upper half-plane [30,31]. Subsequently, 5/18 ORIGINAL UNEDITED MANUSCRIPT the IA A inst (t) and instantaneous phase θ inst (t) are defined as The IF F inst (t) is defined as As the IF is a continuous function, it can reflect the frequency modulation for the definition. By applying the HSA to the IMFs obtained via EMD or EEMD, the IAs and IFs of the IMFs can be derived. Through this procedure, HHT can generate instantaneous timefrequency and time-amplitude representations of the target data.

HHT based on Akima spline interpolation
As discussed in Section 2.1, HHT uses cubic spline interpolation (CS) to generate envelopes based on extrema values. However, interpolated curves rendered through CS may include pseudo oscillations and overshoots. In this study, we extended the EMD by substituting CS with Akima spline interpolation (AS). As the AS can produce more natural curves than CS, the performance of HHT is expected to improve through AS. In this paper, we denote EMD based on CS and AS as CS-EMD and AS-EMD, respectively, and HHT based on CS-EMD and AS-EMD as CS-HHT and AS-HHT, respectively.
The interpolation function of the AS is a set of piecewise polynomials of degree 3, which pass through a set of n data points (x i , y i ) for i = 1, · · · , n on the X-Y plane, similar to the cubic spline (CS). The CS is constructed such that the data points are continuous up to the second derivative of the interpolation function, that is, the CS is a C 2 differentiable function. Making the function of C 2 class even where there are jumps or where the second derivative changes markedly is likely to cause overshoot. To eliminate the overshoot, the AS does not keep the continuity of the second derivative and thus it is C 1 differentiable. The AS performs only simple geometric calculations and, unlike the CS, it does not require calculations with tolerance errors, such as the iterative solution method for second-order derivatives.
In the AS, the slope at each given data point is determined first, and then a set of piecewise polynomials is given so that the slopes are continuous. For n data points (x i , y i ), where x i ≤ x i+1 for 1 ≤ i ≤ n − 1, the slope m i between the i-th and (i + 1)-th data points is defined by and the slope t i of the curve at the i-th data point for 3 ≤ i ≤ n − 2 is given by In that case, t i is given by The slopes at the first two and last two data points are defined by From the definitions, the slope of the curve at each data point depends only on the slopes of the nearest four line segments, whereas CS uses all data points for calculating the curve. The interpolation curve is based on a cubic polynomial. The curve between the i-th and (i + 1)-th data points is expressed by The coefficients can be uniquely determined by the boundary conditions: Consequently, the forms of the coefficients are obtained as follows: We compared the results of interpolations through AS and CS with sampled data of a step function. The actual curves generated using the AS and CS are shown in Fig. 2. The points plotted with black dots represent the sampled data. Observably, the curve rendered with AS is more natural than that rendered with CS. In particular, overshoots arose only as a result of CS. Since the second derivative is also continuous in CS, overshoots likely occur at abrupt changes where the second derivative becomes large. AS suppresses the overshoots because it uses a piecewise cubic polynomial that is continuous up to the first derivative.

Performance Evaluation of AS-EMD and CS-EMD
As presented in this section, we evaluated the performances of AS-EMD and CS-EMD with a model waveform. We used a composite waveform s(t) of two types of cosine waves s 1 (t) 7 where A and B are the parameters to control the magnitude of the difference of amplitude and frequency, respectively. Because the higher-frequency component is decomposed into upper IMF, s 2 (t) should be decomposed into IMF1. Fig. 3 shows the resulted IMF1 of the AS-EMD (red line), and CS-EMD (blue line), with the waveform s 2 (t) (black line). Because the higher-frequency component is decomposed into upper IMF, s 2 (t) should be decomposed into IMF1. Fig. 3 shows the resulted IMF1 of the AS-EMD (red line), and CS-EMD (blue line), with the waveform s 2 (t) (black line). The left side of the image corresponds to the case wherein A = 0.25 and B = 0.75, and the right side corresponds to the case wherein A = 1.75 and B = 2.75. The results indicate that in the case wherein the amplitude of the higher mode is greater than the other modes, both EMDs can decompose successfully. However, in the case wherein the amplitude of the higher mode is less than the other modes, CS-EMD may fail to decompose properly. 8

ORIGINAL UNEDITED MANUSCRIPT
For a quantitative evaluation, we defined the error δ between IMF1 c 1 (t) and s 2 (t) as and the values of δ are estimated in the range 0 ≤ A ≤ 2 and 1 ≤ B ≤ 3 for both EMDs. The distributions of δ for the AS-EMD and CS-EMD are plotted in Fig. 4. The black areas imply that δ > 1 herein. We revealed that the smaller the A and B, the larger will be the δ for both EMDs, and the black area in the result of AS-EMD is smaller than that of CS-EMD. The numbers of data points where δ > 1 were 9,378 out of 40,000 for AS-EMD and 11,022 out of 40,000 for CS-EMD. These results suggest that AS-EMD performs more accurately than CS-EMD, even under severe conditions, such as conditions wherein the higher mode is relatively small. It can be explained as a bad influence of the overshoot problem of CS. Theoretically, the upper and lower envelopes generated in the first sifting step and the mean of them are identical to the lower modes 1 (t), so the higher mode s 2 (t) will remain after the subtraction of the mean from s(t). However, if overshoots exist in the envelopes, the mean is no longer identical to s 1 (t), and the result of the subtraction will have new pseudo extrema. Then, the decomposed IMF1 will not be s 2 (t) but be the mixture of two modes as shown in Fig. 3. In addition, the overshoots affect more seriously the case that the amplitude of the higher mode is small because the overshoots can easily exceed the component and change the sequence of extrema drastically. They follow that CS-EMD is not good at decomposing data whose higher mode is small, and the simulation result (Fig. 4) is consistent with that.

Gravitational Waveform from Late Inspiral, Merger, and Post-merger Phases of BNS and Simulation Data
We described gravitational waveforms based on the late inspiral, merger, and post-merger phases of BNS. Five gravitational waveforms computed by Sekiguchi et al. [5,32] were selected. These gravitational waveforms were computed for NSs using two finite-temperature equations of state (EOS). The first state equation is that of a finite-temperature EOS including the contributions of Λ hyperons (Hyp-EOS) [33]. Conversely, the other state equation represents a purely nucleonic EOS (Shen-EOS). The maximum masses of the zerotemperature spherical NSs for Hyp-EOS and Shen-EOS were approximately 1.8M ⊙ and 2.2M ⊙ , respectively. Gravitational waveforms were computed solely for symmetric binaries with M NS = 1.35 M ⊙ and 1.5 M ⊙ for Hyp-EOS (referred to as H135 and H15, respectively) and Shen-EOS (S135 and S15, respectively) and M NS = 1.6 M ⊙ for Shen-EOS (S16). The waveforms consist of the late inspiral, merger, and post-merger waveforms. We assumed that GWs from the late inspiral phase are identical regardless of the EOS models for the case that the component masses are identical, although the effect of EOS on the inspiral phase of BNS is under investigation [34]. In contrast, the waveforms of the merger and post-merger phases depend on the total mass of the binaries as well as the EOS. In the aforementioned models, a massive neutron star (MNS) is formed after coalescence, and quasi-periodic GWs are emitted from the remnant MNS. As MNSs collapse into a black hole before relaxing to a stationary spheroid, the amplitude of the quasi-periodic GWs is steeply damped at the black hole formation for H135 and H15. The frequency of the GWs  Fig. 5 Examples of the data s(t). The blue lines represent s(t) = h(t) + n(t), and the orange lines represent the injected waveforms h(t) (H135 and S15) at 5 Mpc. t merge denotes approximate merger time.
In this study, we focused on two gravitational waveforms used in previous studies: H135 and S15 [10]. The duration of the emission of quasi-periodic GWs from the MNS is τ MNS = 11 ms for H135 and > 25 ms for S15. The general relativistic numerical simulations and gravitational waveforms used in this study are described in detail in the literature [5,10,32].
Based on the sensitivity curve of Advanced LIGO (zero-detuned, high-power sensitivity curve [35]), we produced simulated Gaussian noise in the frequency domain. The sampling frequency f samp was set to f samp = 1/∆t = 16384 Hz, and the frequency range of the detector noise was set to 20 Hz to 8192 Hz. Time-series noise data, n(t), were generated using the inverse Fourier transform of the simulated noise in the frequency domain. Observational data s(t) were generated by injecting the gravitational waveform h(t) from general-relativistic numerical simulations into the Gaussian noise n(t): s(t) = h(t) + n(t), where h(t) can be expressed in terms of antenna pattern functions. For a simpler calculation, we only considered sources located optimally at the zenith of the detector.
Examples of s(t), where the H135 and S15 gravitational waveforms at a distance of 5 Mpc are injected into simulated noise, are shown in Fig. 5, and the properties of the waveforms are summerized in Table 1. The ρ opt denotes the matched-filter signal-to-noise ratio of the post-merger waveform. The definition of it of a waveform h(t) is whereh(f ) is the Fourier transform of h(t) and S n (f ) is the one-side power spectral density of noise. f low and f high are determined to contain the target phase. In Table 1, the ρ opt of the post-merger waveform in the case of an optimally oriented source at 10 Mpc are listed.

Comparison of IMF, IA, and IF characteristics
We applied both AS-HHT and CS-HHT to the simulated observational data of each EOS model. The values of the EEMD parameters in the AS-HHT and CS-HHT are summarized in Table 2. An example of the results is illustrated in Fig. 6. The injected waveform is H135 at 5 Mpc, and the simulated noise is Gaussian noise based on the design sensitivity of the Advanced 10/18 Table 1 Properties of the waveform models used in this paper. M NS denotes the mass of neutron stars (only symmetric binaries are concerned), τ MNS indicates the duration of the post-merger waveforms from the MNS, and ρ opt is the matched-filter signal-to-noise ratio of the post-merger waveform in the case of an optimally oriented source at 10 Mpc.
Model M NS /M ⊙ τ MNS /ms ρ opt H135 1.35 11 5.4 S15 1.5 > 25 6.2 Table 2 EEMD parameters used in this research. ϵ is the stoppage criterion of the sifting process, σ e and N e are the standard deviation and the number of patterns of white noise sequences added in the 1st step of EEMD. They were adjusted for each combination of the models and methods.  Figure 7 shows HHT maps, or time-frequency maps for HHT, from IMF2 to IMF5, which illustrates the time variation of IF f i (t) and IA a i (t). An HHT map looks similar to a spectrogram, or a time-frequency map, for short-time Fourier transform, wavelet transform, etc. Both of them are graphs with two geometric dimensions, the horizontal and vertical axes of which represent time and frequency, respectively. The amplitude or power of a particular frequency at a particular time is represented by the intensity or color of each point in the image. A spectrogram, in principle, fills the entire area of the 2-dimensional graph, while it may be broken up into chunks in some cases, especially for discrete, sampled data. An HHT map, on the other hand, consists of several curves representing f i (t) with color gradients representing a i (t). The HHT map enables high-resolution time-frequency analysis of a target waveform, even with strong frequency modulation.
The top panels of Fig. 7 are the HHT maps of H135, the bottom panels are the HHT maps of S15, the left panels are the results obtained with AS-HHT, and the right panels are the results of CS-HHT. The time resolution ∆t is 1/16284 s, and the frequency resolution ∆f is 20 Hz. The frequency evolution is evident, and the evolution is in good agreement with that of the injected signal.
Reportedly [5], the frequency of GWs from an MNS increases with time for Hyp-EOS, whereas the frequency is approximately constant for Shen-EOS. the IF of H135 (upper panels) increases overall until ∼ 10 ms, whereas the IFs of S15 (bottom panels) are almost constant until ∼ 15 ms. Therefore, the difference in the frequency evolution between the Hyp-EOS and Shen-EOS models can be distinguished using both AS-HHT and CS-HHT.

Comparison of accuracy of frequency evolution
This section focuses on the GWs emitted during the merger and post-merger phases. As discussed in Section 4.2, EOS can be distinguished according to the characteristics of frequency evolution of the post-merger GWs. In the next step, we evaluated AS-HHT and CS-HHT based on the accuracy with which they could estimate the GW frequencies.
To evaluate the estimation accuracy of the IF of GWs from the MNS, we defined an estimation error δ as follows [10,23], where F signal (t) denotes the IF of the injected GW signal, and WTSS stands for a weighted total sum of squares. As reported in Section 4.2, we determined that the IF is not stable when the IA is small. Therefore, Eqs. (20) and (21) are defined by considering the IA as a weight.
We compared AS-HHT and CS-HHT for cases wherein the source distance d is 5, 10, 15, and 20 Mpc. In the example plotted in Fig. 6, the GW from the merger and post-merger phase is mainly decomposed into IMF3. We calculated δ for a certain IMF into which the GW from the post-merger is decomposed in the time segment of the post-merger phase. The 13/18   target IMFs and ranges of the time segment are summarized in Table 3. For each case, we performed AS-HHT and CS-HHT on 1,000 patterns of simulated observational data, each of which is generated by adding the waveform to the simulated detector noise with a different random seed. The EEMD parameters are identical, as listed in Table 2.
The mean and standard deviations of calculated δ for each source distance case are plotted in Fig. 8. The plots for H135 clearly suggest a lack of significant differences between the results of AS-HHT and CS-HHT. The measurement accuracy δ is less than 6% for both AS-HHT and CS-HHT at 20 Mpc. In the case of S15, δ of the AS-HHT is smaller than that of the CS-HHT. These results show that the AS-HHT can estimate the frequency evolution of GW from an MNS with the same or higher accuracy than CS-HHT. In particular, this trend suggests that AS-HHT may be more suitable than CS-HHT for discriminating the EOS.

Discrimination of two EOSs
We performed tests to discriminate the correct EOS from observed GW data from an MNS by quantifying the frequency evolution with the AS-HHT and CS-HHT. Moreover, we leveraged the covariance between time and IF to quantify the frequency evolution. Assuming Cov, µ IF,3 , and µ t as the covariance, mean of IF of IMF3, and mean of time in the segment, these 14/18 parameters are defined as where t a and t b are the start and end times of the segment of post-merger phase, respectively. a and b are the corresponding element numbers of t a and t b , respectively, in the array containing time data. We compared AS-HHT and CS-HHT for cases wherein the source distance d is 5, 7.5, 10, ..., and 60 Mpc. For each case, we performed AS-HHT and CS-HHT on 1,000 patterns of simulated observational data, each of which is generated by adding the waveform to the simulated detector noise with a different random seed. From a theoretical point of view, considering H135, a positive correlation should exist between the time and IF, whereas no correlation should be present between the time and IF for S15. This implies that the obtained covariance values should distribute around a certain positive value for H135 and around 0 for S15.
The histogram of the covariance values for the case that d = 5 Mpc is exemplified in Fig. 9, and the mean µ Cov and standard deviation σ Cov of the covariance values Cov at d = 5 Mpc are also shown in Table 4. Both results obtained via AS-HHT and CS-HHT are consistent with the theoretical assumption. The standard deviation in the analysis using AS-HHT is smaller than that in the analysis performed using CS-HHT; this result clearly indicates that AS-HHT can analyze more precisely than CS-HHT. This trend is visualized in Fig. 9: the overlap of the histograms obtained via AS-HHT is smaller than that obtained using CS-HHT.
In addition, we conducted a statistical test to determine whether the distributions of the covariance differ with the EOSs or not. By performing this test while increasing the source distance, we evaluated the extent to which the source can be extended for the HHTs to discriminate the EOS. We utilized the Mann-Whitney U test (Wilcoxon rank sum test), which is a non-parametric test, because we cannot assume any distribution models for the  covariance distribution theoretically. We set the null hypothesis that the distributions of the covariance for two EOSs are identical.
The p-values of the test for the results obtained by AS-HHT and CS-HHT are listed in Table 5 and plotted in Fig. 10 as a function of the source distance. When setting the significance level as p = 0.01, the p-value does not exceed the level until the source distance reaches 45 Mpc for AS-HHT and 20 Mpc for CS-HHT. This implies that compared to CS-HHT, AS-HHT can distinguish the EOSs by analyzing observed GW data from sources that are 2.25 times farther. Assuming that the BNS event distribution is homogeneous and isotropic, AS-HHT improves the possibility of observing events for discriminating the EOSs by a factor of 11.4.

Summary
The Hilbert-Huang transform (HHT) consists of the empirical mode decomposition (EMD) and the Hilbert spectral analysis (HSA). In the EMD, the upper and lower envelopes play an important role. Moreover, the envelope estimation is typically implemented with cubic spline interpolation (CS). Although the HHT based on the cubic spline (CS-HHT) yielded various successful analyses, the CS often causes pseudo oscillations and overshoots; these aberrations may exert a negative influence on the preciseness of the EMD.
In this study, we proposed the AS-HHT, in which Akima spline interpolation was used instead of the CS. Moreover, we compared the results of applying AS-HHT and CS-HHT to model waveforms that were a composite of two cosine waves and simulated observational data of gravitational waves resulting from binary neutron star coalescence. The comparison revealed that AS-HHT was superior to CS-HHT.
In future work, we shall apply AS-HHT to waveforms originating from other sources, such as binary neutron star coalescence based on other EOS models, binary black-hole 16/18 Table 5 p-values of the Mann-Whitney U test for the distributions of the obtained covariance obtained using AS-HHT and CS-HHT. coalescence, and supernova explosions. Furthermore, we shall conduct performance tests with real detector noise.