Data_Science_Problem_Sets/Signal_Processing_with_LIGO/Signal_Processing_and_Visualization_with_LIGO.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Science Problem Set 1: \n",
    "## Signal Processing and Data Visualization of a LIGO Binary Black Hole Merger\n",
    "\n",
    "Josh Dillon\n",
    "\n",
    "Based on [LIGO Open Science Center tutorials](https://www.gw-openscience.org/tutorials/). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In 2015, the [Laser Interferometer Gravitational-Wave Observatory (LIGO)](https://en.wikipedia.org/wiki/LIGO) made the [first direction detection](https://en.wikipedia.org/wiki/First_observation_of_gravitational_waves) of gravitational waves---waves in the fabric of spacetime itself propagating at the speed of light. What they saw was the result of two black holes, each roughly 30 times the mass of the sun approximately 400 Mpc (1.3 billion light years) away. That event was so violent that approximately 3 solar masses were converted into energy (via $E = mc^2$) and released as gravitational radiation. \n",
    "\n",
    "![](data/BBH.gif)\n",
    "**Figure: Binary black holes inspiraling radiate away energy as gravitational radiation. Credit: [LIGO Caltech](https://www.ligo.caltech.edu/page/ligos-ifo)**\n",
    "\n",
    "\n",
    "These were detected by two laser interferometers, one in Hanford, Washington and the other in Livingston, Louisiana.\n",
    "\n",
    "![](data/Hanford.jpg)\n",
    "\n",
    "**Figure: Each LIGO detector featues two 4 km arms (LIGO Hanford shown here. Credit: [LIGO Caltech](https://www.ligo.caltech.edu/page/ligos-ifo)**\n",
    "\n",
    "This discovery was recognized with many awards, including the 2017 Nobel Prize in Physics.\n",
    "\n",
    "The purpose of this problem set is to explore topics in data visualization (and sonification) and signal processing using data from that event, known as GW150914 (since it occurred on September 14, 2015)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:18:45.648524Z",
     "start_time": "2020-11-24T22:18:44.734772Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "from scipy import signal\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data:\n",
    "\n",
    "I've already downloaded the relevant data and saved it to a npz file with everything you need. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:18:46.382370Z",
     "start_time": "2020-11-24T22:18:46.372270Z"
    }
   },
   "outputs": [],
   "source": [
    "data = np.load('data/ligo_data.npz')\n",
    "times = data['times'] # Times in seconds since the gravitational wave event. \n",
    "hanford_strain = data['hanford_strain'] # Strain amplitude measured at LIGO Hanford, unitless\n",
    "livingston_strain = data['livingston_strain'] # Strain amplitude measured at LIGO Livingston, unitless"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1: Plotting Gravitational Wave Strain\n",
    "\n",
    "Gravitational waves can be thought of as a stretching and squeezing of spacetime. This is usually described by a **strain**, $h$, which is the fractional change in length as a gravitational wave passes by. \n",
    "\n",
    "$h = \\Delta L / L$\n",
    "\n",
    "Here's an example of what a ring of \"stationary\" particles look like as a gravitational wave passes for a strain of $h = .5$. \n",
    "\n",
    "\n",
    "\n",
    "<a href=https://en.wikipedia.org/wiki/Gravitational_wave> <img src=\"data/strain.gif\" /> </a>\n",
    "\n",
    "You can think of LIGO as sitting at the center of this ring, with mirrors at 12 o'clock and 3 o'clock. LIGO is designed to measure the difference of the lengths of the two arms very precisely as a fuction of time.\n",
    "\n",
    "<a href=https://www.ligo.caltech.edu/page/ligos-ifo> <img src=\"data/ligo_layout.jpg\" style=\"width:600px;\"/> </a>\n",
    "\n",
    "However, instead of strains of order 1, LIGO is looking for is looking for strains $h \\lesssim 10^{-21}$. \n",
    "\n",
    "### Plot the strain amplitudes measured at the two detectors as a function of time. Don't forget to label your axes and include a legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:18:52.180218Z",
     "start_time": "2020-11-24T22:18:52.177890Z"
    }
   },
   "outputs": [],
   "source": [
    "# Your code here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2: Looking at the Amplitude Strain Density\n",
    "\n",
    "You should see that both detectors measured a complicated signal with content at many different frequencies. However, most of that is not the gravitational wave signal. What you're seeing is various sources of noise in the measurement. To understand better the systematics, it is useful to look at the Fourier transform of the signal. \n",
    "\n",
    "The traditional way of doing this in gravitational wave physics is is to look at a quantity called the \"Amplitude Strain Density\" which is a basically the square root of the power spectrum. There are several ways of estimating the spectral content of a signal using a [periodogram](https://en.wikipedia.org/wiki/Periodogram). One useful one, [Welch's method](https://en.wikipedia.org/wiki/Welch%27s_method) is built into `scipy.signal`, which breaks the signal up into chunks, computes the power, and then averages that together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:18:53.583835Z",
     "start_time": "2020-11-24T22:18:53.541384Z"
    }
   },
   "outputs": [],
   "source": [
    "delta_t = (times[1] - times[0])\n",
    "sampling_freq = delta_t**-1\n",
    "\n",
    "# this treats each second as a different \"segment\" or chunk, though you can play around with different nperseg\n",
    "freqs, power_hanford = scipy.signal.welch(hanford_strain, sampling_freq, nperseg=1*sampling_freq, detrend=False)\n",
    "freqs, power_livingston = scipy.signal.welch(livingston_strain, sampling_freq, nperseg=1*sampling_freq, detrend=False)\n",
    "\n",
    "# These interpolators will be useful later.\n",
    "from scipy.interpolate import interp1d\n",
    "hanford_power_interpolator = interp1d(freqs, power_hanford)\n",
    "livingston_power_interpolator = interp1d(freqs, power_livingston)\n",
    "\n",
    "plt.figure()\n",
    "# Note we're taking the square root, since we want something that's directly proportional to the strain amplitude.\n",
    "# If you double the strain, the power quadruples.\n",
    "plt.loglog(freqs, np.sqrt(power_hanford), label='Hanford') \n",
    "plt.loglog(freqs, np.sqrt(power_livingston), label='Livingston')\n",
    "# plt.loglog(freqs, np.sqrt(noise_power_model(freqs)), 'k--', label='Noise Power Model')\n",
    "plt.xlabel('Frequency (Hz)')\n",
    "plt.ylabel('Amplitude Strain Density (strain/$\\sqrt{Hz}$)')\n",
    "plt.xlim(20, 2000)\n",
    "plt.ylim([1e-25, 1e-18])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This produces \"cleaner\" results than simply taking a Fourier transform without any fancy chunking or averaging, but it's still useful to explore that method to get a sense of what's going on.\n",
    "\n",
    "### Plot the Amplitude Strain Density by taking a simple FFT using `np.fft.fft`:\n",
    "\n",
    "***Hints: you'll probably need to multiply by a [window function](https://en.wikipedia.org/wiki/Window_function) like a Blackman-Harris to get sensible results. That window function should be normalized to 1 so as to not lose power. Think about how to properly normalize the FFT to have sensible physical units; you're approximating an integral over time by a sum, so you'll have to think about how the length element $dt$ is accounted for. Also, check out how [`np.fft.fft` is normalized here](https://numpy.org/doc/stable/reference/routines.fft.html#normalization).***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:18:59.109860Z",
     "start_time": "2020-11-24T22:18:59.107663Z"
    }
   },
   "outputs": [],
   "source": [
    "# your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 3: Interlude: Colored Noise, Whitening, and Sonification\n",
    "\n",
    "An important feature of the different sources of noise in the LIGO data is that it is they have different amounts of power at different frequencies. That would be \"white noise.\" Noise with more high frequency power, in analogy with visible light, is \"blue noise.\" Noise with more low frequency power is \"red noise.\" We can simulate noise like that by taking white noise, which is just uncorrelated Gaussian random noise generated with `np.random.randn`, Fourier transforming it, modifying the amplitude in a frequency-dependent way, and the Fourier transforming back. For example, we use the following function to generate some white, red, and blue noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:00.586931Z",
     "start_time": "2020-11-24T22:19:00.568360Z"
    }
   },
   "outputs": [],
   "source": [
    "def colored_noise(times, exponent, min_freq=20):\n",
    "    '''Creates random noise with a spectrum proportional to freq^exponent.'''\n",
    "    white_noise = np.random.randn(len(times))\n",
    "    freqs = np.fft.rfftfreq(len(times), times[1] - times[0]) #rfft is useful when the input signal is real\n",
    "    noise = np.fft.rfft(white_noise)\n",
    "    noise[np.abs(freqs) >= min_freq] *= np.abs(freqs[np.abs(freqs) >= min_freq])**exponent\n",
    "    # Remove power below where you can hear it, avoids division by 0\n",
    "    noise[np.abs(freqs) < min_freq] = 0 \n",
    "    noise = np.fft.irfft(noise)\n",
    "    # Renormalize the noise so it has the same power as white noise\n",
    "    noise /= np.std(noise) \n",
    "    return noise\n",
    "\n",
    "white_noise = colored_noise(times, 0)\n",
    "blue_noise = colored_noise(times, 1)\n",
    "red_noise = colored_noise(times, -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3a: Use `scipy.signal.welch` to plot the amplitude spectral density of white, blue, and red noise, just as we did above with LIGO data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:04.749505Z",
     "start_time": "2020-11-24T22:19:04.747270Z"
    }
   },
   "outputs": [],
   "source": [
    "# your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way of understanding timestreams like these is to transform them into sound. Sound is just pressure waves in the air. What you hear can be expressed as an amplitude as a function of time, which is precisely what a `.wav` file includes. While LIGO doesn't measure sound waves, the frequencies it \"listens\" are in the audible range (roughly 20 Hz to 20 kHz, though the top end generally gets lower as you get older). So for example, we can use the following function to play data vectors, be they simulated noise or LIGO signal, as sound. So turn your sound on and see if you can hear the difference!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:06.253364Z",
     "start_time": "2020-11-24T22:19:06.247753Z"
    }
   },
   "outputs": [],
   "source": [
    "def play_audio(data, times, t_start=None, t_end=None):\n",
    "    '''Convert data to a wav file and play it in the notebook, given times in s.'''\n",
    "    from IPython.display import Audio\n",
    "    from scipy.io import wavfile\n",
    "    sampling_freq = (times[1] - times[0])**-1\n",
    "    if t_start is None:\n",
    "        t_start = times[0]\n",
    "    if t_end is None:\n",
    "        t_end = times[-1]\n",
    "    d = data[(times >= t_start) & (times <= t_end)]\n",
    "    # renornmalize the data so that it uses the dynamic range of a 32-bit wav file\n",
    "    d = np.int32(d / np.max(np.abs(d)) * 2147483647 * 0.99)\n",
    "    wavfile.write('temp.wav', int(sampling_freq), d)\n",
    "    return Audio('temp.wav', autoplay=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:07.417812Z",
     "start_time": "2020-11-24T22:19:07.387324Z"
    }
   },
   "outputs": [],
   "source": [
    "play_audio(white_noise, times, t_start=-2.5, t_end=2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:08.146808Z",
     "start_time": "2020-11-24T22:19:08.106019Z"
    }
   },
   "outputs": [],
   "source": [
    "play_audio(blue_noise, times, t_start=-2.5, t_end=2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:08.697338Z",
     "start_time": "2020-11-24T22:19:08.691303Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "play_audio(red_noise, times, t_start=-2.5, t_end=2.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In signal processing we can use our knowledge the about the statistical properties of the noise to upweight the modes with more signal relative to the noise. In some cases, this is a simple as dividing out by the noise amplitude spectrum in Fourier space. This is the basic idea behind the [Wiener filter](https://en.wikipedia.org/wiki/Wiener_filter) and many more sophisticated techniques. Our goal today is just to get some intuition for this using sonification. Let's pretend for a moment that the following is our \"signal\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:10.547858Z",
     "start_time": "2020-11-24T22:19:10.474943Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from scipy.io import wavfile\n",
    "sampling_rate, song = wavfile.read('data/beatles.wav')\n",
    "song = song / np.max(np.abs(song))\n",
    "ts = np.arange(len(song)) / sampling_rate\n",
    "play_audio(song, ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's add some white noise, some very strong blue noise, and listen again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:11.814049Z",
     "start_time": "2020-11-24T22:19:11.615440Z"
    }
   },
   "outputs": [],
   "source": [
    "blue_noise = colored_noise(ts, 4)\n",
    "white_noise = colored_noise(ts, 0, min_freq=0)\n",
    "# Add very strong, very blue noise to dominate the high frequencies\n",
    "noise = 100 * blue_noise + white_noise\n",
    "noisy_song = 5 * song + noise\n",
    "play_audio(noisy_song, ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The signal is still there, but it's swamped by the noise--just like in LIGO. But if we use the following function to whiten the data given the power spectrum of the noise, we can start teasing out the signal from the noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:13.061127Z",
     "start_time": "2020-11-24T22:19:13.056614Z"
    }
   },
   "outputs": [],
   "source": [
    "def whiten(data, times, noise_power_spectrum_func):\n",
    "    '''Use a function that returns power spectra as a function of frequency\n",
    "    (either interpolated or analytic) to whiten the noise in the data.'''\n",
    "    fft_data = np.fft.rfft(data)\n",
    "    freqs = np.fft.rfftfreq(len(data), d=(times[1] - times[0]))\n",
    "    fft_data /= np.sqrt(noise_power_spectrum_func(freqs))\n",
    "    return np.fft.irfft(fft_data * 1./np.sqrt(1./((times[1] - times[0]*2))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:13.764027Z",
     "start_time": "2020-11-24T22:19:13.700362Z"
    }
   },
   "outputs": [],
   "source": [
    "# calculate new sampling frequency for this audio clip\n",
    "fs = (ts[1] - ts[0])**-1\n",
    "# calculate power spectrum of the song with noise and just the noise\n",
    "freqs, power_noisy_song = scipy.signal.welch(noisy_song, fs, nperseg=fs, detrend=False)\n",
    "freqs, power_noise = scipy.signal.welch(noise, fs, nperseg=fs, detrend=False)\n",
    "# get an interpolated noise power spectrum\n",
    "noise_power_interpolator = interp1d(freqs, power_noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3c: Plot the amplitude spectral density of the the noise, the song clip with the noise, and of the whitened noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:17.791632Z",
     "start_time": "2020-11-24T22:19:17.789528Z"
    }
   },
   "outputs": [],
   "source": [
    "# Your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's still not perfect, but hopefully that sounds _a lot_ better than with the raw noise!\n",
    "\n",
    "And now, back to LIGO..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4: Now let's look at how the LIGO spectrum evolves in time.\n",
    "\n",
    "As we said before, the vast majority of what we see in the spectrum is not the gravitational wave signal. That signal occupies a small fraction of the times and frequencies, so in order to see that we want to further break up the spectrum into different time bins, producing a [spectrogram](https://en.wikipedia.org/wiki/Spectrogram). This is a bit of a funny operation, since it's turning a 1D time series into a 2D plot. This is a achieved by breaking the data into blocks of size `NFFT`, taking the Fast Fourier Transform, then moving on to the next chunk. The larger the chunk, the better the frequency resolution, but the worse the temporal resolution. If this sounds familiar, it's because it's essentially a sort of [uncertainty principle](https://en.wikipedia.org/wiki/Uncertainty_principle#Signal_processing).\n",
    "\n",
    "Better temporal resolution can be achieved by letting subsequent blocks overlap, though this introduces correlations between neighboring times. This makes plots a bit easier to read, but doesn't fundamentally change the information content.\n",
    "\n",
    "Let's start with a simplified example, where we have a sinusoidal signal $sin(\\omega t)$ whose angular frequency (in radians per second) is given by \n",
    "\n",
    "$\\omega = 2 \\pi \\times 4 \\left(\\frac{t}{1 s}\\right)^2 $.\n",
    "\n",
    "Here we show what this signal looks like in time and as a spectrogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:20.756228Z",
     "start_time": "2020-11-24T22:19:20.576697Z"
    }
   },
   "outputs": [],
   "source": [
    "spectrogram_interval = 1 / 8 # seconds\n",
    "NFFT = int(spectrogram_interval / delta_t)\n",
    "noverlap = int(NFFT*63./64) # the closer this value is to NFFT, the higher time resolution, but the more correlated the points are\n",
    "window = np.blackman(NFFT)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10,6))\n",
    "omega = 4.0 * (times / 1)**2\n",
    "variable_tone = np.sin(times * omega * 2 * np.pi)\n",
    "axes[0].plot(times, variable_tone, lw=.25)\n",
    "axes[0].set_xlabel('Time (s)')\n",
    "axes[0].set_ylabel('Amplitude (Unitless)')\n",
    "\n",
    "# note: ax.specgram returns 4 arguments, the last of which is the image data needed for the colorbar\n",
    "_, _, _, im = axes[1].specgram(variable_tone, NFFT=NFFT, Fs=sampling_freq, window=window,\n",
    "                               noverlap=noverlap, xextent=[times[0], times[-1]], cmap='inferno')\n",
    "axes[1].set_xlabel('Time (s)')\n",
    "axes[1].set_ylabel('Frequency (Hz)')\n",
    "plt.colorbar(im, ax=axes[1])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectrogram allows us to clearly see the frequency's quadratic time dependence. Try playing around with `spectrogram_interval`, `noverlap`, and the `window` to see how they affect the spectrogram. \n",
    "\n",
    "### Your turn. Plot a spectrogram of the strain at both Hanford and Livingston using `plt.specgram`. \n",
    "**Hint: You should see lots of noise at all frequencies (just as in the power spectrum), as well as horizontal lines corresponding to the spikes in the power spectrum (these are resonances in the suspension system of the LIGO mirrors). Do you see anything interesting at t$\\approx$0 if you zoom in?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:25.709515Z",
     "start_time": "2020-11-24T22:19:25.707577Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5: Looking for the gravitational wave chirp\n",
    "\n",
    "Now let's put together what we've learned about whitening and spectrograms to plot the gravitational wave signal as cleanly as possible.\n",
    "\n",
    "### Use the interpolated powers from Problem 1, the whitening function from Problem 3, and the `specgram` method you employed above to show the strain for both Hanford and Livingston as a function of time and frequency.\n",
    "**Hint: use `ax.set_ylim` and `ax.set_xlim` to zoom in on the signal at t$\\approx$0 and the frequency range where the signal is detectable over the noise.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:37.756862Z",
     "start_time": "2020-11-24T22:19:37.754587Z"
    },
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-20T01:53:52.763643Z",
     "start_time": "2020-11-20T01:53:52.712293Z"
    }
   },
   "source": [
    "You should see a \"chirp\", a signal that's rising in frequency. *That's the signal!* The frequency is set by the orbital period of the two black holes, which emitting energy in the form of gravitational waves. As they lose energy, their orbit decays or \"inspirals,\" increasing the orbital frequency until they merge.\n",
    "\n",
    "(The reason the whitening method worked so well here is because the power spectra we computed in Problem 1 include many time samples with just noise and only a few times with signal in them. So they are good estimates of the background noise-only, including the bright resonant lines from vibrational modes in the detector.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 6: Filtering and plotting the gravitational wave signal.\n",
    "\n",
    "The spectrogram above has the virtue of showing us different frequencies separately. But if we just want to see the waveform itself, our best bet is not just to whiten, but actually filter off the frequencies we know are just noise. Filtering is a deep topic in signal processing, but the absolute simplest frequency filter is simply to Fourier transform, set some frequencies to 0, and then Fourier transform back. This isn't ideal for many applications, but it'll do for our purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:39.205229Z",
     "start_time": "2020-11-24T22:19:39.200428Z"
    }
   },
   "outputs": [],
   "source": [
    "def bandpass_filter(data, times, fmin=0, fmax=np.inf):\n",
    "    '''This simple filter zeros-out fourier content in the data\n",
    "    outside of the the frequency range fmin to fmax.'''\n",
    "    fft_data = np.fft.rfft(data)\n",
    "    freqs = np.fft.rfftfreq(len(data), d=(times[1] - times[0]))\n",
    "    fft_data[np.abs(freqs) < fmin] = 0\n",
    "    fft_data[np.abs(freqs) > fmax] = 0\n",
    "    return np.fft.irfft(fft_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the whitened and filtered gravitational wave strain as a function of time for both Hanford and Livingston, zooming in on the chirp.\n",
    "**Hints: You'll notice that the filter makes the strain go wild at at the beginning and end of the interval. There are lots of ways of combating this effect, like using window functions, but for now it's not necessary since the signal is at t$\\approx$0. Also, don't worry about the units on the y-axis; we could work them out with a more careful whitening, but that's beyond the scope for this exercise.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:43.733812Z",
     "start_time": "2020-11-24T22:19:43.731554Z"
    }
   },
   "outputs": [],
   "source": [
    "# Your code here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see that the two waves are out of sync with each other. This is because Hanford detected the event about 7 ms before Livingnston. This is due to the fact that gravitational waves travel at the speed of light, and Hanford was slightly closer to the merger at the time the wave passed by Earth. With 3 or more detectors, these time delays can be used to triangulate the position of the source on the sky. This was done a couple of years after this with an event called [GW170817](https://en.wikipedia.org/wiki/GW170817), the merger of two neutron stars, which was also observed as a gamma-ray burst and then followed up in the X-ray, optical, infrared, and radio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 7: \"Hearing\" a graviational wave\n",
    "\n",
    "As we said before, the gravitational wave signal is in the audible range, though it's rather low in that range. \n",
    "\n",
    "### Problem 7a: Use whitening and filtering to play the gravitational wave signal at Hanford as a sound.\n",
    "**Hint: just play the middle few seconds of the filtered and whitened signal to avoid edge effects.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:49.144758Z",
     "start_time": "2020-11-24T22:19:49.142636Z"
    }
   },
   "outputs": [],
   "source": [
    "# your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should hear a faint thump. To make the iconic \"chirp\" that that was shared in widely in the media when the detection was announced in 2016, this signal was artificially frequency shifted up by 400 Hz to sound like a clearer blip to the human ear. You can do that yourself with the following function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:50.432443Z",
     "start_time": "2020-11-24T22:19:50.427090Z"
    }
   },
   "outputs": [],
   "source": [
    "def freq_shift(data, times, fshift=0):\n",
    "    \"\"\"Shift up all frequency components of the signal by fshift.\"\"\"\n",
    "    # calculate frequency resolution\n",
    "    df = 1.0 / (len(data) * (times[1] - times[0]))\n",
    "    # calculate shift amount\n",
    "    shift = int(fshift / df)\n",
    "    # Fourier transform\n",
    "    fft_data = np.fft.rfft(data)\n",
    "    # Shift\n",
    "    shifted_fft_data = np.roll(fft_data.real, shift) + 1.0j * np.roll(fft_data.imag, shift)\n",
    "    # Remove low frequency content (which would have shifted over from high frequency using np.roll)\n",
    "    shifted_fft_data[0:shift] = 0\n",
    "    # Fourier transform back\n",
    "    return np.fft.irfft(shifted_fft_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 7b: Make the famous \"chirp\" by upshifting the Hanford strain frequency by 400 Hz and playing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-24T22:19:54.961030Z",
     "start_time": "2020-11-24T22:19:54.959048Z"
    }
   },
   "outputs": [],
   "source": [
    "# Your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should sound a lot like what you hear in [this news story](https://www.npr.org/sections/thetwo-way/2016/02/11/466286219/in-milestone-scientists-detect-waves-in-space-time-as-black-holes-collide) from the day of the announcement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Acknowledgements:\n",
    "\n",
    "### This work was supported by the National Science Foundation Astronomy and Astrophysics Postdoctoral Fellowship (Award \\#1701536).\n",
    "\n",
    "![](data/nsf.png)"
   ]
  },
  {
   "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.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "310.594px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}