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

In [2]:
frame_size = 223*5
frames = np.fromfile('bepi_colombo_frames.u8', dtype = 'uint8')
frames = frames[:frames.size//frame_size*frame_size].reshape((-1, frame_size))
frames.shape[0]


Out[2]:
1808

CCSDS TM frame CRC is CRC16_CCITT_FALSE from this calculator.


In [3]:
crc_table = [0x0000, 0x1021, 0x2042, 0x3063, 0x4084, 0x50a5, 0x60c6, 0x70e7,
    0x8108, 0x9129, 0xa14a, 0xb16b, 0xc18c, 0xd1ad, 0xe1ce, 0xf1ef,
    0x1231, 0x0210, 0x3273, 0x2252, 0x52b5, 0x4294, 0x72f7, 0x62d6,
    0x9339, 0x8318, 0xb37b, 0xa35a, 0xd3bd, 0xc39c, 0xf3ff, 0xe3de,
    0x2462, 0x3443, 0x0420, 0x1401, 0x64e6, 0x74c7, 0x44a4, 0x5485,
    0xa56a, 0xb54b, 0x8528, 0x9509, 0xe5ee, 0xf5cf, 0xc5ac, 0xd58d,
    0x3653, 0x2672, 0x1611, 0x0630, 0x76d7, 0x66f6, 0x5695, 0x46b4,
    0xb75b, 0xa77a, 0x9719, 0x8738, 0xf7df, 0xe7fe, 0xd79d, 0xc7bc,
    0x48c4, 0x58e5, 0x6886, 0x78a7, 0x0840, 0x1861, 0x2802, 0x3823,
    0xc9cc, 0xd9ed, 0xe98e, 0xf9af, 0x8948, 0x9969, 0xa90a, 0xb92b,
    0x5af5, 0x4ad4, 0x7ab7, 0x6a96, 0x1a71, 0x0a50, 0x3a33, 0x2a12,
    0xdbfd, 0xcbdc, 0xfbbf, 0xeb9e, 0x9b79, 0x8b58, 0xbb3b, 0xab1a,
    0x6ca6, 0x7c87, 0x4ce4, 0x5cc5, 0x2c22, 0x3c03, 0x0c60, 0x1c41,
    0xedae, 0xfd8f, 0xcdec, 0xddcd, 0xad2a, 0xbd0b, 0x8d68, 0x9d49,
    0x7e97, 0x6eb6, 0x5ed5, 0x4ef4, 0x3e13, 0x2e32, 0x1e51, 0x0e70,
    0xff9f, 0xefbe, 0xdfdd, 0xcffc, 0xbf1b, 0xaf3a, 0x9f59, 0x8f78,
    0x9188, 0x81a9, 0xb1ca, 0xa1eb, 0xd10c, 0xc12d, 0xf14e, 0xe16f,
    0x1080, 0x00a1, 0x30c2, 0x20e3, 0x5004, 0x4025, 0x7046, 0x6067,
    0x83b9, 0x9398, 0xa3fb, 0xb3da, 0xc33d, 0xd31c, 0xe37f, 0xf35e,
    0x02b1, 0x1290, 0x22f3, 0x32d2, 0x4235, 0x5214, 0x6277, 0x7256,
    0xb5ea, 0xa5cb, 0x95a8, 0x8589, 0xf56e, 0xe54f, 0xd52c, 0xc50d,
    0x34e2, 0x24c3, 0x14a0, 0x0481, 0x7466, 0x6447, 0x5424, 0x4405,
    0xa7db, 0xb7fa, 0x8799, 0x97b8, 0xe75f, 0xf77e, 0xc71d, 0xd73c,
    0x26d3, 0x36f2, 0x0691, 0x16b0, 0x6657, 0x7676, 0x4615, 0x5634,
    0xd94c, 0xc96d, 0xf90e, 0xe92f, 0x99c8, 0x89e9, 0xb98a, 0xa9ab,
    0x5844, 0x4865, 0x7806, 0x6827, 0x18c0, 0x08e1, 0x3882, 0x28a3,
    0xcb7d, 0xdb5c, 0xeb3f, 0xfb1e, 0x8bf9, 0x9bd8, 0xabbb, 0xbb9a,
    0x4a75, 0x5a54, 0x6a37, 0x7a16, 0x0af1, 0x1ad0, 0x2ab3, 0x3a92,
    0xfd2e, 0xed0f, 0xdd6c, 0xcd4d, 0xbdaa, 0xad8b, 0x9de8, 0x8dc9,
    0x7c26, 0x6c07, 0x5c64, 0x4c45, 0x3ca2, 0x2c83, 0x1ce0, 0x0cc1,
    0xef1f, 0xff3e, 0xcf5d, 0xdf7c, 0xaf9b, 0xbfba, 0x8fd9, 0x9ff8,
    0x6e17, 0x7e36, 0x4e55, 0x5e74, 0x2e93, 0x3eb2, 0x0ed1, 0x1ef0]

def crc16_ccitt_false(data):
    crc = 0xffff
    for d in data:
        tbl_idx = ((crc >> 8) ^ d) & 0xff
        crc = (crc_table[tbl_idx] ^ (crc << 8)) & 0xffff
    return crc & 0xffff

In [4]:
crc_ok = np.array([crc16_ccitt_false(f) for f in frames]) == 0

In [5]:
np.average(crc_ok)


Out[5]:
0.9994469026548672

In [6]:
np.sum(~crc_ok)


Out[6]:
1

In [7]:
plt.figure(figsize = (15,15), facecolor = 'w')
plt.imshow(frames[crc_ok,:], interpolation = 'nearest', aspect = 0.4)


Out[7]:
<matplotlib.image.AxesImage at 0x7f3ed4c7d6a0>

In [8]:
plt.figure(figsize = (30,15), facecolor = 'w')
plt.imshow(frames[crc_ok,:16], interpolation = 'nearest', aspect = 0.1)


Out[8]:
<matplotlib.image.AxesImage at 0x7f3ed4dc68d0>

In [9]:
plt.figure(figsize = (30,15), facecolor = 'w')
plt.imshow(frames[crc_ok,-16:], interpolation = 'nearest', aspect = 0.1)


Out[9]:
<matplotlib.image.AxesImage at 0x7f3ed4da2550>

In [10]:
TMPrimaryHeader = BitStruct('transfer_frame_version_number' / BitsInteger(2),
                            'spacecraft_id' / BitsInteger(10),
                            'virtual_channel_id' / BitsInteger(3),
                            'ocf_flag' / Flag,
                            'master_channel_frame_count' / BitsInteger(8),
                            'virtual_channel_frame_count' / BitsInteger(8),
                            'transfer_frame_secondary_header_flag' / Flag,
                            'synch_flag' / Flag,
                            'packet_order' / Flag,
                            'segment_length_id' / BitsInteger(2),
                            'first_header_pointer' / BitsInteger(11))

TMSecondaryHeaderID = BitStruct('transfer_frame_secondary_version_number' / BitsInteger(2),
                                'transfer_frame_secondary_header_length' / BitsInteger(6))

TMSecondaryHeader = Struct('id' / TMSecondaryHeaderID,
                            'data' / Bytes(this.id.transfer_frame_secondary_header_length + 1))

primary_headers = [TMPrimaryHeader.parse(bytes(f)) for f in frames[crc_ok]]
mcfc = np.array([p.master_channel_frame_count for p in primary_headers])
vcid = np.array([p.virtual_channel_id for p in primary_headers])
vcfc = np.array([p.virtual_channel_frame_count for p in primary_headers])

In [11]:
plt.figure(figsize = (8,4), facecolor = 'w')
plt.plot(mcfc)
plt.title('Master channel frame count')
plt.ylabel('Value')
plt.xlabel('Frame');



In [12]:
plt.figure(figsize = (8,4), facecolor = 'w')
plt.plot(np.diff(mcfc.astype('uint8')))
plt.title('Master channel frame count derivative')
plt.ylabel('Value')
plt.xlabel('Frame');



In [13]:
np.sum(np.diff(mcfc.astype('uint8')) - 1)


Out[13]:
8

In [14]:
vcids = np.unique(vcid)
for v in vcids:
    print(f'Virtual channel ID {v}: {np.sum(vcid == v)} frames')


Virtual channel ID 0: 38 frames
Virtual channel ID 1: 99 frames
Virtual channel ID 2: 39 frames
Virtual channel ID 7: 1631 frames

In [15]:
for v in vcids:
    plt.figure(figsize = (8,4), facecolor = 'w')
    style = '.-' if np.sum(vcid == v) < 100 else '-'
    plt.plot(vcfc[vcid == v], style)
    plt.title(f'Virtual channel {v} frame count')
    plt.ylabel('Value')
    plt.xlabel('Frame')
    plt.figure(figsize = (8,4), facecolor = 'w')
    plt.plot(np.diff(vcfc[vcid == v].astype('uint8')), style)
    plt.title(f'Virtual channel {v} frame count derivative')
    plt.ylabel('Value')
    plt.xlabel('Frame')



In [16]:
for v in vcids:
    plt.figure(figsize = (15,15), facecolor = 'w')
    total = np.sum(vcid == v)
    plt.imshow(frames[crc_ok][vcid == v], interpolation = 'nearest', aspect = 500/total)
    plt.title(f'Virtual channel {v}')


Virtual channel 7: idle data only

The first header pointer in virtual channel 7 frames is always 2046, which means idle data only.


In [17]:
np.unique([p.first_header_pointer for p in primary_headers if p.virtual_channel_id == 7])


Out[17]:
array([2046])

This is a typical virtual channel 7 frame header:


In [18]:
[p for p in primary_headers if p.virtual_channel_id == 7][0]


Out[18]:
Container(transfer_frame_version_number=0, spacecraft_id=816, virtual_channel_id=7, ocf_flag=True, master_channel_frame_count=233, virtual_channel_frame_count=105, transfer_frame_secondary_header_flag=True, synch_flag=False, packet_order=False, segment_length_id=3, first_header_pointer=2046)

We now show the nature of the payload in these idle frames.

First we extract payload bits in each frame in virtual channel 7, leaving gaps (as zeros) for missing frames.


In [19]:
seq_unwrap = np.int32(np.round(np.unwrap(vcfc[vcid == 7]/256*2*np.pi)/(2*np.pi)*256))
vc7_bits = np.zeros((seq_unwrap[-1]-seq_unwrap[0]+1, (frame_size - 16)*8), dtype = 'float32')
for j, frame in enumerate(frames[crc_ok][vcid == v]):
    vc7_bits[seq_unwrap[j]-seq_unwrap[0], :] = 2*np.unpackbits(frame[10:-6]).astype('float32')-1

Circular correlation analysis with FFT.


In [20]:
pn_len = 511
vc7_pn = vc7_bits.ravel()
vc7_pn = vc7_pn[:vc7_pn.size // pn_len * pn_len].reshape((-1,pn_len))
f = np.fft.fft(vc7_pn)
vc7_pn_corr = np.abs(np.fft.ifft(f * f[0,:].conjugate()))
vc7_pn_shift = np.argmax(vc7_pn_corr, axis = 1)
vc7_pn_corr = np.max(vc7_pn_corr, axis = 1)

In [21]:
plt.figure(figsize = (10,6), facecolor = 'w')
plt.plot(vc7_pn_shift)
plt.xlabel('Block')
plt.ylabel('Circular shift(bits)', color = 'C0')
ax2 = plt.gca().twinx()
ax2.plot(vc7_pn_corr, '.', color = 'C1')
plt.title(f'{pn_len} bit blocks circular correlation')
ax2.set_ylabel('Correlation magnitude (bits)', color = 'C1');


A more in detail look at the jumps in the correlation shift, cleaning the problems due to missing frames:


In [22]:
plt.plot(vc7_pn_shift[vc7_pn_corr>200])


Out[22]:
[<matplotlib.lines.Line2D at 0x7f3ed261ab00>]

In [23]:
d_shift = np.diff(vc7_pn_shift[vc7_pn_corr>200])
d_shift[d_shift != 0]


Out[23]:
array([ 308, -203,  308, -203, -203,  308])

In [24]:
308 - 511


Out[24]:
-203

In [25]:
idx_jumps = np.where(np.diff(vc7_pn_shift) % 511 == 308)
np.diff(idx_jumps)


Out[25]:
array([[4405, 4404, 4405, 4404, 4405]])

This motivates (and actually proves) the idea that the PN sequence is reset every 2250752 bits.


In [26]:
long_period = 2250752

We now build the same PN sequence, including the resets.

The PN sequence is generated by an LFSR with the polynomial x^9 + x^5 + 1.


In [27]:
lfsr = np.ones(9, dtype = 'uint8')
lfsr_out = np.empty(511, dtype = 'uint8')
for j in range(lfsr_out.size):
    out = lfsr[8] ^ lfsr[4]
    lfsr_out[j] = lfsr[8]
    lfsr = np.roll(lfsr, 1)
    lfsr[0] = out

This shows that the LFSR output (with an appropriate circular shift) coincides with the PN sequence.


In [28]:
lfsr_out_f = np.fft.fft(2*lfsr_out.astype('float32')-1)
lfsr_corr = np.fft.ifft(lfsr_out_f * f[0].conjugate()).real
np.argmax(lfsr_corr), np.max(lfsr_corr)


Out[28]:
(294, 510.9999999999999)

We create our long LFSR output by repeating and truncating to 2250752 bits.


In [29]:
lfsr_out_long = np.tile(2*lfsr_out.astype('int')-1, int(np.ceil(long_period/511)))[:long_period]

We just need to repeat lfsr_out_long and trim it appropriately so that it lines up with the bits in virtual channel 7.

We know that we need to take trim away 294 bits due to the circular correlation offset, but we also need to find how many integer multiples of 511 bits we need to trim in addition. We derive this from the location where the first PN sequence reset happens in virtual channel 7.


In [30]:
int(np.ceil(long_period/511)) - idx_jumps[0][0] - 2


Out[30]:
1806

So we trim 1806 full 511 bit periods.


In [31]:
lfsr_out_full = np.tile(lfsr_out_long, int(np.ceil(vc7_bits.size/lfsr_out_long.size)))[1806*511 + 294:]

Our LFSR resetting sequence only has -1's and 1's. The vc7_bits sequence has -1's and 1's, and 0's where we lost a frame. To check if they coincide (except for lost frames), we see that the product of both is never -1.


In [32]:
np.unique(lfsr_out_full)


Out[32]:
array([-1,  1])

In [33]:
np.unique(vc7_bits.ravel())


Out[33]:
array([-1.,  0.,  1.], dtype=float32)

In [34]:
np.any(vc7_bits.ravel() * lfsr_out_full[:vc7_bits.size] == -1)


Out[34]:
False

This proves the structure of the idle data in virtual channel 7: A 511 bit PN sequence generated by the polynomial x^9 + x^5 + 1 with initial state 0x1ff and resetting every 2250752 bits.

It turns out that 2250752 bits is exactly 256 frames. Let us see if the sequence resets exactly when the virtual channel frame counter wraps to zero. It suffices to check if the frames with virtual channel frame counter equal to zero have a payload which starts by the PN sequence at shift zero.


In [35]:
wrap_frames = frames[crc_ok][(vcid == v) & (vcfc == 0)]
wrap_frames_start = np.unpackbits(wrap_frames[:,10:], axis = 1)[:,:511]
np.any(lfsr_out ^ wrap_frames_start)


Out[35]:
False

Indeed the sequence resets anytime that the virtual channel frame counter is zero.