```
In [1]:
```%matplotlib inline
import numpy as np
import scipy.signal
import matplotlib.pyplot as plt

```
In [2]:
```x = np.fromfile('/home/daniel/Descargas/DSLWP-B_PI9CAM_2018-08-12T08_32_25_436.4MHz_40ksps_complex.raw', dtype='complex64')

The 250bps GMSK signal is converted down to baseband and lowpass filtered to 800Hz.

```
In [3]:
```fs = 40e3
f = 1400
x = x * np.exp(-1j*2*np.pi*np.arange(x.size)*f/fs).astype('complex64')
h = scipy.signal.firwin(1000, 0.01).astype('float32')
x = scipy.signal.lfilter(h, 1, x).astype('complex64')

Perform arctangent FSK demodulation.

```
In [4]:
```s = np.diff(np.angle(x).astype('float32'))

Correct for phase wrapping.

```
In [5]:
```s[s > np.pi] -= 2*np.pi
s[s < -np.pi] += 2*np.pi

```
In [6]:
``````
#s.tofile('/home/daniel/s.f32')
```

```
In [7]:
``````
#s = np.fromfile('/home/daniel/s.f32', dtype='float32')
```

We extract the soft bits by guessing the correct clock phase and decimating. No pulse shaping matched filtering has been done, and tracking of the clock frequency doesn't seem necessary either. The separation between the bits 1 and 0 is good enough for demodulation without bit errors.

Note that we correct for the clock drift and for frequency offset and drift. This simple open loop clock recovery is enough to get good separation between the bits, even across different GMSK transmissions.

```
In [8]:
```phase = 110
softbits = s[np.int32(np.arange(phase,s.size,160 - 0.0006))]
softbits = softbits + 1.7e-8 * np.arange(softbits.size) - 0.0035 # correction for frequency offset and drift
plt.plot(softbits,'.')
plt.axhline(y = 0, color='green')
plt.ylim([-0.03,0.03]);

```
```

```
In [9]:
```starts = np.arange(0,softbits.size,100000)
ranges = zip(starts[:-1], starts[1:])
for r in ranges:
plt.figure()
plt.plot(np.arange(r[0],r[1]), softbits[r[0]:r[1]], '.')
plt.axhline(y = 0, color='green')
plt.ylim([-0.03,0.03]);

```
```

Soft bits are now converted to hard bits.

```
In [10]:
```bits = (softbits > 0)*1

Now we undo GMSK precoding, since we want to export the OQPSK bits. This part is tricky and use the knowledge of the ASM to correct for phase ambiguities.

As a first step, since the GMSK bits are differential phase, a cumulative sum of them gives phase, which can be converted to the QPSK constellation.

```
In [11]:
```decbits = np.cumsum(np.int32(2*bits-1))%4
decbits[::2] = (decbits[::2] == 1)*1
decbits[1::2] = (decbits[1::2] == 0)*1

```
In [12]:
```decbits_inv = decbits.copy()
decbits_inv[::2] ^= 1

```
In [13]:
```asm = np.unpackbits(np.array([0x03,0x47,0x76,0xC7,0x27,0x28,0x95,0xB0], dtype='uint8'))
asm_straight_corr = scipy.signal.correlate(2*decbits-1, 2*asm.astype('float')-1)
asm_inv_corr = scipy.signal.correlate(2*decbits_inv-1, 2*asm.astype('float')-1)
plt.figure()
plt.plot(asm_straight_corr)
plt.figure()
plt.plot(asm_inv_corr);

```
```

```
In [14]:
```corr_positions = np.sort(np.concatenate((np.where(np.abs(asm_straight_corr) > 40)[0], np.where(np.abs(asm_inv_corr) > 40)[0])))
with open('/tmp/dslwp_bits', 'w') as file:
for st in corr_positions:
if np.abs(asm_straight_corr[st]) > 40:
b = np.sign(asm_straight_corr[st])*(2*decbits[st+1:st+1+3576]-1)
else:
b = np.sign(asm_inv_corr[st])*(2*decbits_inv[st+1:st+1+3576]-1)
b.astype('float32').tofile(file)

Total count of Turbo codewords found:

```
In [15]:
```corr_positions.size

```
Out[15]:
```

Distances between ASMs:

```
In [16]:
```np.diff(corr_positions)

```
Out[16]:
```

All seems OK

We turn to the method of correlating against the precoded ASM. To validate the method above.

We construct the CCSDS ASM used to mark the beginning of each Turbo coded word. The ASM is precoded as indicated by the CCSDS standard. See this post for more information.

```
In [17]:
```asm_diff = asm[:-1] ^ asm[1:]
asm_diff[1::2] ^= 1

```
In [18]:
```asm_corr = scipy.signal.correlate(2*bits-1, 2*asm_diff.astype('float')-1)
plt.plot(asm_corr);

```
```

Note that all the ASMs have correlated without any bit errors.

```
In [19]:
```asm_corr[asm_corr > 40]

```
Out[19]:
```

```
In [20]:
```asm_corr[asm_corr > 40].size

```
Out[20]:
```

We now look at the distances between the ASMs. DSLWP-B uses Turbo codewords of 3576 symbols. Since the ASM is 64 bits long and Turbo codewords are transmitted back-to-back, without any gap between them, we expect a distance of 3640 bits.

Note that before we have stated that the ASM is 63 bits long. This is because the GMSK precoder is differential, so the first bit of the ASM is not defined, as it depends on the preceding data. Thus, we only use the 63 bits of the ASM that are fixed when doing the correlation.

Below we show the distances between the ASMs.

```
In [21]:
```np.diff(np.where(asm_corr > 40)[0])

```
Out[21]:
```

We see that there are no problems with Turbo words being cut short.

Now we turn to the investigation of the frequency jump that happened about bit 20000.

First we extract the relevant piece of the signal.

```
In [22]:
```start,end = 18000, 30000
xx = np.fromfile('/home/daniel/Descargas/DSLWP-B_PI9CAM_2018-08-12T08_32_25_436.4MHz_40ksps_complex.raw', dtype='complex64')[start*160:end*160]
xx = xx * np.exp(-1j*2*np.pi*np.arange(xx.size)*f/fs).astype('complex64')
xx = scipy.signal.lfilter(h, 1, xx).astype('complex64')

Now we plot the phase of the signal. The change in slope shows the frequency jump.

```
In [23]:
```plt.plot(np.arange(start, end, 1/160), scipy.signal.detrend(np.unwrap(np.angle(xx))))
plt.title('Dentrended and unwrapped phase')
plt.ylabel('Phase (rad)')
plt.xlabel('Bit');

```
```