WSJT-X and linear satellites: part I

wsjtx-doppler.py

#!/usr/bin/env python3

import ephem

import numpy as np
import matplotlib.pyplot as plt

from scipy.io import wavfile
from scipy.signal import hilbert

import tempfile
import subprocess
import os

from enum import Enum

WSJTX_PATH = '/home/daniel/wsjt/branches/wsjtx/build/'

class JTMode(Enum):
    jt9a = 1
    qra64a = 2
    qra64c = 3
    qra64e = 4
    ft8 = 5
    jt4g = 6

samp_rate = 12000
    
c = 299792458
f_down = 436.5e6

line1 = "LILACSAT-1"
line2 = "1 42725U 98067ME  17188.90653574  .00007620  00000-0  11662-3 0  9998"
line3 = "2 42725  51.6409 294.7316 0007224  26.7821 333.3542 15.55427837  6725"

sat = ephem.readtle(line1, line2, line3)

qth = ephem.Observer()
qth.lat, qth.lon, qth.elevation = '40.6007', '-3.7080', 700

start = ephem.Date('2017/7/9 03:12:00')
t_step = 1e-3
ts = np.arange(0, 600, t_step)

def jt9sim(msg, span, nsigs, snr):
    subprocess.call([WSJTX_PATH + 'jt9sim', msg, str(span), str(nsigs), '1', str(snr), '1'],
                    stdout = subprocess.DEVNULL)
    filename = '000000_0000.wav'
    _, wav = wavfile.read(filename)
    os.remove(filename)
    return wav

def qra64sim(msg, submode, snr):
    subprocess.call([WSJTX_PATH + 'qra64sim', msg, submode,
                     '1', '0', '0', '1', str(snr)],
                    stdout = subprocess.DEVNULL)
    filename = '000000_0001.wav'
    _, wav = wavfile.read(filename)
    os.remove(filename)
    return wav

def ft8sim(msg, snr):
    subprocess.call([WSJTX_PATH + 'ft8sim', msg, 's', '1500.0', '0.0', '0.0', '0.0', '1', str(snr)],
                    stdout = subprocess.DEVNULL)
    filename = '000000_000001.wav'
    _, wav = wavfile.read(filename)
    os.remove(filename)
    return wav

def jt4sim(msg, submode, snr):
    subprocess.call([WSJTX_PATH + 'jt4sim', msg, submode,
                     '1', '0', '0', '1', str(snr)],
                    stdout = subprocess.DEVNULL)
    filename = '000000_0001.wav'
    _, wav = wavfile.read(filename)
    os.remove(filename)
    return wav

def jt9(mode, f):
    if mode is JTMode.jt9a:
        args = ['-9', '-d', '3']
    elif mode is JTMode.qra64a:
        args = ['-q', '-f', '1000']
    elif mode is JTMode.qra64c:
        args = ['-q', '-b', 'C', '-f', '1000']
    elif mode is JTMode.qra64e:
        args = ['-q', '-b', 'E', '-f', '1000']
    elif mode is JTMode.ft8:
        args = ['-8', '-d', '3']
    elif mode is JTMode.jt4g:
        args = ['-4', '-b', 'G', '-f', '1000']
    subprocess.call([WSJTX_PATH + 'jt9'] + args + [f])

def jtGenerate(mode, snr = None):
    message = 'EA4GPZ M0HXM IO94'
    if mode is JTMode.jt9a:
        if snr == None:
            snr = -24
        return jt9sim(message, 1000, 25, -24)
    if mode is JTMode.qra64a:
        if snr == None:
            snr = -24
        return qra64sim(message, 'A', snr)
    if mode is JTMode.qra64c:
        if snr == None:
            snr = -24
        return qra64sim(message, 'C', snr)
    if mode is JTMode.qra64e:
        if snr == None:
            snr = -20
        return qra64sim(message, 'E', snr)
    if mode is JTMode.ft8:
        if snr == None:
            snr = -18
        return ft8sim(message, snr)
    if mode is JTMode.jt4g:
        if snr == None:
            snr = -20
        return jt4sim(message, 'G', snr)

def trials(mode):
    if mode is JTMode.jt9a:
        return 1
    if mode is JTMode.qra64a:
        return 10
    if mode is JTMode.qra64c:
        return 10
    if mode is JTMode.qra64e:
        return 10
    if mode is JTMode.ft8:
        return 10
    if mode is JTMode.jt4g:
        return 10
    
def jtDecode(mode, wav, doppler, signal_start):
    # Generate complex oscillator for mixing with Doppler shift
    phase = 0.0
    drift_oscillator = np.empty_like(wav, dtype=np.complex)
    for j in range(len(drift_oscillator)):
        index = int((signal_start + j/samp_rate)/t_step)
        freq = doppler[index]
        phase += 2*np.pi*freq/samp_rate
        phase = np.fmod(phase + np.pi, 2*np.pi) - np.pi
        drift_oscillator[j] = np.exp(1j*phase)

    # Mix with Doppler oscillator
    wav_doppler = np.real(drift_oscillator * hilbert(wav))

    f_out = tempfile.NamedTemporaryFile()
    wavfile.write(f_out, samp_rate, wav_doppler.astype(np.int16))

    jt9(mode, f_out.name)
    f_out.close()
    return

def computeDoppler(tles, delay = 0.0, freq = f_down):
    dopplers = np.empty_like(ts)
    for j in range(len(ts)):
        qth.date = start + ts[j]*ephem.second - delay*ephem.second
        tles.compute(qth)
        dopplers[j] = -tles.range_velocity/c*freq
    return dopplers

def dopplerDelays(width, step, signal_start = None, mode = None, snr = None, plot = False):
    plt.figure(1)
    doppler = computeDoppler(sat)
    for delay_param in np.arange(-width/2, width/2, step):
        print('Calculating Doppler with delay', delay_param)
        if not plot:
            wavs = [jtGenerate(mode, snr) for _ in range(trials(mode))]
        residual_doppler = doppler - computeDoppler(sat, delay_param)
        if plot:
            plt.plot(ts, residual_doppler)
        else:
            [jtDecode(mode, wav, residual_doppler, signal_start) for wav in wavs]
    return


plt.figure(1)
plt.plot(ts, computeDoppler(sat))
plt.xlabel('Time (s)')
plt.ylabel('Doppler (Hz)')
plt.title('Downlink Doppler at 436.5MHz starting on 2017/7/9 03:12:00')

plt.figure(2)
dopplerDelays(0.4, 0.01, plot = True)
plt.xlabel('Time (s)')
plt.ylabel('Residual Doppler (Hz)')
plt.title('Residual Dopplers for offsets from -0.2s to 0.19s in 0.01s increments')

plt.figure(3)
dopplerDelays(0.8, 0.02, plot = True) # OK with 0.05 @ -24dB
plt.xlabel('Time (s)')
plt.ylabel('Residual Doppler (Hz)')
plt.title('Residual Dopplers for offsets from -0.4s to 0.38s in 0.02s increments')

print('Low SNR')
dopplerDelays(0.2, 0.01, 320, JTMode.jt9a, -24) # OK with 0.05
dopplerDelays(0.4, 0.01, 320, JTMode.qra64a, -24) # OK with 0.12
dopplerDelays(0.4, 0.01, 320, JTMode.qra64c, -24) # OK with 0.13
dopplerDelays(0.4, 0.01, 320, JTMode.qra64e, -22) # OK with 0.16
dopplerDelays(0.4, 0.01, 320, JTMode.ft8, -18) # OK with 0.17
dopplerDelays(0.4, 0.01, 320, JTMode.jt4g, -22) # OK with 0.12

print('-20dB SNR')
dopplerDelays(0.2, 0.01, 320, JTMode.jt9a, -20) # OK with 0.06
dopplerDelays(0.4, 0.01, 320, JTMode.qra64a, -20) # OK with 0.12
dopplerDelays(0.4, 0.01, 320, JTMode.qra64c, -20) # OK with 0.16
dopplerDelays(0.4, 0.01, 320, JTMode.qra64e, -20) # OK with 0.18
dopplerDelays(0.4, 0.01, 320, JTMode.ft8, -20) # xxxx
dopplerDelays(0.4, 0.01, 320, JTMode.jt4g, -20) # OK with 0.14

print('-15dB SNR')
dopplerDelays(0.2, 0.01, 320, JTMode.jt9a, -15) # OK with 0.06
dopplerDelays(0.4, 0.01, 320, JTMode.qra64a, -15) # OK with 0.12
dopplerDelays(0.4, 0.01, 320, JTMode.qra64c, -15) # OK with 0.17
dopplerDelays(0.5, 0.01, 320, JTMode.qra64e, -15) # OK with 0.21
dopplerDelays(0.6, 0.01, 320, JTMode.ft8, -15) # OK with 0.26
dopplerDelays(0.4, 0.01, 320, JTMode.jt4g, -15) # OK with 0.16

plt.show()