FO-29 FT8 waterfall.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import ephem\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.signal import blackman\n",
    "\n",
    "from multiprocessing import Pool\n",
    "\n",
    "import time\n",
    "\n",
    "# Larger figure size\n",
    "fig_size = [14, 8]\n",
    "plt.rcParams['figure.figsize'] = fig_size\n",
    "\n",
    "c = 299792458\n",
    "\n",
    "line1 = \"FO-29\"\n",
    "line2 = \"1 24278U 96046B   17296.90226955 -.00000046 +00000-0 -99197-5 0  9993\"\n",
    "line3 = \"2 24278 098.5410 162.9058 0350559 003.9820 356.3983 13.53079989046036\"\n",
    "\n",
    "sat = ephem.readtle(line1, line2, line3)\n",
    "\n",
    "ea4gpz_qth = ephem.Observer()\n",
    "ea4gpz_qth.lat, ea4gpz_qth.lon, ea4gpz_qth.elevation = '40.5962', '-3.6963', 700\n",
    "\n",
    "start = ephem.Date('2017/10/23 20:26:00')\n",
    "t_step = 1e-2\n",
    "ts = np.arange(0, 990, t_step)\n",
    "\n",
    "f_down = 435.850e6\n",
    "\n",
    "def computeDoppler(tles, freq, delay = 0.0, qth = ea4gpz_qth):\n",
    "    dopplers = np.empty_like(ts)\n",
    "    for j in range(len(ts)):\n",
    "        qth.date = start + ts[j]*ephem.second - delay*ephem.second\n",
    "        tles.compute(qth)\n",
    "        dopplers[j] =  -tles.range_velocity/c*freq\n",
    "    return dopplers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "downlink_doppler = computeDoppler(sat, f_down)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wav_iq = np.memmap('/home/daniel/fo29ft8.wav', dtype=np.int16, mode='r', offset=0x28)\n",
    "samp_rate = 192000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def oscillator(doppler, length, signal_start = 0):\n",
    "    phase = 0.0\n",
    "    x = np.empty(length, dtype=np.complex)\n",
    "    for j in range(length):\n",
    "        index = int((signal_start + j/samp_rate)/t_step)\n",
    "        freq = -doppler[index]\n",
    "        phase += 2*np.pi*freq/samp_rate\n",
    "        phase = np.fmod(phase + np.pi, 2*np.pi) - np.pi\n",
    "        x[j] = np.exp(1j*phase)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 4096\n",
    "#N = 1024\n",
    "w = blackman(N)\n",
    "\n",
    "def process_fft(chunk):\n",
    "    piece = wav_iq[chunk * N : (chunk + 2) * N]\n",
    "    x = piece[::2] + 1j * piece[1::2]\n",
    "    x = x - np.average(x) # remove DC\n",
    "    o = oscillator(downlink_doppler, length = N, signal_start = chunk * N//2 / samp_rate)\n",
    "    f = np.fft.fftshift(np.fft.fft(x * o * w))\n",
    "    return np.real(f * np.conj(f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_chunks = len(wav_iq) // N - 1\n",
    "lines = 4096\n",
    "#lines = 990\n",
    "averaging = total_chunks // lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def process_line(line):\n",
    "    chunks = range(line*averaging, (line+1)*averaging)\n",
    "    return 10*np.log10(np.average(np.array([process_fft(chunk) for chunk in chunks]), axis = 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1709.3412630781531\n"
     ]
    }
   ],
   "source": [
    "t1 = time.perf_counter()\n",
    "with Pool() as p:\n",
    "        waterfall = np.array(p.map(process_line, range(lines)))\n",
    "t2 = time.perf_counter()\n",
    "print(t2 - t1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.imsave('fo29-waterfall.png', waterfall, cmap = 'viridis', vmin=52, vmax = 80)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python3.4",
   "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.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}