Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
WSJT-X and linear satellites: part I
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment