We present the application of a novel method of time-series analysis, the
Hilbert-Huang Transform, to the search for gravitational waves. This algorithm
is adaptive and does not impose a basis set on the data, and thus the
time-frequency decomposition it provides is not limited by time-frequency
uncertainty spreading. Because of its high time-frequency resolution it has
important applications to both signal detection and instrumental
characterization. Applications to the data analysis of the ground and space
based gravitational wave detectors, LIGO and LISA, are described

Camp, Jordan B.

Cannizzo, John K.

Numata, Kenji

English

arXiv

We present the application of a novel method of time-series analysis, the Hilbert-Huang Transform, to the search for gravitational waves. This algorithm is adaptive and does not impose a basis set on the data, and thus the time-frequency decomposition it provides is not limited by time-frequency uncertainty spreading. Because of its high time-frequency resolution it has important applications to both signal detection and instrumental characterization. Applications to the data analysis of the ground and space based gravitational wave detectors, LIGO and LISA, are described

NASA Technical Reports Server

Source of Acquisiliao 
NASA Cioddard Space Flight Center 
Application of the Hilbert-Huang Transform to the Search for Gravitational PJaves 
Jordan B. Camp,* John K. Cannizzo,t and Kenji ~ u m a t a ~  
Laboratory for Gravitational Physics, Goddard Space Flight Center, Greenbelt, Maryland 20771 
We present the application of a novel method of time-series analysis, the Hilbert-Huang Transform, 
to the search for gravitational waves. This algorithm is adaptive and does not impose a basis set on 
the data, and thus the time-frequency decomposition it provides is not limited by time-frequency 
uncertainty spreading. Because of its high time-frequency resolution it has important applications 
to both signal detection and instrumental characterization. Applications to the data analysis of the 
ground and space based gravitational wave detectors, LIGO and LISA, are described. 
PACS numbers: 04.80.Nn, 07.05.Kf, 95.55.Ym 
Introduction - A vital area of gravitational wave 
(GW) research is the development of data analysis al- 
gorithms that are able to  detect non-linear, transient 
signals. The data analysis effort is especially relevant 
as the Laser Interferometer Gravitational Wave Obser- 
vatory (LIGO) [I, 21 undertakes a one-year science run 
at  full design sensitivity. Also, the Laser Interferometer 
Space Antenna (LISA) [3] is engaged in the development 
of analysis algorithms to search for its candidate sources. 
The Hilbert Huang Transform (HHT) [4] presents a 
fundamentally new approach to the analysis of time se- 
ries data. Its essential feature is the use of an adap- 
tive time-frequency decomposition that does not impose 
a fixed basis set on the data, and therefore it is not lim- 
ited by the time-frequency uncertainty relation charac- 
teristic of Fourier or Wavelet analysis. This leads to a 
highly efficient tool for the investigation of transient and 
non-linear features. Applications of the HHT include ma- 
terials damage detection [5] and biomedical nlonitoring 
16, 71. Because General Relativity is an inherently non- 
iinear theory, and because LIGO, LISA, and other GW 
detectors produce a great variety of non-linear and tran- 
sient signals, the HHT has the promise of being a power- 
ful new tool in the search for gravitational waves. This 
article describes the application of the HHT to GW data 
analysis. 
The HHT proceeds in two steps. First, the process of 
Empirical Mode Decomposition (EMD) reduces the time- 
series under analysis into components, known as Intrinsic 
Mode Functions (IMFs) , thereby "sifting" or separating 
out the different frequency scales of the data. The IMFs 
when summed reproduce the original time series. The 
sifting is done adaptively, with no a priori structure im- 
posed on the data. 
The IMFs have a vertically symmetric and narrowband 
form that allow the second step of the HHT to be ap- 
plied: the Hilbert transform of each IMF. As explained 
'Electronic address: campQlheapop. gsf c .nasa. gov 
t ~ l s o  at Physics Department, University of Maryland, Baltimore 
County, Baltimore, Maryland 21250 
t ~ l s o  at  Department of Astronomy, University of Maryland, Col- 
lege Park, Maryland 20742 
below, the Hilbert transform obtains the best fit of a si- 
nusoid to each IMF a t  every point in time, identifying an 
instantaneous frequency (IF), along with its associated 
instantaneous amplitude (IA). The I F  and IA provide a 
time-frequency decomposition of the data that is highly 
effective a t  resolving non-linear and transient features 
Instantaneous Frequency - The IF is generally ob- 
tained from the phase of a complex signal z ( t )  which is 
constructed by analytical continuation of the real signal 
x(t) onto the complex plane 171. By definition, the ana- 
lytic signal is 
where y(t) is given by the Hilbert Transform: 
(Here P denotes the principal Cauchy value.) The am- 
plitude and phase of the analytic signal are defined in the 
usual manner: a(t)  = Iz(t) 1 and Q(t) = arg z ( l )  .The ana- 
lytic signal represents the time-series as a slowly varyicg 
amplitude envelope modulating a faster varying phase 
function 181. The I F  is then given by w ( t )  = dQ( t ) / d - l ,  
while the IA is a(t). We emphasize that the IF, a function 
of time, has a very different meaning from the Fourier 
frequency, which is constant across the data record being 
transformed. Indeed, as the IF is a continuous f~nction, 
it may express a modulation of a base frequency over a 
small fraction of the base wave-cycle. Conditions on the 
form of x(t)  that allow the application of this procedure 
are: 1) vertical symmetry with respect to the local zero 
mean, and 2) the same number of zero crossings and ex- 
trema [4]. 
A time-frequency analysis of a frequency modulated 
signal serves to illustrate the concept of instantaneous 
frequency. An important GMT source for LISA is the 
inspiral of supermassive black holes (SMBH). The GW 
waveform from this source has a base frequency which 
chirps as the masses approach, and is also phase and 
amplitude modulated due to spin-orbit coupling IS]. The 
time evolution of the base frequency of the GW from the 
https://ntrs.nasa.gov/search.jsp?R=20070034819 2019-08-30T01:50:28+00:00Z
inspiral is given to first order by [lo]: 
where p = m1rn2/(ml +m2) is the system reduced mass, 
iWt is the system total mass, G is Newton's constant, c 
is the speed of light, and to is the coalescence time. 
Fig.1-(a) shows the waveform from an SMBH binary 
inspiral with masses ml and ma of 1 0 h n d  lo5 so- 
lar mass, and spins sl and $2 of 0.95m: and 0.95mz. 
This waveform is used as input to the HHT and FFT- 
based spectrogram to perform the time-frequency anal- 
ysis, shown in Figs.1-(b), -(c), -(d). It is clear that the 
extraction of the I F  and IA provides high time-frequency 
resolution of the waveform, showing detail of the chirp as 
predicted by Eq.(3) (dotted line in Fig.1-(b)), and also 
amplitude modulation due to the spin-orbit coupling. In 
contrast, the decomposition of the chirp through the FFT 
is seen to be limited by the time-frequency uncertainty re- 
lation: htA f - 1, where At and A f are the time interval 
and (Fourier) frequency bandwidth. The decoupling of 
the modulations of the I F  and IA provided by the HHT 
allows a detailed look a t  the amplitude and frequency 
evolution of the waveform, and will aid significantly in 
determining the SMBH binary parameters. 
Empirical Mode Decomposition - The waveform of the 
previous example has an approximate symmetry about 
zero that satisfies the conditions listed above and allows 
meaningful application of the Hilbert ~ransform for ex- 
~raction of the IF. To express an arbitrary time-series in 
this form the process of EMD is applied [4]. 
The EMD process consists of forming an envelope 
about the extrema of the data by cubic spline fitting, 
and then subtracting the average of the envelope from 
-the data. The extrema of remainder are then fitted, and 
the process is repeated as many times as necessary to 
obtain a waveform that is symmetric about zero mean 
(within a predetermined tolerance). Once this has oc- 
curred, the waveform is labeled IMFl and is subtracted 
from the original time series, removing the highest fre- 
quency content, and allowing the remainder to be sifted 
again to obtain the next IMF. This procedure is illus- 
trated in Fig.2. 
- The sifting identifies and removes components of the 
data first at the highest frequencies, then down in fre- 
quency to the lowest trends. For white noise the EMD 
acts as a dyadic filter, so that the central frequency as- 
sociated with lMFn is of order fs,,,/2n where f,,,, is 
the data sampling frequency [ll]. A consequence of the 
sifting is .that the IMFs are vertically symmetric about 
zero mean (although both the frequency and amplitude 
of an IMF vary in time.) The narrowband and symmet- 
ric form of the IMFs are essential to allow the Hilbert 
Tra~isform to be applied in a meaningful way [4]. 
EMD an,d the Hilbert Transform: The HHT- To illus- 
crate the application of the full HHT to GW analysis, we 
look at the identification of a signal in white noise (Fig.3). 
Time [sec] 
FIG. 1: Comparison of IF and FFT-based spectrogram of 
GW signal from SMBH inspiral with spin-orbit coupling. (a) 
GW waveform. (b),(c) IF and IA show sharp det+l of fre- 
quency and amplitude modulation. The dotted line on (b) 
plots Eq.(3). (d) FFT spectrogram is blurred because of time- 
frequency spreading and spurious harmonics. 
As an example relevant to LIGO analysis, where strong 
signals are not expected, we consider here a weak signal 
relative to the noise. An issue for the HHT is that low 
signal to noise tends to degrade the IF and IA identifica- 
tion, as the signal extrema becomes distorted. A rough 
measure of this tendency was determined from simula- 
tions by injecting sine waves of varied signal to noise, 
frequency and duration into white noise. This showed 
Af /  f - 2 J=/SN, where A f is the extracted 
frequency width of a periodic wave of frequency f and 
duration 6t  in white noise, with signal to noise SN [12]. 
Thus in searching for signals with SN < 20 the HHT 
works best for 6t < 20msec. 
As an example we inject a 20 solar mass black hole 
4.2 a,< 4.S 4.8 5.2 
Time (5) 
FIG. 2: Illustration of Empirical Mode Decomposition (from 
Ref.[4]). (a) time series. (b) average of envelope formed by 
fitting extrema. (c) subtraction of average from time-series. 
Because the result is not vertically symmetric, the waveform 
will be re-sifted until an IMF is identified. 
20 25 30 35 40 45 50x10'~ 
Time [sec] 
FIG. 3: HHT of 20 solar mass black hole binary merger sig- 
nal in white noise. (a) merger signal. (b) signal in noise 
with SN=10. (c),(d),(e) IMFs 2,3,4. (f),(g) instantaneous 
frequency and power derived from the Hilbert Transform of 
IMF3 show the time-frequency-power structure of the merger. 
binary merger and ringdown signal 1131 wit11 SS=lO and 
a time duration of 5 msec into white noise a.l 16 kHz 
sampling rate. The merger signal is shown in Fig.3-(a) 
and the time-series of signal in noise is shown in Fig.3- 
(b). Figures 3-(c),-(d),-(e) show the sifted IMFs; in the 
3rd IMF the signal can be seen, largely separate from the 
noise. 
Figures 3-(f) and 3-(g) show the instantaneous fre- 
quency and power derived from the Hilbert Transform 
of IMF3. The signal frequency and power are clearly 
visible in the Hilbert Transform of IMF3: the frequency 
shows a ramp during the merger and levels off to a con- 
stant during the ringdown, while the power shows the ex- 
pected rise and decay. This level of detail will a.id signal 
identification, allowing comparison of signals from mul- 
tiple detectors. This may be contrasted with the Fourier 
Transform of the merger, which would give the power of 
only 1 point with a 200 Hz frequency spread for the en- 
tire 5 msec. Finally, we note that further noise reduction 
through averaging of the I F  and IA over time is possible, 
as they are oversampled: the plots of Fig.3 a.re sampled 
a t  the LIGO data rate of 16 kHz while the frequencies of 
interest for LIGO analysis are typically below 1 1tNz. 
In analogy with an FFT-based search algorithm that 
looks for excess power in the Fourier power spectrum 
[14], the HHT can be used to search for excess power in 
the. time domain. For a given time record, the summed 
power is computed, and compared to the background 
level; thus clear evidence of a signal is seen in Fig.3- 
(g), the power associated with IMF3. In comparison to 
an FFT-based search, the HHT is most effective for short, 
signals (< 20msec), where FFT analysis tends to lose sen- 
sitivity from time-frequency spreading. Thus the HHT is 
a useful tool, for example, in searching for black hole 
binary merger signals up to 100 solar masses. 
Analysis of Spectral Shapes -To facilitate the analysis 
of data spectral shapes, we use the HHT a.nalogy of tlie 
power spectrum, called the "marginal spectrum." The 
marginal spectrum is defined as: 
where ai(w, t) is the amplitude as a function of the in- 
stantaneous frequency at a given point in time and i is 
an IMF summation index. Thus the marginal spectrum 
is a measure of the total power present a t  a frequency w 
over the time interval T. 
Figure 4 shows the marginal'spectrum of a time-series 
of length 1/16 second with the following Fourier pourer 
spectral shapes: flat, proportional to frequency, and in- 
versely proportional to frequency. The marginal spectra 
are seen to have the expected shapes. Since the sha,pe of 
the marginal spectrum does not depend on the integra- 
tion time, it may be used to examine the detector noise 
stationarity in fine detail. 
Application of HHT to G W  detector characteri,zation 
- The LIGO detector includes a large number of me- 
chanical resonators, servos, cavity optica,l resonances, 
and light and sideband frequencies with differing Q's and 
transient excitation levels [15]. They can produce signals 
which can interact with each other in very transient and 
non-linear ways, for example upconverting noise at low 
frequencies into the signal band [16]. The HHT can be 
used to  examine the detailed dynamics of their interac- 
tions in a way that is unavailable to Fourier decomposi- 
tion. 
We show the capability of the HHT to resolve the fol- 
lowing non-linear, transient waveform: 
10' I o2 I o3 I o4 
Frequency [Hz] 
FIG. 4: HHT marginal spectrum of 1/16sec time:series with 
Fourier spectral shapes that are flat, proportional to fre- 
cjuency, and proportional to l/frequency. 
y = cos [ 2 ~ 6 0 t  + sin (2~40t)I  (5) 
which contains a 60 Hz carrier frequency, phase mod- 
1 0  
M ulated a t  40 Hz. This could be produced by a laser 
. o 5 frequency variation a t  60 Hz coupled t o  seismic isola- 
(Y 
00 tion stack motion a t  40 Hz, for example. This waveform, 
01 which lasts only 40 msec, is shown added to  white noise 
. -0 5 
I- in Fig.&(a). 
-1 0 
0.00 0.02 0.04 0.06 0.08 0.10 0.12 
Time [sec] 
FIG. 5: Time-frequency analysis of simulated instrumental 
artifact (Eq.(5)). a) Time-series of frequency modulated, 
40 msec transienr; in white noise. b) IF and IA from HHT: 
the modulation of the 60 Hz carrier from 20 Hz to 100 Hz 
is clearly visible. c)  FFT-based spectrogram: time-frequency 
spreading blurs the modulation so that detail of its time struc- 
ture is lost. 
After application of the HHT, the modulation of the 
carrier frequency from 20 Hz to  100 Hz is clearly visi- 
ble in Fig.5-(b), where the IA is shown as a function of 
color. In Fig.5-(b) the process of the upconversion of 
noise components a t  40 and 60 Hz (outside the LIGO 
signal band), into components a t  100 Hz (within the sig- 
nal band), is clearly revealed. Resolution a t  this level 
will aid in understanding transient experimental artifacts 
and distinguishing them from real signals. In contrast, 
Fig.5-(c) shows the FFT-based spectrogram, where time- 
frequency spreading blurs the modulation. In this figure 
the detail of the time structure of the modulation is lost. 
Summary - We have presented the application of the 
Hilbert-Huang Transform to  the analysis of GW data. 
The high time-frequency resolution of the HHT, which 
is not limited by time-frequency spreading, will have im- 
portant application to both signal detection and detector 
characterization. 
We acknowledge helpful comments on this manuscript 
by Malik Rakhmanov, Peter Shawhan, and Steve 
Ritz, and early programming assistance from Sean 
McWilliams. 
[I] B. Barish and R. Weiss, Phys. Today 52, 44 (1999). (2004). 
[2] A. Abramovici et al., Science 256, 325 (1992). [6] V. Iqovak et al., BioMedical Engineering Online 3, 39 
[3] K. Darrzmann et al., LISA Pre-Phase A Report - Second (2004). 
Edztzon, MPQ233, July (1998). [7] Hilbert-Huang Transform and its Applications, N .  Huang [a] N. Huang et al., Proc. Roy. Soc. Lond. 454, 903 (1998). et al. (Eds.), World Scientific Publishing, 2005. 
[5] J. Ya,g et al., ASCE J. Engineering Mechanics 130, 85 [8] Time-frequency Analysis, L. Cohen, Prentice Hall, En- 
glewood Cliffs, KJ, pg. 35-36, 1995. Fourier transform of the signal and Sj~vj 1s ihe ipecsrai 
[9] A. Vecchio, Phys. Rev. D 70, 042001 (2004). density of the noise. 
[lo] L. Blanchet, Living Reviews in Relativity, 1131 J .  Baker et al., Phys. Rev. Lett. 96, 111102 (2006). 
www.livingreviews.org/lrr-2006-4, p.91. [14] W. Anderson et al., Phys. Rev. D 63,  042003 (2001). 
[ll] N. Huang et al., Proc. Roy. Soc. Lond. A 460, 1597 [15] B. Abbott et al., Nucl. Instrum. Meth. A 517, (2004) 
(2004). [16] J. Camp et al., J. Opt. Soc. Am. A 17, 120 (2000). 
[12] ( S N ) ~  = (1/27r) j'Tw(lV(w)/2/S(w))dw where V(w) is the 


Application of the Hilbert-Huang Transform to the Search for Gravitational Waves

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server