voyager/voyager_notebook_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Voyager 2\n",
    "\n",
    "Example data taken on 2018-10-22 during MARS receiver testing, using the Breakthrough Listen backend.\n",
    "\n",
    "Data recorded over full bandwidth of MARS receiver, here we have extracted a small bandwidth around the Voyager2 telemetry signal.\n",
    "\n",
    "Spectral data product has 2.79 Hz resolution, 18.25 s time integrations.\n",
    "\n",
    "*Change `datadir` to match your local config.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = \"/scratch/bluse/PKS\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import blimpy as bl\n",
    "import pylab as plt\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams['font.size'] = 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamic spectra data\n",
    "\n",
    "Data is stored in filterbank format, which can be loaded using [blimpy](https://github.com/ucberkleyseti/blimpy):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = os.path.join(datadir, 'voyager2_hires_2018.10.22.h5')\n",
    "a = bl.Waterfall(filename)\n",
    "a.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting and showing dynamic spectra, we can see the signal drifting due to doppler acceleration (LO does not correct for LSRK)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "a.plot_spectrum()\n",
    "\n",
    "plt.figure(figsize=(14, 6))\n",
    "a.plot_waterfall()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract sidebands and zoom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc_cw   = 8575.34545  # Center freq of transmission\n",
    "sb_sep  = 0.0225      # Sideband separation from carrier in MHz\n",
    "sb_bw   = 0.001       # Sideband bandwidth\n",
    "plot_bw = sb_bw / 2   # Plotting bandwidth\n",
    "t = a.timestamps\n",
    "\n",
    "# Calculate center freq of sidebands\n",
    "fc_lsb = fc_cw - sb_sep\n",
    "fc_usb = fc_cw + sb_sep\n",
    "\n",
    "# Extract data\n",
    "f_cw,   d_cw   = a.grab_data(f_start=fc_cw-plot_bw, f_stop=fc_cw+plot_bw)\n",
    "f_lsb,  d_lsb  = a.grab_data(f_start=fc_lsb-plot_bw, f_stop=fc_lsb+plot_bw)\n",
    "f_usb,  d_usb  = a.grab_data(f_start=fc_usb-plot_bw, f_stop=fc_usb+plot_bw)\n",
    "\n",
    "# Plot\n",
    "def plot_waterfall(t, f, d):\n",
    "    t_elapsed = (t[-1] - t[0]) * 86400\n",
    "    plt.imshow(d[::-1], aspect='auto', extent=(f[-1], f[0], 0, t_elapsed))\n",
    "    plt.xticks(rotation=30)\n",
    "    plt.xlabel(\"Frequency [MHz]\")\n",
    "    \n",
    "plt.figure(figsize=(14, 5))\n",
    "plt.subplot(1,3,1)\n",
    "plt.ylabel(\"Elapsed time [s]\")\n",
    "plot_waterfall(t, f_lsb, d_lsb)\n",
    "plt.subplot(1,3,2)\n",
    "plot_waterfall(t, f_cw, d_cw)\n",
    "plt.subplot(1,3,3)\n",
    "plot_waterfall(t, f_usb, d_usb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Raw voltage data\n",
    "\n",
    "Raw voltage data for a 2.92 MHz subband was recorded about the Voyager signal. It is stored in Guppi RAW format, as 8-bit integers. This can also be loaded using blimpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_filename = os.path.join(datadir, 'voyager_2018.10.22.raw')\n",
    "raw = bl.GuppiRaw(raw_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_header0, raw_data0 = raw.read_next_data_block()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data shape is (1, 524288, 2), and dtype is presented as complex64 (but stored as 8-bit in file)\n",
    "\n",
    "Axes are (channel, time sample, polarization). You can iterate over multiple blocks of data to get more time samples if required, until the end of the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_data0.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the spectrum from these data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = raw_data0.squeeze()[:, 0]\n",
    "dy = raw_data0.squeeze()[:, 1]\n",
    "\n",
    "n_chans = 8192\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(1,2,1)\n",
    "Dx = np.abs(np.fft.fftshift(np.fft.fft(dx)))\n",
    "Dx = Dx.reshape([n_chans, -1]).mean(axis=-1)\n",
    "plt.plot(Dx)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "Dy = np.abs(np.fft.fftshift(np.fft.fft(dy)))\n",
    "Dy = Dy.reshape([n_chans, -1]).mean(axis=-1)\n",
    "plt.plot(Dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use multiple blocks to integrate the spectrum down:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_chans = 8192\n",
    "n_ints  = 64\n",
    "\n",
    "raw.reset_index()\n",
    "\n",
    "x, y = np.zeros(n_chans), np.zeros(n_chans)\n",
    "for idx in range(n_ints):\n",
    "    h, d = raw.read_next_data_block()\n",
    "    dx = d.squeeze()[:, 0]\n",
    "    dy = d.squeeze()[:, 1]\n",
    "    Dx = np.abs(np.fft.fftshift(np.fft.fft(dx)))\n",
    "    Dx = Dx.reshape([n_chans, -1]).mean(axis=-1)\n",
    "    Dy = np.abs(np.fft.fftshift(np.fft.fft(dy)))\n",
    "    Dy = Dy.reshape([n_chans, -1]).mean(axis=-1)\n",
    "    x += Dx\n",
    "    y += Dy\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(x)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "BLUSE Python 3.8",
   "language": "python",
   "name": "bluse_python38"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}