analyse-ligo-data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieve the LIGO data for the first GW detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pycbc.catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merger = pycbc.catalog.Merger('GW150914')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detectors = [ 'L1', 'H1' ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "strain = {}\n",
    "for detector in detectors:\n",
    "    strain[detector] = merger.strain(detector)\n",
    "    print(\"Loaded {} strain.\".format(detector))\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figstyle_timeseries = {\n",
    "    'figsize': (15, 3)\n",
    "}\n",
    "t_zoom = (merger.time - 0.2, merger.time + 0.1)\n",
    "plt.figure(**figstyle_timeseries)\n",
    "for detector in detectors:\n",
    "    # Zoom around the merger\n",
    "    zoom = strain[detector].time_slice(*t_zoom)\n",
    "    plt.plot(zoom.sample_times, zoom, label=detector)\n",
    "plt.title('Fig. 1: GW150914 data')\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Strain\");\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate the detector sensitivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psd = {}\n",
    "figstyle_psd = {\n",
    "    'figsize': (9, 6)\n",
    "}\n",
    "plt.figure(**figstyle_psd)\n",
    "for detector in detectors:\n",
    "    # Compute the detector PSD\n",
    "    psd[detector] = strain[detector].psd(4)\n",
    "    plt.loglog(psd[detector].sample_frequencies, np.sqrt(psd[detector]), label=detector)\n",
    "plt.xlim(1e1, 1e3)\n",
    "plt.legend()\n",
    "plt.title('Fig. 2: Detector sensitivity')\n",
    "plt.ylabel(\"Amplitude spectral density\")\n",
    "plt.xlabel(\"Frequency [Hz]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Whiten and bandpass the strain data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.psd import interpolate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def whiten(strain, psd):\n",
    "    psd = interpolate(psd, strain.delta_f)\n",
    "    return (strain.to_frequencyseries() / np.sqrt(psd)).to_timeseries()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "strain_whitened = {}\n",
    "plt.figure(**figstyle_timeseries)\n",
    "for detector in detectors:\n",
    "    # Whiten the strain data\n",
    "    strain_whitened[detector] = whiten(strain[detector], psd[detector])\n",
    "    zoom = strain_whitened[detector].time_slice(*t_zoom)\n",
    "    plt.plot(zoom.sample_times, zoom, label=detector)\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Strain\")\n",
    "plt.title(\"Fig. 3: Whitened data\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(**figstyle_timeseries)\n",
    "strain_bandpassed = {}\n",
    "for detector in detectors:\n",
    "    # Bandpass the whitened strain data\n",
    "    strain_bandpassed[detector] = strain_whitened[detector].highpass_fir(frequency=30, order=512).lowpass_fir(frequency=250, order=512)\n",
    "    zoom = strain_bandpassed[detector].time_slice(*t_zoom)\n",
    "    plt.plot(zoom.sample_times, zoom, label=detector)\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Strain\")\n",
    "plt.title(\"Fig. 4: Whitened and bandpassed data\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matched filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.psd import inverse_spectrum_truncation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just apply a highpass filter and crop to make sure we get no boundary artifacts\n",
    "strain_preconditioned = strain['H1'].highpass_fir(15.0, 512).crop(2, 2)\n",
    "\n",
    "# Compute the PSD\n",
    "psd_preconditioned = strain_preconditioned.psd(4)\n",
    "# Interpolate to the strain frequency sampling\n",
    "psd_preconditioned = interpolate(psd_preconditioned, strain_preconditioned.delta_f)\n",
    "# Cut the PSD where we applied the highpass filter\n",
    "psd_preconditioned = inverse_spectrum_truncation(psd_preconditioned, 4 * strain_preconditioned.sample_rate, low_frequency_cutoff=15)\n",
    "\n",
    "plt.figure(**figstyle_timeseries)\n",
    "zoom = whiten(strain_preconditioned, psd_preconditioned).time_slice(*t_zoom)\n",
    "plt.plot(zoom.sample_times, zoom, label=detector)\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Strain\")\n",
    "plt.title(\"Fig. 5: Data for matched filtering\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.waveform import get_td_waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 30 # Component masses of the black-hole binary (in solar masses)\n",
    "template_hplus, template_hcross = get_td_waveform(\n",
    "    approximant=\"SEOBNRv4\",\n",
    "    mass1=M,\n",
    "    mass2=M,\n",
    "    delta_t=strain_preconditioned.delta_t,\n",
    "    f_lower=20\n",
    ")\n",
    "\n",
    "plt.figure(**figstyle_timeseries)\n",
    "plt.title(\"Fig. 6: An equal-mass $M={:.1f}M_\\odot$ BBH waveform template\".format(M))\n",
    "plt.plot(template_hplus.sample_times, template_hplus);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.filter import (matched_filter as pycbc_matched_filter, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matched_filter(template, data, psd):\n",
    "    # Shift template around so that it has the same length as data, but keep its merger at time 0.\n",
    "    template.resize(len(strain_preconditioned))\n",
    "    template = template.cyclic_time_shift(template.start_time)\n",
    "\n",
    "    # Perform the matched filtering\n",
    "    matched_snr = pycbc_matched_filter(\n",
    "        template=template,\n",
    "        data=data,\n",
    "        psd=psd,\n",
    "        low_frequency_cutoff=20,\n",
    "    ).crop(8, 4)\n",
    "\n",
    "    # Analyse the SNR\n",
    "    peak_index = abs(matched_snr).numpy().argmax()\n",
    "    peak_snr = matched_snr[peak_index]\n",
    "    merger_time = matched_snr.sample_times[peak_index]\n",
    "    \n",
    "    # Shift the template to the merger time\n",
    "    template_aligned = template.cyclic_time_shift(merger_time - data.start_time)\n",
    "\n",
    "    # Scale the template so that it would have SNR 1 in this data\n",
    "    template_aligned /= sigma(template_aligned, psd=psd, low_frequency_cutoff=20.0)\n",
    "\n",
    "    # Scale the template amplitude and phase to the peak value\n",
    "    template_aligned = (template_aligned.to_frequencyseries() * peak_snr).to_timeseries()\n",
    "\n",
    "    # Align the start time of the template\n",
    "    template_aligned.start_time = data.start_time\n",
    "    \n",
    "    return template_aligned, matched_snr, peak_snr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matched_template, matched_snr, _ = matched_filter(template_hplus, strain_preconditioned, psd_preconditioned)\n",
    "\n",
    "plt.figure(**figstyle_timeseries)\n",
    "plt.plot(matched_snr.sample_times, abs(matched_snr))\n",
    "plt.ylabel('Signal-to-noise ratio')\n",
    "plt.title(\"Fig. 7: Matched-filtering SNR\")\n",
    "plt.xlabel('Time [s]');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(**figstyle_timeseries)\n",
    "\n",
    "zoom_data = strain_bandpassed['H1'].time_slice(*t_zoom)\n",
    "plt.plot(zoom_data.sample_times, zoom_data, label=\"Data\")\n",
    "\n",
    "zoom_template = whiten(matched_template, psd_preconditioned).time_slice(*t_zoom)\n",
    "plt.plot(zoom_template.sample_times, zoom_template, label=\"Matched template\", color='green')\n",
    "\n",
    "plt.ylabel('Strain')\n",
    "plt.xlabel('Time [s]')\n",
    "plt.title(\"Fig. 8: Matched waveform template aligned with the data\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(**figstyle_timeseries)\n",
    "\n",
    "zoom_data = strain_bandpassed['H1'].time_slice(*t_zoom)\n",
    "plt.plot(zoom_data.sample_times, zoom_data, label=\"Data\")\n",
    "\n",
    "snrs = []\n",
    "\n",
    "masses = np.linspace(25, 45, 20)\n",
    "for mass in masses:\n",
    "    template, _ = get_td_waveform(\n",
    "        approximant=\"SEOBNRv4\",\n",
    "        mass1=mass,\n",
    "        mass2=mass,\n",
    "        delta_t=strain_preconditioned.delta_t,\n",
    "        f_lower=20\n",
    "    )\n",
    "    matched_template, _, peak_snr = matched_filter(template, strain_preconditioned, psd_preconditioned)\n",
    "    snrs.append(abs(peak_snr))\n",
    "    \n",
    "    zoom_template = whiten(matched_template, psd_preconditioned).time_slice(*t_zoom)\n",
    "    plt.plot(zoom_template.sample_times, zoom_template)\n",
    "\n",
    "plt.ylabel('Strain')\n",
    "plt.xlabel('Time [s]')\n",
    "plt.title(\"Fig. 9: Parameter estimation with many waveform templates\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(masses, snrs, marker='.', ls='dotted')\n",
    "plt.xlabel(\"Mass\")\n",
    "plt.title(\"Fig. 10: Parameter estimation with many waveform templates\")\n",
    "plt.ylabel(\"SNR\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}