In [1]:

    

%matplotlib inline

import ephem

import numpy as np
import matplotlib.pyplot as plt

from scipy.signal import blackman

from multiprocessing import Pool

import time

# Larger figure size
fig_size = [14, 8]
plt.rcParams['figure.figsize'] = fig_size

c = 299792458

line1 = "FO-29"
line2 = "1 24278U 96046B   17296.90226955 -.00000046 +00000-0 -99197-5 0  9993"
line3 = "2 24278 098.5410 162.9058 0350559 003.9820 356.3983 13.53079989046036"

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

ea4gpz_qth = ephem.Observer()
ea4gpz_qth.lat, ea4gpz_qth.lon, ea4gpz_qth.elevation = '40.5962', '-3.6963', 700

start = ephem.Date('2017/10/23 20:26:00')
t_step = 1e-2
ts = np.arange(0, 990, t_step)

f_down = 435.850e6

def computeDoppler(tles, freq, delay = 0.0, qth = ea4gpz_qth):
    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


    

In [2]:

    

downlink_doppler = computeDoppler(sat, f_down)


    

In [3]:

    

wav_iq = np.memmap('/home/daniel/fo29ft8.wav', dtype=np.int16, mode='r', offset=0x28)
samp_rate = 192000


    

In [4]:

    

def oscillator(doppler, length, signal_start = 0):
    phase = 0.0
    x = np.empty(length, dtype=np.complex)
    for j in range(length):
        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
        x[j] = np.exp(1j*phase)
    return x


    

In [10]:

    

N = 4096
#N = 1024
w = blackman(N)

def process_fft(chunk):
    piece = wav_iq[chunk * N : (chunk + 2) * N]
    x = piece[::2] + 1j * piece[1::2]
    x = x - np.average(x) # remove DC
    o = oscillator(downlink_doppler, length = N, signal_start = chunk * N//2 / samp_rate)
    f = np.fft.fftshift(np.fft.fft(x * o * w))
    return np.real(f * np.conj(f))


    

In [11]:

    

total_chunks = len(wav_iq) // N - 1
lines = 4096
#lines = 990
averaging = total_chunks // lines


    

In [12]:

    

def process_line(line):
    chunks = range(line*averaging, (line+1)*averaging)
    return 10*np.log10(np.average(np.array([process_fft(chunk) for chunk in chunks]), axis = 0))


    

In [13]:

    

t1 = time.perf_counter()
with Pool() as p:
        waterfall = np.array(p.map(process_line, range(lines)))
t2 = time.perf_counter()
print(t2 - t1)


    

    




1709.3412630781531

In [16]:

    

plt.imsave('fo29-waterfall.png', waterfall, cmap = 'viridis', vmin=52, vmax = 80)