<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Detect subframe preamble by majority voting amongst the possible offsets.
Step1:  Most subframes do not have valid parity, as shown below. We use a weaker heuristic, where only parity of TLM and HOW words are required to be valid. With this criterion, all the subframes are valid, except the first few and last subframes.
Step2:  Analysis of the TOW in valid frames. It ranges between 298338 and 301728, corresponding to 10
Step3:  The subframe ID in the HOW word cycles as usual.
Step4:  Alert and anti-spoofing flags in the HOW word are not set.
Step5:  For subframes 1 (WN, health and clock), 2 and 3 (ephemeris), a filler of alternating 1's and 0's is transmitted in all the words after the HOW (including the parity bits). This makes parity invalid for these subframes.
Step6:  On the other hand, the parity for subframes 4 and 5 (almanacs) is correct.
Step7:  Data ID field for subframes 4 and 5 has the nominal value 01.
Step8:  The SVID in subframes 4 and 5 follows the nominal schedule, except that SVID 4 has been replaced with 0 to indicate dummy SV. This is normal, since PRN 4 is not currently assigned.
Step9:  For subframe 5, we omit the study of pages 1 through 24, which contain almanac data and we assume to be valid. We study page 25, which is marked by SVID 51 and contains SV health.
Step10:  The t_oa and WN_a for page 25 correspond to times near the beginning and end of GPS week 2059.
Step11:  SV health in subframe 5 page 25 indicates that all SV except SV 4 are healthy.
Step12:  The anti-spoofing and SV configurations flags in subframe 4 page 25 indicate that AS is on for all SVs and different signal capabilities for different SVs.
Step13:  The health flags in subframe 4 page 25 indicate that SV 25 to 32 are all healthy.
Step14:  Below we show t_oa for almanac entries in subframes 4 and 5.
<ASSISTANT_TASK:>
Python Code:
preamble = np.array([1,0,0,0,1,0,1,1], dtype = 'uint8')
preamble_detect = np.where(np.abs(np.correlate(2*bits.astype('int')-1, 2*preamble.astype('int')-1)) == 8)[0]
preamble_offset = np.argmax(np.histogram(preamble_detect % subframe_size, bins = np.arange(0,subframe_size))[0])
subframes = bits[preamble_offset:]
subframes = subframes[:subframes.size//subframe_size*subframe_size].reshape((-1,subframe_size))
words = subframes.reshape((-1,word_size))
# Last bits from previous word, used for parity calculations
words_last = np.roll(words[:,-1], 1)
words_prelast = np.roll(words[:,-2], 1)
# Correct data using last bit from previous word
words_data = words[:, :word_data_size] ^ words_last.reshape((-1,1))
subframes_data = words_data.reshape((-1,subframe_data_size))
# Parity checks for each of the bits (0 means valid)
parity0 = np.bitwise_xor.reduce(words_data[:, np.array([1,2,3,5,6,10,11,12,13,14,17,18,20,23])-1], axis = 1) ^ words_prelast ^ words[:,word_data_size]
parity1 = np.bitwise_xor.reduce(words_data[:, np.array([2,3,4,6,7,11,12,13,14,15,18,19,21,24])-1], axis = 1) ^ words_last ^ words[:,word_data_size+1]
parity2 = np.bitwise_xor.reduce(words_data[:, np.array([1,3,4,5,7,8,12,13,14,15,16,19,20,22])-1], axis = 1) ^ words_prelast ^ words[:,word_data_size+2]
parity3 = np.bitwise_xor.reduce(words_data[:, np.array([2,4,5,6,8,9,13,14,15,16,17,20,21,23])-1], axis = 1) ^ words_last ^ words[:,word_data_size+3]
parity4 = np.bitwise_xor.reduce(words_data[:, np.array([1,3,5,6,7,9,10,14,15,16,17,18,21,22,24])-1], axis = 1) ^ words_last ^ words[:,word_data_size+4]
parity5 = np.bitwise_xor.reduce(words_data[:, np.array([3,5,6,8,9,10,11,13,15,19,22,23,24])-1], axis = 1) ^ words_prelast ^ words[:,word_data_size+5]
# Parity check for word
parity = parity0 | parity1 | parity2 | parity3 | parity4 | parity5
# Parity check for subframe
parity_subframe = np.any(parity.reshape((-1,subframe_size//word_size)), axis = 1)
Explanation: Detect subframe preamble by majority voting amongst the possible offsets.
End of explanation
parity_subframe
correct_frames = (parity[::10] == 0) & (parity[1::10] == 0)
plt.plot(correct_frames)
Explanation: Most subframes do not have valid parity, as shown below. We use a weaker heuristic, where only parity of TLM and HOW words are required to be valid. With this criterion, all the subframes are valid, except the first few and last subframes.
End of explanation
tow = np.sum(words_data[1::10,:17].astype('int') * 2**np.arange(16,-1,-1), axis = 1) * 6
plt.plot(np.arange(tow.size)[correct_frames], tow[correct_frames])
Explanation: Analysis of the TOW in valid frames. It ranges between 298338 and 301728, corresponding to 10:52:18 and 11:48:48 on Wednesday (2019-06-26).
End of explanation
subframe_id = np.packbits(words_data[1::10,19:22], axis = 1).ravel() >> 5
subframe_id[correct_frames]
Explanation: The subframe ID in the HOW word cycles as usual.
End of explanation
np.any(words_data[1::10,17:19][correct_frames])
Explanation: Alert and anti-spoofing flags in the HOW word are not set.
End of explanation
filler_subframe = subframes[correct_frames & (subframe_id == 1), 60:][0,:]
filler_subframe
np.any(subframes[correct_frames & (subframe_id <= 3), 60:] ^ filler_subframe, axis = 1)
np.all(parity_subframe[correct_frames & (subframe_id <= 3)])
Explanation: For subframes 1 (WN, health and clock), 2 and 3 (ephemeris), a filler of alternating 1's and 0's is transmitted in all the words after the HOW (including the parity bits). This makes parity invalid for these subframes.
End of explanation
np.all(parity_subframe[correct_frames & (subframe_id >= 4)])
Explanation: On the other hand, the parity for subframes 4 and 5 (almanacs) is correct.
End of explanation
np.any(subframes_data[correct_frames & (subframe_id >= 4), 2*word_data_size:2*word_data_size+2] ^ np.array([0,1]))
Explanation: Data ID field for subframes 4 and 5 has the nominal value 01.
End of explanation
svid_subframe4 = np.packbits(subframes_data[correct_frames & (subframe_id == 4), 2*word_data_size+2:2*word_data_size+8], axis = 1).ravel() >> 2
svid_subframe5 = np.packbits(subframes_data[correct_frames & (subframe_id == 5), 2*word_data_size+2:2*word_data_size+8], axis = 1).ravel() >> 2
svid_subframe4
svid_subframe5
Explanation: The SVID in subframes 4 and 5 follows the nominal schedule, except that SVID 4 has been replaced with 0 to indicate dummy SV. This is normal, since PRN 4 is not currently assigned.
End of explanation
subframe5_page25 = subframes_data[correct_frames & (subframe_id == 5), :][svid_subframe5 == 51, :]
Explanation: For subframe 5, we omit the study of pages 1 through 24, which contain almanac data and we assume to be valid. We study page 25, which is marked by SVID 51 and contains SV health.
End of explanation
toa = np.packbits(subframe5_page25[:,2*word_data_size+8:2*word_data_size+16], axis = 1).ravel().astype('int') * 2**12
wna = np.packbits(subframe5_page25[:,2*word_data_size+16:2*word_data_size+24], axis = 1).ravel() + 2048
toa
wna
Explanation: The t_oa and WN_a for page 25 correspond to times near the beginning and end of GPS week 2059.
End of explanation
subframe5_page25[:,2*word_data_size+24:][:,:6*6*4:].reshape((-1,6*4,6))
Explanation: SV health in subframe 5 page 25 indicates that all SV except SV 4 are healthy.
End of explanation
subframe4_page25 = subframes_data[correct_frames & (subframe_id == 4), :][svid_subframe4 == 63, :]
anti_spoofing = subframe4_page25[:,2*word_data_size+8:][:,:32*4].reshape((-1,32,4))
anti_spoofing
Explanation: The anti-spoofing and SV configurations flags in subframe 4 page 25 indicate that AS is on for all SVs and different signal capabilities for different SVs.
End of explanation
health = subframe4_page25[:,2*word_data_size+8+32*4+2:][:,:6*8].reshape((-1,8,6))
health
Explanation: The health flags in subframe 4 page 25 indicate that SV 25 to 32 are all healthy.
End of explanation
np.packbits(subframes_data[correct_frames & (subframe_id == 5), :][svid_subframe5 <= 24, 3*word_data_size:3*word_data_size+8], axis = 1) * 2**12
np.packbits(subframes_data[correct_frames & (subframe_id == 4), :][svid_subframe4 <= 32, 3*word_data_size:3*word_data_size+8], axis = 1) * 2**12
Explanation: Below we show t_oa for almanac entries in subframes 4 and 5.
End of explanation
<END_TASK>