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

```
In [27]:
```fs = 40000 # sample rate

```
In [28]:
```skip = 3200 * fs * 8 # first transmission seems to have a missing samples from the recorder
skip = 5550 * fs * 8
length = 140 * fs
with open('/home/daniel/Descargas/DSLWP-B_PI9CAM_2019-02-27T07_26_01_436.4MHz_40ksps_complex.raw') as f:
f.seek(skip)
x = np.fromfile(f, dtype='complex64', count = length)

The 500bps GMSK signal is converted down to baseband and lowpass filtered to 1600Hz.

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

```
```

Perform arctangent FSK demodulation.

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

Correct for phase wrapping.

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

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 is not necessary either. The separation between the bits 1 and 0 is good enough for demodulation without bit errors.

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

```
In [8]:
```plt.plot(s)

```
Out[8]:
```

```
In [9]:
```phase = 40
softbits = s[np.int32(np.arange(phase,s.size, 80))]
softbits = softbits - 5e-8 * np.arange(softbits.size) - 1e-3 # correction for frequency offset and drift
plt.plot(softbits,'.')
plt.axhline(y = 0, color='green')
plt.ylim([-0.05,0.05]);

```
```

Soft bits are now converted to hard bits.

```
In [11]:
```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 [12]:
```decbits = np.cumsum(np.int32(2*bits-1))%4
decbits[::2] = (decbits[::2] == 1)*1
decbits[1::2] = (decbits[1::2] == 0)*1

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

```
In [14]:
```asm = np.unpackbits(np.array([0x03,0x47,0x76,0xC7,0x27,0x28,0x95,0xB0, 0xFC, 0xB8, 0x89, 0x38, 0xD8, 0xD7, 0x6A, 0x4F], 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 [26]:
```corr_positions = np.sort(np.concatenate((np.where(np.abs(asm_straight_corr) > 100)[0], np.where(np.abs(asm_inv_corr) > 100)[0])))
with open('/tmp/dslwp_bits', 'w') as file:
for st in corr_positions:
if np.abs(asm_straight_corr[st]) > 100:
b = np.sign(asm_straight_corr[st])*(2*decbits[st+1:st+1+7152]-1)
else:
b = np.sign(asm_inv_corr[st])*(2*decbits_inv[st+1:st+1+7152]-1)
b.astype('float32').tofile(file)

```
In [24]:
```phase = np.unwrap(np.angle(x[1000000:3000000]).astype('float32'))
plt.plot(scipy.signal.detrend(phase))

```
Out[24]:
```