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

```
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))

```
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')

```
```

```
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)

```
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]

`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');

```
```

```
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');

```
```