Created
August 12, 2017 14:19
-
-
Save daniestevez/1eeb6270d0fbe7070f761d351573aeb0 to your computer and use it in GitHub Desktop.
Revisions
-
daniestevez created this gist
Aug 12, 2017 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,221 @@ #!/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()