A time–frequency analysis of gravitational wave signals with non-harmonic analysis

...................................................................................................................Weapplynon-harmonicanalysis(NHA)toanalyzegravitationalwavesignalsfromcompact binarycoalescences.AfterinvestigatingthegeneralpropertiesofNHAbyusingseveraltest signalswithamplitudeandfrequencymodulation,weanalyzegravitationalwavesignalsfrom theinspiralingofcompactbinarycoalescences.Weﬁndthatthetime–frequencypropertyof thesignalcanbereproducedaccurately.Weanalyzethesignalsfrominspiralingbinaryblack holesinjectedintosimulatedGaussiannoise.Weﬁndthatwecanidentifythesignalevenifthe signal-to-noiseratiointhesenseofamatchedﬁlterisaslowas10.Weinvestigatetheaccuracy oftheinstantaneousamplitudeandfrequencyofthesignalrecoveredwithNHA.Weanalyze gravitationalwaveeventsdetectedbyLIGOandVirgo,andﬁndthatgravitationalwavesignals withstrongamplitudelikeGW150914canbeshownclearlyonthetime–frequencymapby NHA.Wealsoﬁndthateveninthecaseofalow-amplitudesignallikeGW170817,thesignal canbeclearlyseenbychoosinganappropriatewindowsizethatdependsonthefrequencyband. .


Introduction
After considerable experimental effort to construct high-sensitivity gravitational wave detectors, gravitational wave (GW) signals were finally detected in 2015 during the first observing run of the advanced LIGO detectors [1]. Advanced Virgo [2,3] started observations in 2017. At the time of writing, 11 GW signals have been observed [4]. GW170817 was a GW from the coalescence of binary neutron stars, and others are GWs produced by the coalescence of binary black holes (BBHs).
Within a few years, KAGRA [5][6][7] will also be operational. There will also be a LIGO-India detector in the future. The international network of laser interferometric GW detectors composed of these detectors will be a very powerful instrument to observe astronomical phenomena related to compact objects such as black holes (BHs) and neutron stars (NSs). From the event rate of the coalescence of BBHs, if the sensitivity of the LIGO and Virgo detectors are improved as predicted PTEP 2019, 063F01 K. Yanagisawa et al. in Ref. [8] then it is expected that several tens of BBH events will be detected in the next few years, along with several binary neutron star events. The coalescence of NS-BH binaries may also be detected in the near future. There are other possible sources, including core collapse supernovae, rotating neutron stars, giant flares from magnetars, etc. [9].
Several methods can be used to analyze gravitational wave signals. The most widely used method when a waveform is known theoretically is matched filtering, in which data from detectors are crosscorrelated with theoretical waveform (templates). However, if such templates are not available, we cannot use matched filtering. One method to analyze data in such a case is the excess power statistics, in which the excess of power of the data on the time-frequency plane is searched by assuming the time duration and frequency band of signals. Since this is a kind of time-frequency analysis, the choice of the basis used to represent the data in time-frequency space is important. Although the short-time Fourier transform (STFT) was used initially, in the current pipeline of LIGO [10] the Wilson-Daubechies-Meyer transform is used [11]. This method was used to detect LIGO's first event, GW150914 [12].
There is another time-frequency method that is used to analyze gravitational wave signals: the Hilbert-Huang transform (HHT) [13][14][15][16][17][18][19][20]. It consists of an empirical mode decomposition followed by Hilbert spectral analysis to evaluate the instantaneous frequency and amplitude of each mode.
In discrete Fourier transformation (DFT) analysis, the frequency resolution is limited to f = 1/T (T is the length of the data). Thus, the presence of spectral components whose frequencies are not integer multiples of f results in spectral leakage, and thus hampers the accuracy of the frequency evaluation. Even though the choice of the window function may allow some side lobe control, spectral leakage is unavoidable in this method. On the other hand, NHA can estimate spectral components with non-integer multiples of f . Thus, the spectral leakage is weaker, which is beneficial when we treat data with a short time length. These make it possible to simultaneously achieve good frequency and time resolutions. This suggests that NHA can be useful when we perform time-frequency analysis as an alternative method to those based on STFT.
In this paper, we discuss the possibilities of applying NHA to gravitational wave data. In particular, we discuss the NHA algorithm first introduced by Yoshizawa et al. [22], which we hereafter refer to as "the NHA method." An initial attempt to analyze gravitational wave signals in noisy data from a gravitational wave detector was performed in Ref. [31]. We extend the analysis and investigate the accuracy of evaluating the instantaneous frequency with NHA. Since NHA is a kind of timefrequency analysis, it does not assume particular signal waveforms before the analysis. Thus, the most interesting application of NHA to GW data analysis may be the analysis of signals whose templates are not available. However, as a first step we investigate the application of NHA to GWs from inspiraling BBHs whose waveforms are approximately given analytically. One issue we encounter when we apply NHA to the analysis of GWs is that the noise level of GW detectors is much higher than the previous industrial applications. Thus, it is important to investigate whether we can recover the time-frequency evolution of GW signals accurately in noisy conditions. In particular, since there is no mathematical proof on the uniqueness of the time-frequency representation with NHA, it is important to investigate the performance of the NHA method in GW data analysis. By using GWs from inspiraling BBHs whose frequency evolution is approximately given analytically, it becomes 2/25 PTEP 2019, 063F01 K. Yanagisawa et al. easy to evaluate the accuracy of the instantaneous frequency derived with the NHA method. This is the reason for using BBH signals in this paper.
Finally, we apply NHA to real gravitational wave data from LIGO and Virgo. In this analysis, in order to optimize the time-frequency representation of real data we introduce a method which uses multiple window sizes depending on the frequency band [32]. With this method, we find that we can obtain good time-frequency representation of real gravitational wave signals detected by LIGO and Virgo.
This paper is organized as follows. We discuss the NHA method in Sect. 2. In Sect. 3, we compare NHA with short-time Fourier transforms and Hilbert spectral analysis by using simple test signals. In Sect. 4, we analyze GWs from inspiraling BBHs, and compare the results with other methods. We also consider GWs in Gaussian noise, and evaluate the accuracy of recovering the frequency evolution. In Sect. 5, GW signals detected by LIGO and Virgo are analyzed and the time-frequency maps are shown. Section 6 gives a summary and discussion.

Non-harmonic analysis
In this section we present an overview of the NHA method. We consider a time series data set, equally sampled with an interval t, written as {x 0 , x 1 , . . . , x M −1 }, where M is the number of data points. From this data set, we pick a data point, x c , and a short data segment with N data points around In order to realize this, x c must be located within x N /2 ≤ x c ≤ x M −N /2 . The time length of this data segment is T = N t, which we call the window size of NHA. In order to simplify the notation, we rewrite the index of this time series as In NHA, we decompose the data with a linear combination of cosine waves: where A k , f k , and φ k are the amplitude, frequency, and phase, respectively. Contrary to the discrete Fourier transformation (DFT), the frequencies f k are not limited to integer multiples of the frequency interval, f = 1/T . Instead, any frequency is allowed for each cosine wave. We call the expression in Eq. (1) the NHA decomposition. The full details of the computing algorithm of NHA are described in Ref. [22]. Here, we explain it briefly. We assume that the amplitudes of Eq. (1) are ordered such that A 1 > A 2 > · · · > A k max . We determine the modes with order ( x j by minimizing the fitting error, defined as Since there may be several local minima in this function, multiple steps are used to derive the best-fit value of (A, f , φ). First, the discrete Fourier transform of the data, X , is computed using the fast Fourier transform (FFT). The initial value of A is obtained by the maximum of |X |. Initial guesses for f and φ are computed from X m , where m is the frequency index which realizes the maximum of |X k |. These initial values are expressed as (A 1,0 , f 1,0 , φ 1,0 ). Now, we use an iteration. Suppose (A 1,m , f 1,m , φ 1,m ) are the values obtained by the mth iteration. Then the second stage of NHA is described as follows: (i) By fixing A 1,m , and by adopting (f 1,m , φ 1,m ) as the initial value, the steepest descent method is used to find the value of (f 1,m+1 , φ 1,m+1 ) until F(A, f , φ) converges to a certain level. (ii) By fixing (f 1,m+1 , φ 1,m+1 ), the steepest descent method is used to find A (1,m+1) that realizes the minimum of F(A, f , φ). (iii) Steps (i) and (ii) are repeated until the value of F(A, f , φ) converges to a certain level.
At the third stage, we use Newton's method as follows. Suppose again that (A 1,m , f 1,m , φ 1,m ) are the values obtained by the mth iteration.
(iv) By fixing A 1,m , and by adopting (f 1,m , φ 1,m ) as the initial values, Newton's method is used to find the value of (f 1,m+1 , φ 1,m+1 ) until F(A, f , φ) converges to a certain level. (v) By fixing (f 1,m+1 , φ 1,m+1 ), the steepest descent method is used to find A (1,m+1) that realizes the minimum of F(A, f , φ). (vi) Steps (iv) and (v) are repeated until the value of F(A, f , φ) converges to a certain level. When The third stage is necessary to obtain (A 1 , f 1 , φ 1 ) accurately. The final step is subtraction: (vii) Once we obtain the best-fit value, (A 1 , f 1 , φ 1 ), we subtract it from the data and define the residual data: We repeat this process and derive (A k , f k , φ k ) until the residual data {x (k+1) j } becomes sufficiently small. If we are interested only in a few spectrum components with strong amplitude, this can be stopped after a desired iteration time.
In order to make a time-frequency map, we pick all data points x j within N /2 ≤ j ≤ M − N /2 in Eq. (1). At each sampling point, t = j t, we plot points for each frequency component f k . The color of each point is determined by the amplitude, {A k }, of the frequency component. We use both grayscale and color.
The above procedure has been applied for various data, and it is known that it works empirically. However, it is difficult to prove mathematically that we can always obtain the global minimum of F(A, f , φ). It is usual for there to be multiple local minima in F(A, f , φ), which makes it difficult to determine the global minimum. In the case of gravitational wave data, this effect is expected to be large since the amplitude of the signals is not very large compared to that of the noise. The GW signal is strongly affected by the noise, which may create a lot of local minima. Thus, it is important to check whether NHA can be used for the time-frequency representation of GW signals. In the subsequent sections we examine how accurately we can extract the instantaneous frequency of BBH signals with NHA.

Performance of NHA for test signals
In this section we first apply NHA to a model signal: s(t) = cos φ(t). The sampling frequency is set to 10000 Hz and the duration of signal is set to 1 s. Thus, the number of data points is N sample = 10000. The window size for NHA is set to 25.6 ms (256 points).  For the frequency evolution, we consider the two cases shown in Figs. 1 and 2: Signal 1: For comparison, we compute the time-frequency maps by using STFT and compare the results with those derived with NHA. In STFT, the Hamming window is used as a window function. The window size is set to 1024 points (∼0.1 s), and the time shift length between adjacent segments is one sampling interval (10 −4 s).
As a reference, we also show the time-frequency maps computed with Hilbert spectral analysis (HSA). This is based on the analytic representation of a signal s(t) [33] defined as where where PV denotes the Cauchy principal value. The instantaneous amplitude IA(t) and the instantaneous phase θ(t) are given as and The instantaneous frequency IF(t) is defined as the derivative of θ(t): Note, however, that it is well known that it is possible to obtain a reliable instantaneous frequency with HSA only in the case when the data satisfy very special conditions. In particular, it is usually difficult to use HSA for noisy data. 1 Thus, the results with HSA are shown only for noiseless cases. Figure 1 shows the results of the analysis of Signal 1. Figure 1(c) shows the results from the STFT analysis: the linear frequency evolution is well represented, but the spectral lines are broadened; both time and frequency resolutions are therefore low. Figure 1(d) shows the results obtained from the analysis of Signal 1 with HSA. Except for some distortion at the beginning and at the end of the analysis interval, the frequency evolution is well represented. The results obtained with NHA are shown in Fig. 1(e). We can see that NHA can accurately capture the frequency evolution with much smaller side lobes (spectral dispersion and distortion) than the other two methods. captured with moderate accuracy, several spikes appear along the instantaneous frequency trajectory. Figure 2(e) shows the NHA result. We find that although there are small distortions around 20 Hz, the frequency evolution is represented clearly.

Inspiraling compact binares
In this section we consider gravitational waves from inspiraling BBHs. We use gravitational waveforms calculated up to the third post-Newtonian order beyond the lowest order [34]. We assume that each black hole is non-rotating. We also assume that the orbit is circular.
The gravitational waveform in the time domain is given as where and M = m 1 + m 2 (m 1 , m 2 are the masses of the black holes), η = m 1 m 2 /M 2 , t c and c are the time and phase at coalescence, and r, G, and c are the distance to the source, Newton's gravitational constant, and the speed of light, respectively. Here, the formula for the phase, orb (t), contains only up to second post-Newtonian order. The full formula up to the third post-Newtonian order is very long and is not shown here-it can be found in Eq. (317) of Ref. [34]. In the amplitude of this waveform, only the quadrupole mode is included. In this paper, we only consider the contribution of the lowest post-Newtonian order effects to the amplitude. Thus, in the amplitude, we use the formula for the orbital angular frequency, which is given at the lowest order as From Eqs. (14) and (17), we find that the amplitude evolves as ∝ (t c − t) −1/4 . Since the frequency of gravitational waves, f , is twice the orbital frequency, it is given as f = orb /π. Thus, from Eq. (17), the instantaneous frequency of gravitational waves evolves approximately as We notice that the amplitude and the frequency diverge at t = t c . This suggests the breakdown of the post-Newtonian approximation. The waveform is therefore terminated before t = t c when the orbital angular frequency reaches that of the innermost stable circular orbit (ISCO), which is roughly given as We quantify the strength of signals with the matched filter signal-to-noise ratio (SNR), ρ c , defined as where S n (f ) is the one-sided power spectrum density of the noise.

Analysis of the inspiraling binary black hole waveform
Here we analyze the signal from a coalescing binary black hole without noise with three analysis methods: STFT, HSA, and NHA. We use the signal with masses m 1 = m 2 = 15 M . The signal starts about 2 s before the coalescing time t c . The sampling frequency is set to f s = 16384 Hz. At the start time, the frequency of the signal is 23 Hz. The ISCO frequency is 146.6 Hz. Since there is no noise, the amplitude of the signal is arbitrary. In the STFT analysis, a data length of 62.5 ms is used to perform the STFT, and the Hamming window is used as a window function. The STFT is performed at each time step to obtain the time-frequency map. In the NHA analysis, a window size of 31.25 ms is used for each time step. Figure 3 shows the results. We find in Fig. 3(b) that, in the STFT analysis, although the frequency evolution can be seen in the time-frequency plane, it is difficult to find the detailed spectral structure because the signal is spread broadly. A better result can be obtained with HSA, as shown in Fig. 3(c). Although some fluctuations can be seen at the beginning and end of the signal, the spectral characteristics of the signal can be adequately determined with HSA. NHA also works for this signal, as shown in Fig. 3(d). We find that the spectral structure can be easily determined with high frequency and time resolutions, although some side lobes in frequency components appear around the end of the signal. From Fig. 3, HSA seems to be the best of the three methods for the noiseless case. However, as discussed in the previous section, it is difficult to use HSA for noisy data containing multiple frequency components. Thus, hereafter, we do not discuss HSA.
In order to investigate the effect of a different window size in NHA, Fig. 4 shows the results with window sizes t w = 31.25 ms, 62.5 ms, and 125 ms. We find that the result with 62.5 ms is similar to that with 31.25 ms. However, at the end of the signal, the very fast increase of the frequency creates some side lobes on the time-frequency plane in the 62.5 ms case. This feature is even more pronounced in the 125 ms case. These features can also be seen in Figs. 5 and 6, in which we show the difference between the peak frequency and amplitude evaluated by NHA and the theoretical value as a function of time. Near ISCO, the deviation of the value with NHA from the theoretical one starts earlier for the 125 ms case. In NHA, the obtained frequency and amplitude are average values within each window. When the window size is larger, since the frequency evolution in a window becomes larger, it becomes difficult to obtain accurate values for the frequency and amplitude. In that case, it also becomes difficult to represent the signal with fewer cosine waves. The residual of the fitting produces multiple side lobes, which is manifested in the 62.5 ms and 125 ms cases. We find that the frequency intervals of the side lobes are approximately 1/t w . Since the NHA itself does not have any feature related to this frequency interval, this should be related to the FFT used for the initial guess of the NHA algorithm. The frequencies of the side lobes are the local minima of Eq. (2) near the frequency bin of the FFT. On the other hand, for the 31.25 ms case, since it is possible to obtain a good fit with a single cosine wave, the residual is very small. Thus, significant side lobes are not produced.  These results suggest that it is necessary to keep the window size small when we analyze signals whose frequency evolves rapidly. Note, however, that the differences for both frequency and amplitude are very small until just before ISCO. The difference in the frequency stays within a few Hz, except near ISCO.

Analysis of the inspiraling binary black hole waveform in Gaussian noise
In this section we apply the NHA method to gravitational wave signals from coalescing binary black holes embedded in the simulated noise of a gravitational wave detector. The noise data is produced by assuming the KAGRA detector [5]. We use a design sensitivity of KAGRA, called VRSE(B), 11  whose theoretical noise power spectrum density is publicly available at KAGRA's website [35]. By using this, we first generate random Gaussian noise in the frequency domain, and then apply an inverse FFT to obtain time series data. The data we analyze are obtained by adding the BBH signals, as discussed in Sect. 4.1. We set the minimum frequency to 10 Hz when we compute the matched filter SNR in Eq. (19).
In order to reduce the computation time of NHA, we apply a low-pass filter (LPF) [36] to the data as a pre-processing step. The LPF is implemented with an infinite impulse response filter to minimize the computation required. In order to reduce the ripple effect in the passband, a type II Chebyshev filter is used which has fairly steep cutoff characteristics, and no ripple in the passband [37]. The LPF cutoff frequency is set to 300 Hz.
We first consider the case when the matched filter SNR of the BBH signal is 100. Figure 7 shows the results obtained with NHA with the different window sizes of 31.25 ms, 62.5 ms, and 125 ms. The spectral line in the vicinity of 150 Hz is due to noise present in the noise curve of KAGRA. We find that, at this SNR, we can clearly identify the chirp signal for all window sizes. At earlier times (t < 1.5 s) when the frequency is low, oscillation of the frequency is observed for the 31.25 ms window which does not exist in the noiseless case in Fig. 4(b). This oscillation is reduced as the   Figure 8 shows the results of NHA analysis for SNR = 10, 20, 30, and 50 with window sizes t w = 31.25, 62.5, and 125 ms. Figures 8(a), (e), and (i) show the cases for SNR = 10. We find that it is possible to identify the signal even for the case of SNR = 10. When t w = 31.25 ms, although the frequency resolution is not so good, the final frequency evolution around ∼100 Hz can be seen vaguely. For t w = 62.5 ms, the frequency resolution is clearly better than that for t w = 31.25 ms for t < 1.5 s. However, for t > 1.5 s (Fig. 8(e)), there are several lines that appear due to noise. The recovered signal also seems to be distorted for several times. It is thus difficult to identify the signal in this case. Similar features can be observed for t w = 125 ms (Fig. 8(i)). The frequency resolution is good for t < 1.5 s, and many of the noise components which appear in Figs. 8(a) and (e) disappear in Fig. 8(i). However, in the final part of the signal the frequency evolution is not PTEP 2019, 063F01 K. Yanagisawa et al.   Fig. 7. In order to clearly see the effect of noise, the density scale is represented as 20 log 10 A for amplitude A in this figure. traced properly, and some side lobes also appear. Oscillation of the frequency at t < 1.5 s appears in all cases; the oscillation is strongest for t w = 31.25 ms and weakest for t w = 125 ms. This is because the fluctuation of the instantaneous frequency due to noise is captured in the case of the small window size. If the window size is larger, the evaluated frequency is averaged in the window, and the fluctuation becomes smaller.
As the SNR of the signal increases, it becomes easier to identify the signal. However, the abovementioned features remain the same: when the window size is smaller the oscillation of the frequency becomes strong at t < 1.5 s, but the final rapid increase of the frequency can be captured accurately. In the case of t w = 125 ms, the frequency fluctuation at t < 1.5 s becomes weaker as the SNR becomes larger. In the cases of t w = 31.25 and 62.5 ms, the rapid increase of the frequency near t ≈ 2 s can be traced better as the SNR becomes larger.
To quantify the statistical fluctuation of the instantaneous frequency determined by NHA, we compute the distribution of the difference between the signal-plus-noise case and the signal-only case of the frequency determined by NHA. We identify the signal for the signal-plus-noise cases as the frequency at which the maximum amplitude is realized at each time step. Figure 9 shows the evolution of this frequency for the data in Figs. 7 and 8. We find that the fluctuation of the frequency PTEP 2019, 063F01 K. Yanagisawa et al. of the signal is largest for a window size of 31.25 ms. The fluctuation of the frequency becomes weaker when the window size becomes longer.
In Fig. 10 we plot the histogram of the difference of the frequency of signals between the signalplus-noise case and the signal-only case. Figure 11 shows the standard deviation of Fig. 10 as a function of SNR for each window size. We find, as expected, that when the SNR is smaller, the standard deviation is larger. From SNR = 100 to 10, the standard deviation is about 2 Hz ∼ 6.5 Hz for t w = 31.25 ms, about 1 Hz ∼ 4 Hz for t w = 62.5 ms, and about 0.5 Hz ∼ 3 Hz for t w = 125 ms.
We again conclude that the capability to follow the frequency evolution is improved when we use a short window size for NHA. However, when we use a short window size, the influence of the noise becomes larger. When we use a long window size, although the time resolution becomes worse and it becomes difficult to follow the rapid time evolution of the frequency, the influence of the noise becomes smaller. Therefore, to accurately visualize the frequency evolution of gravitational wave signals, an appropriate window size needs to be used in the NHA analysis. 16 Fig. 10. Histogram of the difference between the frequency evaluated by NHA (Fig. 9) and the analytical frequency.

NHA time-frequency map of LIGO/Virgo events
To demonstrate the NHA method, in this section we visualize the gravitational wave signals detected by LIGO and Virgo on the time-frequency map using NHA method. We use the data available at the LIGO Open Science Center [38,39]; the data files used are listed in Table 1. The data are whitened, band-limited to 20-500 Hz by using an ideal bandpass filter, and down-sampled from the original 4096 Hz to 1024 Hz.
As discussed in the previous sections, the choice of the NHA window length is important to capture rapid frequency evolution and to avoid the influence of the noise. We thus introduce multiple window lengths into the analysis [32]. We use a longer window length for low-frequency bands, and a shorter 17   window length for high-frequency bands. We divide the data into several frequency bands, as shown in Table 2. We first set the maximum frequency to be 500 Hz. The lower frequency border of the band is decided by multiplying by 2/3 (half octave), so we have 500 × 2/3 = 333 Hz. Each frequency band is determined in this way. When we analyze the frequency band with maximum frequency f m the window length is determined as α/f m , where α is an adjustable parameter which is determined empirically. We compute 30 components of the NHA decomposition for each frequency band. Thus, in total, 390 components are plotted in each map.
The time-frequency maps for LIGO/Virgo events are shown in Figs. 12-14. For these cases we used α = 1.75. In each figure, the time is set to 0 at the coalescence time of each data.
We find that GW150914 can be visualized most clearly on the NHA time-frequency map. This is because the performance of NHA is sensitive to the instantaneous amplitude signal-to-noise ratio, and the matched filter signal-to-noise ratio is not necessarily a good indicator. GW150914 has the largest PTEP 2019, 063F01 K. Yanagisawa et al. amplitude relative to the noise among the events analyzed in this paper. The matched filter SNR for GW150914 is about 19 (Livingston) and 13 (Hanford), while the duration is only 200 ms from 30 Hz to coalescence. For a given matched filter SNR, a signal of shorter duration is advantageous for time-frequency analysis. For comparison, we show the spectrogram with STFT in figure 15. GW170104 can also be visualized clearly, particularly for Livingston data. Other high-mass BBH signals, LVT151012 and GW170608, are not very clear but we can roughly identify the signal on the time-frequency map. However, it is difficult to identify GW151226. The matched filter SNR is only about 10 (Livingston) and 8 (Hanford), and the duration from 35 Hz is about 1 s. The low matched filter SNR and longer duration make it difficult to identify signals on the time-frequency map. Figure 16 shows the results for GW170817. Since this signal is low-mass binary neutron star coalescence, it continues for about 60 s from 30 Hz, which is much longer than the binary black hole cases above. This means that although the matched filter SNR is not so small, the amplitude of this signal is very small. Thus, it is difficult to find the signal in Fig. 16. We try to optimize the window length so that we can identify the signal. Figure 17 shows the results using the longer window size given in Table 3, which corresponds to α = 25. We find that by using a longer window size we can clearly identify the signal for the LIGO-Livingston and LIGO-Hanford cases. On the other hand, it is still difficult to find the signal for the Virgo case because of the very low amplitude of the signal for this case.

Summary and discussion
We have discussed the possibility of using the NHA method as a time-frequency analysis tool to analyze gravitational wave signals. First we analyzed simple waveforms whose frequency evolves with time and compared the results derived by STFT, HSA, and NHA, finding that the frequency evolution can be evaluated clearly with NHA, while the other two methods have some difficulties in evaluating the frequency evolution. We also compared three methods by using gravitational waves from 15 M inspiraling binary black holes lasting for 2 s. We found that NHA can represent the frequency evolution clearly. In this case, HSA also works well.  In order to investigate different choices for the NHA window length, we compared the results of the analysis of a BBH signal with different window lengths. We found that while the window length is short, it is possible to track the frequency evolution accurately. However, when window length is longer, many side lobes of spectral lines appear when the frequency evolution becomes very rapid just before the coalescence time. These spectral lines do not exist in the original signal.
Next, we added simulated Gaussian noise to the BBH signal, and compared various SNRs and NHA window lengths. We found that in order to represent the rapid frequency evolution before the coalescence, it is useful to use a shorter window size. However, if the window size is shorter, the frequency fluctuation induced by the noise becomes larger even for the low-frequency region where the frequency evolution is not very rapid. If we use a longer window size, the frequency fluctuation becomes smaller. As expected, the effect of noise is larger when the SNR is smaller. These facts show 20  the importance of the choice of NHA window size to accurately represent the frequency evolution while minimizing the effect due to noise.
We quantitatively evaluated the accuracy of the instantaneous frequency of the signal, finding that in the case of SNR = 10, we can evaluate the instantaneous frequency of the BBH signal with an accuracy of 3-7 Hz, depending on the window size. When SNR = 100, the accuracy becomes 0.5-2 Hz. These results were derived for a BBH signal with m 1 = m 2 = 15 M . These results are very promising for the usefulness of NHA in evaluating the instantaneous frequency of gravitational wave signals.
Based on these investigations, we analyzed the data of LIGO/Virgo events by introducing multiple NHA window sizes. This means that for the same data segment we divided the frequency band into small frequency bands, and used a different window size for each frequency band. We used a shorter window size for the high-frequency region, and a longer window size for the low-frequency region. The choice of window size has so far been empirically determined. With this method, we find that we can obtain very a clear time-frequency map of signals with large amplitude, like GW150914. Even for an event like GW170817 whose duration is very long and whose amplitude is small, we can obtain a clear time-frequency representation of the signal by choosing the window size appropriately.
Since the NHA assumes that the data consist of many sinusoidal waves for short data segments, NHA seems to be useful for analyzing signals which oscillate many cycles and whose frequency evolution is not so rapid. Thus, NHA might be useful for the analysis of post-bounce gravitational wave signals of core collapse supernovae [40]. It may also be useful to analyze the various line noises that exist in the interferometer data. In fact, in the KAGRA detector NHA was already used 21   as a tool to monitor the line noise during the first test operation of the interferometer with a simple configuration [7].
One issue with NHA is the computation time. Since we need to minimize the function in Eq. (2) with respect to three parameters, it takes much more time than FFT. To analyze 32 s of LIGO/Virgo    Table 3. In these figures, the signal is clearly visualized for the Hanford and Livingston data. The figure for the Virgo data is not shown since it is still difficult to identify the signal in the Virgo data. data with the window lengths in Table 2, it took about 5 min using a modern CPU. However, if we use the window lengths in Table 3, the computation time is about 144 hr. On the other hand, to produce a time-frequency map with STFT, Fig. 15, it takes only about 0.48 s. The NHA computation time could be reduced by using a graphics processing unit; this will be discussed in the future. 23