<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Recovering SSDV image 0xAB, transmitted on 2019-02-27 08
Step1:  We load a file containing the relevant GMSK transmission. The recording was done at the Dwingeloo radiotelescope and can be obtained here. Remember to edit the path below to point to the correct file.
Step2:  The 500bps GMSK signal is converted down to baseband and lowpass filtered to 1600Hz.
Step3:  Perform arctangent FSK demodulation.
Step4:  Correct for phase wrapping.
Step5:  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.
Step6:  Soft bits are now converted to hard bits.
Step7:  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.
Step8:  However we need to consider another branch in which I and Q are swapped, so the sign of one of them is inverted (consider swapping I and Q versus multiplying by 1j).
Step9:  We correlate both branches against the uncoded ASM. Note the correlation can happen on either branch and have either sign.
Step10:  We take note of all correlation positions, the branch where they happen and their sign to correct the phase and output packets that can be fed into the GNU Radio decoder. Note that we're always dealing with hard bits for simplicity.
Step11:  Now we plot the phase of the signal in the problematic zone. The change in slope shows the frequency jump.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
fs = 40000 # sample rate
Explanation: Recovering SSDV image 0xAB, transmitted on 2019-02-27 08:20
End of explanation
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)
Explanation: We load a file containing the relevant GMSK transmission. The recording was done at the Dwingeloo radiotelescope and can be obtained here. Remember to edit the path below to point to the correct file.
End of explanation
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')
Explanation: The 500bps GMSK signal is converted down to baseband and lowpass filtered to 1600Hz.
End of explanation
s = np.diff(np.angle(x).astype('float32'))
Explanation: Perform arctangent FSK demodulation.
End of explanation
s[s > np.pi] -= 2*np.pi
s[s < -np.pi] += 2*np.pi
Explanation: Correct for phase wrapping.
End of explanation
plt.plot(s)
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]);
Explanation: 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.
End of explanation
bits = (softbits > 0)*1
Explanation: Soft bits are now converted to hard bits.
End of explanation
decbits = np.cumsum(np.int32(2*bits-1))%4
decbits[::2] = (decbits[::2] == 1)*1
decbits[1::2] = (decbits[1::2] == 0)*1
Explanation: 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.
End of explanation
decbits_inv = decbits.copy()
decbits_inv[::2] ^= 1
Explanation: However we need to consider another branch in which I and Q are swapped, so the sign of one of them is inverted (consider swapping I and Q versus multiplying by 1j).
End of explanation
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);
Explanation: We correlate both branches against the uncoded ASM. Note the correlation can happen on either branch and have either sign.
End of explanation
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)
Explanation: We take note of all correlation positions, the branch where they happen and their sign to correct the phase and output packets that can be fed into the GNU Radio decoder. Note that we're always dealing with hard bits for simplicity.
End of explanation
phase = np.unwrap(np.angle(x[1000000:3000000]).astype('float32'))
plt.plot(scipy.signal.detrend(phase))
Explanation: Now we plot the phase of the signal in the problematic zone. The change in slope shows the frequency jump.
End of explanation
<END_TASK>