![](\"BHmerger.png\")
\n",
"![](\"LIGO_detector.png\")
\n",
"![](\"fft1.png\")
\n",
"and the inverse transform is defined as follows\n",
"![](\"fft2.png\")
\n",
"![](\"fft3.png\")
\n",
"\n", "Following https://arxiv.org/pdf/1710.04635.pdf, we model the gravitational wave signal by making analogies to electromagnetic waves.\n", "\n", "Assume that two masses orbit their common center of mass. The orbiting system emits gravitational waves, which carry energy away from the system, leading to a decay of the orbits, that is, there is an “in-spiral” of the masses.\n", "\n", "Starting from an expression that indicates how the orbital frequency changes as the binary system radiates\n", "energy, we can find how the distance between two masses (denoted as $r$) changes with time:\n", "\n", "$$ \\frac{dr}{dt} = \\frac{-\\eta N c}{4}\\Big( \\frac{r_s}{r} \\Big)^3 $$\n", "\n", "where $r_s$ is the Schwarzschild radius, the distance from the center of a spherically symmetric mass at which the escape speed, calculated from Newtonian mechanics in this case, is equal to the speed of light. For two binary masses, $r_s = 2GM/c^2$. $\\eta = m_1 m_2/(m_1 + m_2)^2$ is a dimensionless raio of the masses, and $N$ is a numerical factor which depends on the details of the theory used to derive the result. Here, $\\eta = 0.247, N = 6.4, M = 60.5M_{\\mathrm{sun}}$.\n", "\n", "We can now find the waveform with the following calculational steps:\n", "\n", "1. Let $r_i = 5.3r_s$ (initial separation of two masses) and $t_i$ = time_template[64805] (16.26s since Sep 14 9:50:29 GMT 2015). Assume $t_f = 16.4$. Use $10000$ steps, so $\\Delta t = (t_f-t_i)/10000$\n", "\n", "2. Solve the above differential equation in three different ways:\n", "\n", " (1) Analytically: $r(t) = \\big[r_i^4 - N\\eta r_s^3 c (t-t_i)\\big]^{1/4}$\n", "\n", " (2) Euler's Method (numerically): Let $f(r) = \\frac{dr}{dt}$. Then, $r(t+\\Delta t) = \\Delta t \\cdot f(r)$ for small $\\Delta t$.\n", "\n", " (3) Runge-Kutta Method (numerically): Let $f(r) = \\frac{dr}{dt}$. $k_1 = \\Delta t \\cdot f(r), k_2 = \\Delta t \\cdot f(r+\\frac{1}{2}k_1)$. Then, $r(t+\\Delta t) = r(t) + k_2$ for small $\\Delta t$.\n", "\n", " Yes, it is silly to solve this DEQ numerically since its analytic solution can be found very easily. But if you have nonlinear equations, (2) and (3) can be very useful.\n", "\n", "3. Calculate $r$ for the current time value. ($\\Delta t$ should be very small.)\n", "\n", "4. Use Kepler's third law to find $w(t)$ (which shows orbital frequency) and calculate $w$ for the current time value.\n", "\n", "$$ w^2 = \\frac{GM}{r^3} $$.\n", "\n", "5. Use that value of $w$ and calculate the interferometer strain signal $h(t)$:\n", "\n", "$$ h(t) = \\frac{\\Delta L_x(t) - \\Delta L_y(t)}{L} \\propto w^{2/3}\\mathrm{cos}(2w(t-t_i)) $$\n", "\n", "5. For the next time value, repeat step 3-5 to find $r, w, h$.\n", "\n", "6. Stop when $t$ reaches 16.4 ($t_f = 16.4$).\n", "\n", "7. Scale $h(t)$ so that $h(t_i)$ = template_p[64805]. (Remember that $t_i$ = time_template[64805].) \n", "\n", " 5. Solving the given differential equation in three different ways (analytically, using Euler's method, using Runge-Kutta Method), find three different versions of $h(t)$. Then, show how our estimate of the waveform ($h(t)$) matches with the actual waveform (\"template_p\"); first plot \"template_p\" as a function of \"time_template\" and then plot three different versions of $h(t)$ in the same figure. Try different $M = 55M_{\\mathrm{sun}}, 60.5M_{\\mathrm{sun}}, 65M_{\\mathrm{sun}}$ and show that $M = 60.5M_{\\mathrm{sun}}$ fits best to the given waveform (for $M = 55M_{\\mathrm{sun}}, 65M_{\\mathrm{sun}}$, plot $h(t)$ using only the analytic solution. For $M = 60.5M_{\\mathrm{sun}}$, try all three methods to plot $h(t)$.) Label all plots. set the x limits of the axes as plt.xlim(ti, 16.45).
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "G = 6.67408e-11\n", "Msun = 1.99e30\n", "c = 2.99792e8\n", "\n", "N = 32/5\n", "eta = 0.247\n", "\n", "ti = time_template[64805]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Method 1: Analytic solution\n", "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Method 2: Euler's method\n", "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Method 3: Runge-Kutta Method\n", "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "scrolled": false }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should find that with $M = 60.5M_{\\mathrm{sun}}$, $h(t)$ agrees reasonably well with the actual waveform which uses numerical GR. Our estimate breaks down as we get closer to the merger event. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "From the ASD in Part 4, we can see that noise fluctuations are much larger at low and high frequencies and near spectral lines, reaching a roughly flat (\"white\") minimum in the band around 80 to 300 Hz.\n", "\n", "We can \"whiten\" the data (dividing it by the noise amplitude spectrum, in the fourier domain), suppressing the extra noise at low frequencies and at the spectral lines, to better see the weak signals in the most sensitive band. Note that only the data are needed in this process.\n", "\n", "To get rid of remaining high frequency noise, we will also bandpass the data.\n", "\n", "Now, define a function called \"whiten\" which takes the strain data (\"strain\"), interpolated psd (\"interp_psd\"), and dt (1.0 / fs = 1.0 / 4096).\n", "\n", " Note:\n", "\n", " 1. For L1 and H1, strain = strain_L1 and strain = strain_H1 respectively.\n", "\n", " 2. Let Pxx_H1 be the psd values you obtained in Part 4 for H1. Then, interpolate this as \"psd_H1 = interp1d(freqs, Pxx_H1)\" (https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.interpolate.interp1d.html). psd_H1 is a function; to get the psd values at the needed frequencies (say you have an array \"freqs\"), you should do psd_H1(freqs).\n", "\n", " 3. So \"whiten(strain_H1,psd_H1,1.0/fs)\" would whiten the strain data for H1.\n", "\n", " 6. Define a function \"whiten\". Then, whiten the H1 and L1 strain data. (Replace ellipsis)
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# transform to freq domain, divide by asd, then transform back\n", "def whiten(strain, interp_psd, dt):\n", " # Given strain data, first transform to frequency domain\n", " # Use \"np.fft.rfftfreq\" to get the frequencies and \"np.fft.rfft\" to get the FFT of the strain\n", " freqs = ...\n", " hf = ...\n", " \n", " # Normalization factor\n", " norm = 1./np.sqrt(1./(dt*2))\n", " \n", " # Divide by the ASD and multiply by the normalization factor\n", " white_hf = hf / np.sqrt(interp_psd(freqs)) * norm\n", " \n", " # Transform back to the time domain using \"np.fft.irfft\" (Hint: n = len(strain)) \n", " # https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.irfft.html\n", " white_ht = ...\n", " return white_ht " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "dt = 1.0 / fs\n", "\n", "# Interpolate the PSD values\n", "psd_H1 = interp1d( ... )\n", "psd_L1 = interp1d( ... )\n", "\n", "strain_H1_whiten = whiten(strain_H1,psd_H1,dt)\n", "strain_L1_whiten = whiten(strain_L1,psd_L1,dt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", " Remember that the sampling rate (\"fs\") is 4096 Hz. What is the nyquist frequency?
\n", "(Note: The data cannot capture frequency content above the Nyquist frequency, but that's OK, because GW150914 only has detectable frequency content in the range 20 Hz - 300 Hz.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "nyq = ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, suppress the high frequency noise (no signal there) with some bandpassing. (We bandpass from 43 to 300 Hz, since we are looking for a signal in that frequency band and the detector noise is pretty low there. This will keep power in this frequency range, but suppress power away from this frequency.) We use scipy.signal.butter (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.butter.html) and scipy.signal.filtfilt (https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.signal.filtfilt.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# bandpass filter coefficients \n", "order = 4\n", "lowcut= 43\n", "highcut= 360\n", "low = lowcut / nyq\n", "high = highcut / nyq\n", "\n", "bb, ab = butter(order, [low, high], btype='band')\n", "normalization = np.sqrt((highcut-lowcut)/nyq)\n", "strain_H1_whitenbp = filtfilt(bb, ab, strain_H1_whiten) / normalization\n", "strain_L1_whitenbp = filtfilt(bb, ab, strain_L1_whiten) / normalization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Shift L1 event by 7 ms (Signal arrived 7 ms earlier at L1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "strain_L1_shift = -np.roll(strain_L1_whitenbp,int(0.007*fs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Plot \"strain_L1_shift\" (whitened L1 strain shifted by 7ms) and \"strain_H1_whitenbp\" (whitened H1 strain) as a function of time. Set \"plt.xlim(ti, 16.45); plt.ylim(-10, 10)\". Do you see a signal?
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "The FFT input signal is inherently truncated. This truncation can be modelled as multiplication of an infinite signal with a rectangular window function. In the spectral domain this multiplication becomes convolution of the signal spectrum with the window function spectrum, being of form sin(x)/x. This convolution is the cause of an effect called spectral leakage (https://en.wikipedia.org/wiki/Spectral_leakage)\n", "\n", "Hence, spectral leakage is caused by discontinuities in the original, noninteger number of periods in a signal. Windowing the signal with a dedicated window function helps mitigate spectral leakage by reducing the amplitude of the discontinuities at the boundaries\n", "\n", "
![](\"window.png\")
\n",
"\n",
"In this project, we use the Blackman window (https://docs.scipy.org/doc/numpy/reference/generated/numpy.blackman.html#numpy.blackman). \n",
"\n",
" 7. Recall that in Part 2 we found the FFT of $f(t)$ using np.fft.rfft. Now, define a window function \"w = np.blackman(N)\" and multiply $f(t)$ by w and do the FFT. $\\big($np.fft.rfft(w*f(t))$\\big)$. Plot the FFT of $f(t)$ with and without the window function. Set the y-axis as log-scale (plt.semilogy). Show that windowing reduces spectral leakage. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Matched filtering is the optimal way to find a known signal buried in stationary, Gaussian noise. It is the standard technique used by the gravitational wave community to find GW signals from compact binary mergers in noisy detector data. LIGO scientists use matched filtering to find such \"hidden\" signals. A matched filter works by compressing the entire signal into one time bin (by convention, the \"end time\" of the waveform).\n", "\n", "Here, we use only one template (the one identified in the full search as being a good match to the data). Assuming that the data around this event is fairly Gaussian and stationary, we'll use this simple method to identify the signal (matching the template) in our 32 second stretch of data. \n", "\n", "\n", " 8. Do the matched filtering. First, start with the L1 strain data.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# To calculate the PSD of the data, choose an overlap and a window \n", "fs = event['fs']\n", "dt = 1.0 / fs\n", "NFFT = 4*fs\n", "psd_window = np.blackman(NFFT)\n", "# and a 50% overlap:\n", "NOVL = NFFT/2\n", "\n", "# Complex waveform template\n", "template = (template_p + template_c*1.j) \n", "\n", "# the length and sampling rate of the template MUST match that of the data.\n", "datafreq = np.fft.fftfreq(template.size, d = 1./fs)\n", "df = np.abs(datafreq[1] - datafreq[0])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", " Window the data. (Replace ellipsis in the cell below)
\n", "\n", "1. Define the window: \"np.blackman(template.size)\"\n", "\n", "2. Do the FFT of the template with window using \"np.fft.fft\": multiply \"template\" by \"dwindow\" and then do FFT. The normalization factor is \"1/fs\" (i.e. multiply the FFT of the template by 1/fs)\n", "\n", "3. Do the FFT of the L1 strain data with window using \"np.fft.fft\": multiply the L1 strain data by \"dwindow\" and then do FFT. The normalization factor is \"1/fs\" (i.e. multiply the FFT of the template by 1/fs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# window the data (Replace ellipsis)\n", "\n", "# 1. Define the window\n", "dwindow = ...\n", "\n", "# 2. Do the FFT of the template (with dwindow)\n", "template_fft = ...\n", "\n", "# 3. Take the Fourier Transform (FFT) of the data (with dwindow)\n", "data = strain_L1.copy()\n", "data_fft = ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", " Plot the blackman window as a function of time. Do you think using this window is useful in our analysis? Note that the merger event occurred around 16s since Sep 14 9:50:29 GMT 2015, and the data stretches over 32s.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", " Plot the FFT of the template with and without the window (as a function of frequency). Only use positive frequencies. Plot the squared magnitude of the FFT (power spectrum). Label all plots. Set the range of x-axis between 10 and 2000. Plot both axes in log scale.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need an estimate of the noise power in each FFT bin." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# First, calculate the PSD of the data. Also use an overlap, and window:\n", "data_psd, freqs = mlab.psd(data, Fs = fs, NFFT = NFFT, window=psd_window, noverlap=NOVL)\n", "# Interpolate to get the PSD values at the needed frequencies. This defines the noise power in each frequency bin.\n", "power_vec = np.interp(np.abs(datafreq), freqs, data_psd)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", " Plot the noise power spectral density (defined by \"power_vec\"). Set the range of x-axis between 10 and 2000. Plot both axes in log scale.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Calculate the matched filter output in the time domain. (Replace ellipsis in the cell below)
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
![](\"optimalfitering.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"1. Multiply the Fourier Space template and data, and divide by the noise power in each frequency bin: i.e. Multiply the FFT of the data with the conjugate (using \".conjugate()\") of the FFT of the template and divide by the noise power (\"power_vec\")\n",
"\n",
"2. Taking the Inverse Fourier Transform (IFFT) of the filter output (take it back to the time domain): i.e. Take the inverse FFT (using np.fft.ifft) of the step 1 output and multiply by the factor \"2*fs\"\n",
"\n",
"Hence, the result will be plotted as a function of time off-set between the template and the data:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# 1. Multiply the Fourier Space template and data, and divide by the noise power in each frequency bin.\n",
"optimal = ...\n",
"\n",
"# 2. Taking the Inverse Fourier Transform (IFFT) of the filter output (take it back to the time domain)\n",
"optimal_time = ..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Normalize the matched filter output. (Replace ellipsis in the cell below) \n", "\n", "1. Find $\\sigma = \\sqrt{\\big<\\frac{\\tilde{h}(f)\\tilde{h}^*(f)}{S_n(f)}\\big>}$\n", "(Hint: First multiply the FFT of the template by its conjugate (using \".conjugate()\") and then divide by the noise PSD. Basically you did $\\frac{\\tilde{h}(f)\\tilde{h}^*(f)}{S_n(f)}$. You have an array then (each array element correspond to each element in the frequency array \"datafreq\" you defined earlier). Now, sum array elements and multiply by \"df\" which you defined earlier.\n", "\n", "2. Normalize the matched filter output so that we expect a value of 1 at times of just noise: i.e. Divide the absolute value of \"optimal_time\" by sigma." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 1. Find $\\sigma$\n", "sigmasq = ...\n", "sigma = np.sqrt(sigmasq.real)\n", "\n", "# 2. Normalize the matched filter output so that we expect a value of 1 at times of just noise\n", "SNR = ...\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Apply time offset, phase to template. Finally, whiten the template and band-pass it. (All routine given - nothing to do here, just run the cell.)
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# shift the SNR vector by the template length so that the peak is at the END of the template\n", "peaksample = int(data.size / 2) # location of peak in the template\n", "SNR_complex = optimal_time/sigma\n", "SNR_complex = np.roll(SNR_complex,peaksample)\n", "SNR = abs(SNR_complex)\n", "\n", "# apply time offset, phase to template\n", "# find the time and SNR value at maximum:\n", "indmax = np.argmax(SNR)\n", "timemax = time[indmax]\n", "SNRmax = SNR[indmax]\n", "\n", "d_eff = sigma / SNRmax\n", "# Calculate optimal horizon distnace\n", "horizon = sigma/8\n", "\n", "# Extract time offset and phase at peak\n", "phase = np.angle(SNR_complex[indmax])\n", "offset = (indmax-peaksample)\n", "\n", "# apply time offset, phase, and d_eff to template \n", "template_phaseshifted = np.real(template*np.exp(1j*phase)) # phase shift the template\n", "template_rolled = np.roll(template_phaseshifted,offset) / d_eff # Apply time offset and scale amplitude\n", "\n", "# Whiten and band-pass the template for plotting\n", "template_whitened = whiten(template_rolled,interp1d(freqs, data_psd),dt) # whiten the template\n", "template_match = filtfilt(bb, ab, template_whitened) / normalization # Band-pass the template\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Plot the original, phase-shifted template (\"template_rolled\") and the whitened template (\"template_whitened\"). Show the plot between the time range (15.6, 16.45) since Sep 14 9:50:29 GMT 2015. You will see that the waveform before ~ 16.3s is largely suppressed. That is because we divide the FFT of the template by the noise PSD in the Fourier space during whitening. Now, calculate FFT(template)/Noise PSD$^{1/2}$ and plot it as a function of frequency in the log-log plot. Set the x-axis range between 5 and 2000 Hz.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Plot SNR as a function of time. (Remember that the GW signal arrived 7 ms earlier at L1. So shift the L1 time accordingly - this is done conventionally.) Set plt.xlim(16.3, 16.5).
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "scrolled": false }, "outputs": [], "source": [ "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Shift L1 event by 7 ms (Signal arrived 7 ms earlier at L1)\n", "template_match_shift_L1 = -np.roll(template_match,int(0.007*fs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " At what time does SNR peak? Print the result. (Remember that the GW signal arrived 7 ms earlier at L1. So shift the L1 time accordingly - this is done conventionally.)
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " First, plot the whitened strain data (from Part 6) as a function of time. Then, plot the matched waveform (\"template_match_shift_L1\") in the same figure. Show that the whitened strain data match well to the waveform predicted by GR.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " 9. Repeat Part 8 for the H1 strain data and do: 1. Plot SNR vs. time, 2. plot the whitened strain data and matched waveform vs. time, 3. find the time at which SNR peaked.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " 10. First, compare the SNR peak time for both H1 and L1 data. (Print the results) Then, plot the whitened L1 and H1 strain in the same plot. Plot the waveform matched to the H1 strain data in the same plot as well.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we conclude that the whitened strain data match well to the waveform predicted by General Relativity; we have detected the gravitational wave predicted by GR!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Extra (all the routine given; just run the cell!):\n", "\n", "Make wav (sound) files from the filtered, downsampled data, +-2s around the event." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "eventname = 'GW150914' \n", "\n", "Pxx = (1.e-22*(18./(0.1+freqs))**2)**2+0.7e-23**2+((freqs/2000.)*4.e-23)**2\n", "psd_smooth = interp1d(freqs, Pxx)\n", "\n", "template_offset = 16.\n", "\n", "# function to keep the data within integer limits, and write to wavfile:\n", "def write_wavfile(filename,fs,data):\n", " d = np.int16(data/np.max(np.abs(data)) * 32767 * 0.9)\n", " wavfile.write(filename,int(fs), d)\n", "\n", "deltat_sound = 2. # seconds around the event\n", "\n", "# index into the strain time series for this time interval:\n", "indxd = np.where((time >= tevent-deltat_sound) & (time < tevent+deltat_sound))\n", "\n", "# write the files:\n", "write_wavfile(eventname+\"_H1_whitenbp.wav\",int(fs), strain_H1_whitenbp[indxd])\n", "write_wavfile(eventname+\"_L1_whitenbp.wav\",int(fs), strain_L1_whitenbp[indxd])\n", "\n", "# re-whiten the template using the smoothed PSD; it sounds better!\n", "template_p_smooth = whiten(template_p,psd_smooth,dt)\n", "\n", "# and the template, sooming in on [-3,+1] seconds around the merger:\n", "indxt = np.where((time >= (time[0]+template_offset-deltat_sound)) & (time < (time[0]+template_offset+deltat_sound)))\n", "write_wavfile(eventname+\"_template_whiten.wav\",int(fs), template_p_smooth[indxt])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With good headphones, you may be able to hear a faint thump in the middle; that's our signal!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "fna = eventname+\"_template_whiten.wav\"\n", "print(fna)\n", "Audio(fna)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "fna = eventname+\"_H1_whitenbp.wav\"\n", "print(fna)\n", "Audio(fna)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## To Submit\n", "Execute the following cell to submit.\n", "If you make changes, execute the cell again to resubmit the final copy of the notebook, they do not get updated automatically.
\n", "__We recommend that all the above cells should be executed (their output visible) in the notebook at the time of submission.__
\n", "Only the final submission before the deadline will be graded. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "_ = ok.submit()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }