On 2018-02-23, Wei Mingchuan BG2BHC published a time synchronized recording of LilacSat-2 made from groundstations at Chongqing and Harbin (baseline of 2500km).
The recording and Wei's MATLAB files to process it can be found in amateur-vlbi.
Here I process the recording using some techniques from GNSS signal processing.
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
import scipy.io
import ephem
c = 299792458 # speed of light in m/s
The recording is stored in a MATLAB file. The path below should be adjusted to the location of the file.
In [2]:
matfile = scipy.io.loadmat('/home/daniel/amateur-vlbi/lilacsat2_sync_20180223_slice.mat')
The contents of the MATLAB file are as follows:
In [3]:
matfile
Out[3]:
We get the sample rate fs and the two channel IQ signal from the MATLAB file.
In [4]:
fs = matfile['rx_rate'][0][0]
signal_chongqing = matfile['signal_sync_chongqing_slice'][0,:]
signal_harbin = matfile['signal_sync_harbin_slice'][0,:]
signal = np.vstack((signal_chongqing, signal_harbin))
The time difference between both recordings is on the order of 1us (or 300m), much smaller that the sample period, so we can ignore it for now.
In [5]:
clock_difference = (matfile['rx_time_chongqing'] - matfile['rx_time_harbin'])[0][0]
clock_difference
Out[5]:
Now we plot the spectrum of both recordings. They show several 9k6 GMSK packets transmitted by LilacSat-2 at 437.225MHz.
Note that the frequency of the signal is seen higher in Chongqing than in Harbin due to the Doppler effect. LilacSat-2 is moving away from Harbin and towards Chongqing, since this is a north to south pass over China.
In [6]:
def spectrum_plot(x, title = None):
f, t, Sxx = scipy.signal.spectrogram(x, fs, return_onesided=False)
f = np.fft.fftshift(f)
plt.imshow(np.fft.fftshift(np.log10(Sxx), axes=0), extent = [t[0],t[-1],f[-1],f[0]], aspect='auto', cmap='viridis')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
if title:
plt.title(title)
In [7]:
spectrum_plot(signal[0,:], 'Chongqing recording')
In [8]:
spectrum_plot(signal[1,:], 'Harbin recording')
Now we will convert both signals to baseband. The frequency seen at Chongqing is approximately -44kHz and the frequency seen at Harbin is approximately -60.5kHz. We use these frequencies to generate a "local oscillator" (or NCO) that brings both signals near 0Hz.
In [9]:
f_lo_chongqing = -44000
f_lo_harbin = -60500
f_lo = np.vstack((f_lo_chongqing, f_lo_harbin))
lo = np.exp(-1j*2*np.pi*np.arange(signal.shape[1])*f_lo/fs)
We will use a FIR filter to filter the signal. The filter is designed for a cutoff at a normalized frequency of 0.1, or 25kHz.
In [10]:
lowpass = scipy.signal.firwin(128, 0.1)
freqz = scipy.signal.freqz(lowpass)
plt.plot(freqz[0]/np.pi, 20*np.log10(np.abs(freqz[1])))
plt.title('Lowpass filter frequency response')
plt.xlabel('Normalized frequency')
plt.ylabel('Gain (dB)');
Now we apply the LO and filter to the signal.
In [11]:
signal_bb = scipy.signal.lfilter(lowpass, 1, signal*lo)
The results of frequency translation and filtering are as follows:
In [12]:
spectrum_plot(signal_bb[0,:], 'Chongqing recording (baseband)')
In [13]:
spectrum_plot(signal_bb[1,:], 'Harbin recording (baseband)')
The following function try to estimate the position of a correlation peak, given "prompt", "early" and "late" correlators. The idea is taken from GNSS signal processing. The functions can be used as small correction (smaller than 0.5 samples) to obtain a resolution for the delay better than one sample when doing the correlation.
There are two functions depending on the shape of the correlation peak. The triangle function assumes that the correlation peak is an isosceles triangle (which is the case for a rectangular pulse shaping filter). The parabola function assumes that the correlation peak is a parabola (which is an approximation of what happens for a continuous pulse shaping filter). Here we use the parabola, as the correlation peaks are well matched by a parabola.
In [215]:
def peak_estimator_triangle(x):
return np.clip(0.5*(x[2]-x[0])/(x[1]-min(x[0],x[2])), -0.5, 0.5)
def peak_estimator_parabola(x):
return np.clip(0.5*(x[2]-x[0])/(2*x[1]-x[0]-x[2]), -0.5, 0.5)
peak_estimator = peak_estimator_parabola
This is the main open-loop correlation algorithm. The algorithm works by using an FFT of size N in contiguous non-overlapping blocks.
Correlation is done writing the convolution in time domain as a product in the frequency domain. Search in Doppler is done by circular shifting one of the signals in the frequency domain, which amounts multiplying by a complex exponential. Search in Doppler only provides a very coarse frequency estimate and it is only performed here to maximize the correlation peak, since the correlation has a sinc()-like response in frequency. Fine frequency will be calculated later using phase techniques.
For each block of N samples, the Doppler bin producing the largest correlation is found, and the delay where the correlation peak occurs is refined using the peak_estimator() function above. We also store the complex correlation to use it later for phase computations.
In [216]:
N = 2**14
hz_per_bin = fs/N
freq_search_hz = 200
bin_search = int(freq_search_hz/hz_per_bin)
steps = signal_bb.shape[1]//N
best_freq_bin = np.empty(steps)
best_time_bin = np.empty(steps)
complex_corr = np.empty(steps, dtype='complex64')
for step in range(steps):
f = np.fft.fft(signal_bb[:,step*N:(step+1)*N])
best_corr = -np.inf
for freq_bin in range(-bin_search,bin_search):
corr = np.fft.ifft(f[0,:]*np.conj(np.roll(f[1,:], freq_bin)))
corr_max = np.max(np.abs(corr))
if corr_max > best_corr:
best_corr = corr_max
best_freq_bin[step] = freq_bin
j = np.argmax(np.abs(corr))
best_time_bin[step] = j + peak_estimator(np.abs(corr[j-1:j+2]))
complex_corr[step] = corr[j]
The vector t below represents the times of measurement of each correlation. t = 0 corresponds to the start of the recording. The time for each correlation is taken as the time for the midpoint of the correlation interval (hence the 0.5 below).
In [217]:
t = (0.5 + np.arange(steps))*N/fs
Since the signal is continuous but packetised, a meaningful correlation is only produced during packet transmissions. Also, the beginning and end of each packet contain a periodic preamble which causes many self-correlation peaks and confuses the cross-correlation. Therefore, the correlation is only valid when the data is transmitted. Since the data is scrambled, it provides a good cross-correlation with a single peak.
We can plot the magnitue of the correlation to see clearly when the packets are transmitted.
In [218]:
correlating = np.abs(complex_corr) > 1e-3
plt.plot(t, np.abs(complex_corr))
plt.ylabel('Correlation')
plt.xlabel('Time (s)')
plt.title('Correlation magnitude');
The correlation delay can be converted directly to a physical magnitude: delta-range in meters. When doing so, we also incorporate the clock offset between receivers.
In [219]:
delta_range = (best_time_bin / fs + clock_difference)*c
plt.plot(t, delta_range)
plt.ylim([200e3, 600e3])
plt.xlabel('Time (s)')
plt.ylabel('Delta-range (m)')
plt.title('Delta-range');