{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# GWpy Tutorial\n\nGWpy is a collaboration-driven Python package providing tools for studying data from ground-based gravitational-wave detectors.\n\nGWpy provides a user-friendly, intuitive interface to the common time-domain and frequency-domain data produced by the LIGO and Virgo instruments and their analysis, with easy-to-follow tutorials at each step.\n\n- https://gwpy.github.io/docs/v0.1/timeseries/index.html\n- https://www.kaggle.com/mistag/data-preprocessing-with-gwpy","metadata":{"execution":{"iopub.status.busy":"2021-08-10T14:44:20.645292Z","iopub.execute_input":"2021-08-10T14:44:20.645816Z","iopub.status.idle":"2021-08-10T14:44:20.656387Z","shell.execute_reply.started":"2021-08-10T14:44:20.645735Z","shell.execute_reply":"2021-08-10T14:44:20.655028Z"}}},{"cell_type":"code","source":"!python -m pip install gwpy\n!pip install astropy==4.2.1","metadata":{"execution":{"iopub.status.busy":"2023-03-23T17:10:10.599491Z","iopub.execute_input":"2023-03-23T17:10:10.599882Z","iopub.status.idle":"2023-03-23T17:10:42.592002Z","shell.execute_reply.started":"2023-03-23T17:10:10.599849Z","shell.execute_reply":"2023-03-23T17:10:42.590564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"GWpy provides the TimeSeries object as a way of representing time series data. The TimeSeries is built on top of the numpy.ndarray, and so many methods and applications of this object should be familiar to numpy users.","metadata":{}},{"cell_type":"code","source":"import numpy as np\n\nfrom gwpy.timeseries import TimeSeries\n\nfrom gwpy.plot import Plot\nfrom scipy import signal\nfrom sklearn.preprocessing import MinMaxScaler\nfrom PIL import Image\nfrom matplotlib import pyplot as plt\nfrom sklearn.preprocessing import normalize\nfrom sklearn.preprocessing import StandardScaler","metadata":{"execution":{"iopub.status.busy":"2023-03-23T17:22:59.589262Z","iopub.execute_input":"2023-03-23T17:22:59.589713Z","iopub.status.idle":"2023-03-23T17:22:59.732059Z","shell.execute_reply.started":"2023-03-23T17:22:59.589676Z","shell.execute_reply":"2023-03-23T17:22:59.730710Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# TimeSeries.fetch_open_data()\n\nthe TimeSeries.fetch_open_data() method as an interface to the GWOSC data archive, requiring users to provide a minimum of information in order to access data.\n\nFor example, to fetch 30 seconds of strain data around the first ever gravitational-wave detection (GW150914), you need to give the prefix of the relevant observatory ('H1' for the LIGO Hanford Observatory, 'L1' for LIGO Livingston), and the start and end times of your query.","metadata":{}},{"cell_type":"markdown","source":"Let's consider data from the 4km detector at LIGO Livingston Observatory (LLO), known as **L1** and call TimeSeries.fetch_open_data() to download the calibrated GW strain data in that interval:","metadata":{}},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\n#dimensionless\ndata = TimeSeries.fetch_open_data('L1', 1126259446, 1126259478) \ndata\nplot = data.plot()\nplot.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-23T17:23:02.567763Z","iopub.execute_input":"2023-03-23T17:23:02.568240Z","iopub.status.idle":"2023-03-23T17:23:10.363558Z","shell.execute_reply.started":"2023-03-23T17:23:02.568195Z","shell.execute_reply":"2023-03-23T17:23:10.362436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Whitnening a TimeSeries\n \nMost data recorded from a gravitational-wave interferometer carry information across a wide band of frequencies, typically up to a few kiloHertz, but often it is the case that the low-frequency amplitude dwarfs that of the high-frequency content, making discerning high-frequency features difficult.\n\nWe employ a technique called ‘whitening’ to normalize the power at all frequencies so that excess power at any frequency is more obvious.\n\nWe demonstrate below with an auxiliary signal recording transmitted power in one of the interferometer arms, which recorded two large glitches with a frequency of around 5-50Hz.","metadata":{}},{"cell_type":"code","source":"def read_file(fname):\n    data = np.load(fname)\n    data = rescale(data)\n    d1 = TimeSeries(data[0,:], sample_rate=2048)\n    d2 = TimeSeries(data[1,:], sample_rate=2048)\n    d3 = TimeSeries(data[2,:], sample_rate=2048)  \n    return d1, d2, d3\n\n\ndef rescale(data):\n    # Normalize to [-1,1]\n    data = (data - np.mean(data, axis=1).reshape(3,-1)) / (np.max(data, axis=1) - np.min(data, axis=1)).reshape(3,-1) * 2\n    return data\n\n    \ndef plot_time_data(d1, d2, d3, title=''):\n    plot = Plot(d1, d2, d3, separate=True, sharex=True, figsize=[12, 8])\n    ax = plot.gca()   \n    ax.set_title(title)\n    ax.set_xlim(0,2)\n    ax.set_xlabel('Time [s]')\n    plot.show()\n        \n        \ndef preprocess(d1, d2, d3, bandpass=False, normalize=True, lf=35, hf=350):\n    white_d1 = d1.whiten(window=(\"tukey\",0.2))\n    white_d2 = d2.whiten(window=(\"tukey\",0.2))\n    white_d3 = d3.whiten(window=(\"tukey\",0.2))\n        \n    if bandpass: # bandpass filter Design a band-pass filter for the given cutoff frequencies\n        bp_d1 = white_d1.bandpass(lf, hf) \n        bp_d2 = white_d2.bandpass(lf, hf)\n        bp_d3 = white_d3.bandpass(lf, hf)\n        return bp_d1, bp_d2, bp_d3\n    else: # only whiten\n        return white_d1, white_d2, white_d3\n    \nsample_signals    = read_file('/kaggle/input/g2net-gravitational-wave-detection/train/0/0/0/00000e74ad.npy')\nsample_no_signals = read_file('/kaggle/input/g2net-gravitational-wave-detection/train/0/0/0/00001f4945.npy')","metadata":{"execution":{"iopub.status.busy":"2023-03-23T15:35:23.070009Z","iopub.execute_input":"2023-03-23T15:35:23.070229Z","iopub.status.idle":"2023-03-23T15:35:23.115065Z","shell.execute_reply.started":"2023-03-23T15:35:23.070206Z","shell.execute_reply":"2023-03-23T15:35:23.114496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_time_data(*sample_signals, title='raw')","metadata":{"execution":{"iopub.status.busy":"2023-03-23T15:35:23.115926Z","iopub.execute_input":"2023-03-23T15:35:23.116237Z","iopub.status.idle":"2023-03-23T15:35:26.260885Z","shell.execute_reply.started":"2023-03-23T15:35:23.116207Z","shell.execute_reply":"2023-03-23T15:35:26.259455Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_time_data(*preprocess(*sample_signals, bandpass=False), title='whiten')","metadata":{"execution":{"iopub.status.busy":"2023-03-23T15:35:26.263193Z","iopub.execute_input":"2023-03-23T15:35:26.263738Z","iopub.status.idle":"2023-03-23T15:35:29.609787Z","shell.execute_reply.started":"2023-03-23T15:35:26.263712Z","shell.execute_reply":"2023-03-23T15:35:29.608770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_time_data(*preprocess(*sample_signals, bandpass=True),  title='bandpass')","metadata":{"execution":{"iopub.status.busy":"2023-03-23T15:35:29.610831Z","iopub.execute_input":"2023-03-23T15:35:29.611144Z","iopub.status.idle":"2023-03-23T15:35:32.650119Z","shell.execute_reply.started":"2023-03-23T15:35:29.611121Z","shell.execute_reply":"2023-03-23T15:35:32.649118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.array( sample_signals[0].to_value().shape ) ","metadata":{"execution":{"iopub.status.busy":"2023-03-23T15:35:32.651541Z","iopub.execute_input":"2023-03-23T15:35:32.651829Z","iopub.status.idle":"2023-03-23T15:35:32.660177Z","shell.execute_reply.started":"2023-03-23T15:35:32.651805Z","shell.execute_reply":"2023-03-23T15:35:32.657837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Q-Transform\n\nOne of the most useful tools for filtering and visualising short-duration features in a TimeSeries is the Q-transform. This is regularly used by the Detector Characterization working groups of the LIGO Scientific Collaboration and the Virgo Collaboration to produce high-resolution time-frequency maps of transient noise (glitches) and potential gravitational-wave signals.\n","metadata":{}},{"cell_type":"code","source":"# Constant-Q transform (CQT) here refers to a technique that transforms a time-domain signal x(n) \n# into the time- frequency domain so that the center frequencies of the frequency bins are \n# geometrically spaced and their Q-factors are all equal\n\nr1, r2, r3 = read_file('../input/g2net-gravitational-wave-detection/train/0/0/0/000a5b6e5c.npy') # this signal has target=1\np1, p2, p3 = preprocess(r1, r2, r3)\nhq = p2.q_transform(qrange=(16,32), frange=(30,400), logf=True, whiten=False)\nfig4 = hq.plot(figsize=[12, 10])\nax = fig4.gca()\nfig4.colorbar(label=\"Normalised energy\")\nax.grid(False)\nax.set_yscale('log')\nax.set_xlabel('Time [s]');","metadata":{"execution":{"iopub.status.busy":"2023-03-23T15:35:32.661876Z","iopub.execute_input":"2023-03-23T15:35:32.662136Z","iopub.status.idle":"2023-03-23T15:35:34.673871Z","shell.execute_reply.started":"2023-03-23T15:35:32.662111Z","shell.execute_reply":"2023-03-23T15:35:34.673095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This algorithm was used to visualise the first ever gravitational-wave detection GW150914, so we can reproduce that result (bottom panel of figure 1) here.","metadata":{}},{"cell_type":"code","source":"from gwpy.timeseries import TimeSeries\ndata = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)","metadata":{"execution":{"iopub.status.busy":"2023-03-23T17:17:45.675056Z","iopub.execute_input":"2023-03-23T17:17:45.675481Z","iopub.status.idle":"2023-03-23T17:17:54.521797Z","shell.execute_reply.started":"2023-03-23T17:17:45.675438Z","shell.execute_reply":"2023-03-23T17:17:54.520529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"qspecgram = data.q_transform(outseg=(1126259462.2, 1126259462.5))","metadata":{"execution":{"iopub.status.busy":"2023-03-23T17:18:03.234934Z","iopub.execute_input":"2023-03-23T17:18:03.235352Z","iopub.status.idle":"2023-03-23T17:18:03.988324Z","shell.execute_reply.started":"2023-03-23T17:18:03.235312Z","shell.execute_reply":"2023-03-23T17:18:03.987214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot = qspecgram.plot(figsize=[8, 4])\nax = plot.gca()\nax.set_xscale('seconds')\nax.set_yscale('log')\nax.set_ylim(20, 500)\nax.set_ylabel('Frequency [Hz]')\nax.grid(True, axis='y', which='both')\nax.colorbar(cmap='viridis', label='Normalized energy')\nplot.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-23T17:18:13.231002Z","iopub.execute_input":"2023-03-23T17:18:13.231435Z","iopub.status.idle":"2023-03-23T17:18:17.648495Z","shell.execute_reply.started":"2023-03-23T17:18:13.231399Z","shell.execute_reply":"2023-03-23T17:18:17.647350Z"},"trusted":true},"execution_count":null,"outputs":[]}]}