8
$\begingroup$

Has anyone looked at the LIGO data using Mathematica or know how to do this? I'm not familiar with the hdf5 format but I've been able to import the data, consisting primarily of 131,072 floating point numbers representing the strain at different time points.

h1url = "https://losc.ligo.org/s/events/GW150914/H-H1_LOSC_4_V1-1126259446-32.hdf5"
Import[h1url]

shows a number of datasets including "/strain/Strain"

h1losc = Import[h1url, {"Datasets", "/strain/Strain"}];

This yields the main block of data, as a list with 131,072 floating point elements.

There appears to be no other large data in the file. This is confirmed by the output of h5dump. However, the Python tutorial I'm following manages to also extract 131,072 time values from the .hdf5 files from the two detectors. In fact it seems that there are 3 types of information associated with each of the 131,072 parcels of data rather than the single floating point number I can see.

The data can be plotted already using the index as horizontal axis:

ListLogLogPlot[h1losc]

output of previous command as a graphic

My first question is: is it possible that the import is discarding the timestamp data? And if so, (though with apologies this question is more about hdf5dump than Mathematica), why does this not show up in the output of hdf5dump?

Secondly, how to generate the timestamps for the data so they can be analysed, transformed and plotted? The start and finish times can be found with

start = Import[h1loc, {"Datasets", "/meta/GPSstart"}]
duration = Import[h1url, {"Datasets", "/meta/Duration"}]
finish = start + duration

Mathematica graphics

tevent = 1126259462.422`16

The main event at was only a few seconds long so I want to extract the data points that correspond to a 5-second window around that event. How can I do this?

$\endgroup$
  • 2
    $\begingroup$ Since this is a Mathematica forum, you should provide the code that you have actually tried so that we don't have to guess what data you're talking about. Otherwise, the question will be closed due to lack of Mathematica-specific detail. The linked web page has numerous hdf5 files, all with the same prefix quoted in the title. $\endgroup$ – Jens Feb 20 '16 at 18:32
  • 1
    $\begingroup$ sorry, Jens, you are quite right. Code now supplied. $\endgroup$ – fairflow Feb 20 '16 at 23:13
  • 2
    $\begingroup$ Since posting this I found that the samples are evenly distributed in the 32 second data window so it is possible to recover the timestamps. I can post an answer to this and hopefully this will mean my question won't be closed. $\endgroup$ – fairflow Feb 21 '16 at 12:23
  • $\begingroup$ Have you tried TemporalData? $\endgroup$ – FredrikD Feb 21 '16 at 12:46
11
$\begingroup$

My first question is: is it possible that the import is discarding the timestamp data?

The relevant timing data is contained in the attributes of the /Strain/strain dataset. These can be extracted using:

H1url = "https://losc.ligo.org/s/events/GW150914/H-H1_LOSC_4_V1-1126259446-32.hdf5";
strainH1 = Import[H1url, {"Datasets",   "/strain/Strain"}];
attrsH1  = Import[H1url, {"Attributes", "/strain/Strain"}]

Then attrsH1 is a list of replacement rules:

{"Npoints" -> 131072, "Xlabel" -> "GPS time", 
 "Xspacing" -> 0.000244141, "Xstart" -> 1126259446, "Xunits" -> "second",
 "Ylabel" -> "Strain", "Yunits" -> "\.00"}

Secondly, how to generate the timestamps for the data so they can be analysed, transformed and plotted?

Assigning the relevant attributes to the following variables...

{t0, dt, n} = {"Xstart", "Xspacing", "Npoints"} /. attrsH1;

... you can then generate the list of absolute times corresponding to each sample in strainH1.

time = t0 + (Range[n] - 1) dt;

The main event at was only a few seconds long so I want to extract the data points that correspond to a 5-second window around that event. How can I do this?

tEvent = 1126259462.39;  (* 2015-09-14T09:50:45 GMT *)
tWindow = 5;
tFilter = Flatten[Position[tWindow - Abs[time - tEvent], _Real?Positive]];

The list tFilter contains the indices of strainH1 and time that correspond to the interval 5 seconds before and after tEvent.

Note: Here I have used the value for tEvent quoted in the data release for event GW150914, and not that in the Python tutorial you mentioned. (The latter is 30ms later than the former.)

This can be used to index the lists, for example:

ListPlot[
 Transpose[{time[[tFilter]] - tEvent, strainH1[[tFilter]]}],
 Joined -> True,
 Frame -> True,
 FrameLabel -> {"time (s) since " <> ToString[AccountingForm[tEvent, 12]] <>
 " (GPS time)", "strain"}
]

This produces the plot of the raw data shown below.

Raw data from event GW150914

To see the gravitational wave signal, you can repeat the steps in the Python tutorial, using the equivalent Mathematica code to perform the signal whitening, band-pass filtering, etc. For example, PeriodogramArray[strainH1, Round[1/dt], 1, HannWindow] will give you something like the power spectrum Pxx_H1 in that tutorial (shown below for the 4,096 samples/second files).

Amplitude spectral density of LOSC_4_V1-1126257414-4096 data from http://losc.ligo.org.

$\endgroup$
  • 1
    $\begingroup$ Great stuff Russell, just what I needed. You also answered some other questions I was planning to ask! I needed to do h1peri = PeriodogramArray[h1losc, Round[1/dt], 1, HannWindow]; to get PeriodogramArray to work correctly; then I did a ListLogLogPlot $\endgroup$ – fairflow Feb 29 '16 at 12:37
  • $\begingroup$ Thanks, @fairflow. The partition length (equal to the sample rate in this example, i.e. each partition is of 1 second duratiion) must be an integer. I've modified my example PeriogramArray call to reflect this. $\endgroup$ – Russell Anderson Mar 1 '16 at 21:04
6
$\begingroup$

TemporalData has a straightforward way of adding timestamps to data. Using your variables:

td = TemporalData[hlosc, {start, finish}]

enter image description here

Extracting a 5 second interval around tevent:

ListLogPlot[td, PlotRange -> {{tevent - 5, tevent + 5}, Automatic}]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ That helps, thanks! My solution without TemporalData was simple enough but it now looks redundant. $\endgroup$ – fairflow Feb 21 '16 at 14:00
  • $\begingroup$ I'm still puzzled about the apparently missing data though. I've asked losc@ligo.org about this and await a reply. $\endgroup$ – fairflow Feb 21 '16 at 20:10
  • 1
    $\begingroup$ I think that this method is fine, but the solution below is closer to what I wanted! Thanks @FredrikD $\endgroup$ – fairflow Feb 26 '16 at 15:53
  • $\begingroup$ No problem. Much more hdf5 info in Russels answer. $\endgroup$ – FredrikD Feb 27 '16 at 13:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.