Encoding/Decoding CCSDS (255,223) messages with reedsolomon

[Jazzinghen]

I have been trying to use the library to encode/decode messages using the CCSDS standard for the last couple days and I cannot seem to be able to get results in line with a reference implementation. I am already handling conventional to double-base conversion, but it seems impossible to set the proper parameters for the encoder.

The parameters for the CCSDS encoder are the following:

• Symbol size: 8
• Message length: 255
• Nsym: 32
• Generator polynomial: 0x187 (x^8 + x^7 + x^2 + x + 1)
• First consecutive root: 112
• Primitive element: 11

I tried using the library with this configuration but it doesn't seem to generate the correct output:

encoder = crs.RSCodec(32,255,112,0x187,11)

Which has the following members:

encoder.gen -> 0x01,0x92,0x7a,0xf7,0x77,0x13,0x05,0xfe,0x90,0xec,0xdc,0xfd,0x2a,0xa4,0x24,0xd8,0x93,0xd8,0x24,0xa4,0x2a,0xfd,0xdc,0xec,0x90,0xfe,0x05,0x13,0x77,0xf7,0x7a,0x92,0x01
encoder.generator -> 11
encoder.prim -> 391

I have been also trying to debug the issue by checking all the steps, but due to my lack of knowledge on Reed-Solomon and the differences in implementation between the two libraries it is almost impossible for me.

I also tried to use detect_reedsolomon_parameters from pyFileFixty, and it cannot find the proper parameters:

print(detect_reedsolomon_parameters(test_data,test_message,gen_list=[2,3,5,11]))
Found closest set of parameters, with Hamming distance 30:
gen_nb=2 prim=285(0x11d) fcr=38
gen_nb=2 prim=285(0x11d) fcr=43
gen_nb=2 prim=285(0x11d) fcr=167
gen_nb=2 prim=299(0x12b) fcr=59
gen_nb=2 prim=299(0x12b) fcr=168
gen_nb=2 prim=299(0x12b) fcr=182
gen_nb=2 prim=333(0x14d) fcr=98
gen_nb=2 prim=333(0x14d) fcr=180
gen_nb=2 prim=333(0x14d) fcr=200
gen_nb=2 prim=351(0x15f) fcr=227
gen_nb=2 prim=357(0x165) fcr=84
gen_nb=2 prim=357(0x165) fcr=94
gen_nb=2 prim=357(0x165) fcr=213
gen_nb=2 prim=391(0x187) fcr=72
gen_nb=2 prim=391(0x187) fcr=208
gen_nb=2 prim=391(0x187) fcr=233
gen_nb=2 prim=397(0x18d) fcr=250
gen_nb=2 prim=425(0x1a9) fcr=4
gen_nb=2 prim=425(0x1a9) fcr=71
gen_nb=2 prim=451(0x1c3) fcr=3
gen_nb=2 prim=451(0x1c3) fcr=51
gen_nb=2 prim=451(0x1c3) fcr=113
gen_nb=2 prim=451(0x1c3) fcr=155
gen_nb=2 prim=463(0x1cf) fcr=91
gen_nb=2 prim=487(0x1e7) fcr=20
gen_nb=2 prim=501(0x1f5) fcr=9
gen_nb=3 prim=333(0x14d) fcr=98
gen_nb=3 prim=333(0x14d) fcr=180
gen_nb=3 prim=333(0x14d) fcr=200
gen_nb=3 prim=351(0x15f) fcr=227
gen_nb=3 prim=357(0x165) fcr=84
gen_nb=3 prim=357(0x165) fcr=94
gen_nb=3 prim=357(0x165) fcr=213
gen_nb=3 prim=397(0x18d) fcr=250
gen_nb=3 prim=419(0x1a3) fcr=60
gen_nb=3 prim=419(0x1a3) fcr=145
gen_nb=3 prim=419(0x1a3) fcr=230
gen_nb=3 prim=425(0x1a9) fcr=4
gen_nb=3 prim=425(0x1a9) fcr=71
gen_nb=3 prim=451(0x1c3) fcr=3
gen_nb=3 prim=451(0x1c3) fcr=51
gen_nb=3 prim=451(0x1c3) fcr=113
gen_nb=3 prim=451(0x1c3) fcr=155
gen_nb=3 prim=501(0x1f5) fcr=9
gen_nb=5 prim=333(0x14d) fcr=98
gen_nb=5 prim=333(0x14d) fcr=180
gen_nb=5 prim=333(0x14d) fcr=200
gen_nb=5 prim=351(0x15f) fcr=227
gen_nb=5 prim=357(0x165) fcr=84
gen_nb=5 prim=357(0x165) fcr=94
gen_nb=5 prim=357(0x165) fcr=213
gen_nb=5 prim=397(0x18d) fcr=250
gen_nb=5 prim=419(0x1a3) fcr=60
gen_nb=5 prim=419(0x1a3) fcr=145
gen_nb=5 prim=419(0x1a3) fcr=230
gen_nb=5 prim=425(0x1a9) fcr=4
gen_nb=5 prim=425(0x1a9) fcr=71
gen_nb=5 prim=451(0x1c3) fcr=3
gen_nb=5 prim=451(0x1c3) fcr=51
gen_nb=5 prim=451(0x1c3) fcr=113
gen_nb=5 prim=451(0x1c3) fcr=155
gen_nb=5 prim=501(0x1f5) fcr=9
gen_nb=11 prim=351(0x15f) fcr=227
gen_nb=11 prim=357(0x165) fcr=84
gen_nb=11 prim=357(0x165) fcr=94
gen_nb=11 prim=357(0x165) fcr=213
gen_nb=11 prim=391(0x187) fcr=72
gen_nb=11 prim=391(0x187) fcr=208
gen_nb=11 prim=391(0x187) fcr=233
gen_nb=11 prim=463(0x1cf) fcr=91
gen_nb=11 prim=477(0x1dd) fcr=39
gen_nb=11 prim=477(0x1dd) fcr=124
gen_nb=11 prim=477(0x1dd) fcr=209

As an example, this is the output of an application using libfec to compute the CCSDS parity:

Conventional Message:
	0x00,0xcc,0xac,0x60,0x79,0xb5,0xd5,0x19,0xf0,0x3c,0x5c,0x90,0x89,0x45,0x25,0xe9,
	0xfd,0x31,0x51,0x9d,0x84,0x48,0x28,0xe4,0x0d,0xc1,0xa1,0x6d,0x74,0xb8,0xd8,0x14,
	0x2e,0xe2,0x82,0x4e,0x57,0x9b,0xfb,0x37,0xde,0x12,0x72,0xbe,0xa7,0x6b,0x0b,0xc7,
	0xd3,0x1f,0x7f,0xb3,0xaa,0x66,0x06,0xca,0x23,0xef,0x8f,0x43,0x5a,0x96,0xf6,0x3a,
	0x42,0x8e,0xee,0x22,0x3b,0xf7,0x97,0x5b,0xb2,0x7e,0x1e,0xd2,0xcb,0x07,0x67,0xab,
	0xbf,0x73,0x13,0xdf,0xc6,0x0a,0x6a,0xa6,0x4f,0x83,0xe3,0x2f,0x36,0xfa,0x9a,0x56,
	0x6c,0xa0,0xc0,0x0c,0x15,0xd9,0xb9,0x75,0x9c,0x50,0x30,0xfc,0xe5,0x29,0x49,0x85,
	0x91,0x5d,0x3d,0xf1,0xe8,0x24,0x44,0x88,0x61,0xad,0xcd,0x01,0x18,0xd4,0xb4,0x78,
	0xc5,0x09,0x69,0xa5,0xbc,0x70,0x10,0xdc,0x35,0xf9,0x99,0x55,0x4c,0x80,0xe0,0x2c,
	0x38,0xf4,0x94,0x58,0x41,0x8d,0xed,0x21,0xc8,0x04,0x64,0xa8,0xb1,0x7d,0x1d,0xd1,
	0xeb,0x27,0x47,0x8b,0x92,0x5e,0x3e,0xf2,0x1b,0xd7,0xb7,0x7b,0x62,0xae,0xce,0x02,
	0x16,0xda,0xba,0x76,0x6f,0xa3,0xc3,0x0f,0xe6,0x2a,0x4a,0x86,0x9f,0x53,0x33,0xff,
	0x87,0x4b,0x2b,0xe7,0xfe,0x32,0x52,0x9e,0x77,0xbb,0xdb,0x17,0x0e,0xc2,0xa2,0x6e,
	0x7a,0xb6,0xd6,0x1a,0x03,0xcf,0xaf,0x63,0x8a,0x46,0x26,0xea,0xf3,0x3f,0x5f

Conventional Parity:
	0xab,0xc4,0x94,0x3f,0x0a,0xb4,0x52,0x78,0x37,0xc4,0xf9,0x69,0x6e,0x4f,0xa4,0x11,
	0xac,0x99,0xb6,0xe4,0xdd,0x40,0xfc,0x37,0x58,0x7a,0x8e,0x35,0xfb,0xa6,0x10,0x73

This is the output of reedsolo:

Conventional Message:
	0x00,0xcc,0xac,0x60,0x79,0xb5,0xd5,0x19,0xf0,0x3c,0x5c,0x90,0x89,0x45,0x25,0xe9,
	0xfd,0x31,0x51,0x9d,0x84,0x48,0x28,0xe4,0x0d,0xc1,0xa1,0x6d,0x74,0xb8,0xd8,0x14,
	0x2e,0xe2,0x82,0x4e,0x57,0x9b,0xfb,0x37,0xde,0x12,0x72,0xbe,0xa7,0x6b,0x0b,0xc7,
	0xd3,0x1f,0x7f,0xb3,0xaa,0x66,0x06,0xca,0x23,0xef,0x8f,0x43,0x5a,0x96,0xf6,0x3a,
	0x42,0x8e,0xee,0x22,0x3b,0xf7,0x97,0x5b,0xb2,0x7e,0x1e,0xd2,0xcb,0x07,0x67,0xab,
	0xbf,0x73,0x13,0xdf,0xc6,0x0a,0x6a,0xa6,0x4f,0x83,0xe3,0x2f,0x36,0xfa,0x9a,0x56,
	0x6c,0xa0,0xc0,0x0c,0x15,0xd9,0xb9,0x75,0x9c,0x50,0x30,0xfc,0xe5,0x29,0x49,0x85,
	0x91,0x5d,0x3d,0xf1,0xe8,0x24,0x44,0x88,0x61,0xad,0xcd,0x01,0x18,0xd4,0xb4,0x78,
	0xc5,0x09,0x69,0xa5,0xbc,0x70,0x10,0xdc,0x35,0xf9,0x99,0x55,0x4c,0x80,0xe0,0x2c,
	0x38,0xf4,0x94,0x58,0x41,0x8d,0xed,0x21,0xc8,0x04,0x64,0xa8,0xb1,0x7d,0x1d,0xd1,
	0xeb,0x27,0x47,0x8b,0x92,0x5e,0x3e,0xf2,0x1b,0xd7,0xb7,0x7b,0x62,0xae,0xce,0x02,
	0x16,0xda,0xba,0x76,0x6f,0xa3,0xc3,0x0f,0xe6,0x2a,0x4a,0x86,0x9f,0x53,0x33,0xff,
	0x87,0x4b,0x2b,0xe7,0xfe,0x32,0x52,0x9e,0x77,0xbb,0xdb,0x17,0x0e,0xc2,0xa2,0x6e,
	0x7a,0xb6,0xd6,0x1a,0x03,0xcf,0xaf,0x63,0x8a,0x46,0x26,0xea,0xf3,0x3f,0x5f

Conventional Parity:
	0x04,0x01,0x32,0x09,0xc2,0xe5,0x3e,0xd5,0xfc,0x4c,0xac,0x30,0xe7,0x8e,0x8c,0xe0
        0xef,0x9a,0xff,0x4a,0x3c,0xff,0xed,0x3a,0x58,0xbc,0xd4,0x63,0x4f,0xbb,0x6f,0xf9

Did I miss something? It is entirely possible given that I am not an expert in this field.

[lrq3000]

Thank you for your extensive detailing of this issue!

Could you also please provide the python code you wrote to run these tests? One possibility is that the encoding is mangling the input message somehow? Also are you using the low level functions? Because otherwise if you use the high level RSCodec object then it automatically pads input messages so this may be the culprit.

There can also be a difference in the implementation which is not obvious and may change results a lot.

[lrq3000]

With some code i will be able to reproduce the issue and test faster.

[Jazzinghen]

Roger. Just a second.
I wrote a lot of stuff in iPython, so I need to move it to a file.

[Jazzinghen]

Here you go, thanks for waiting! (For some I had to rename it to txt)

rs_test.txt

[Jazzinghen]

If you are interested this is the code I wrote to generate the data using libfec:

#include "fec.h"
#include <memory.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define MSG_SIZE 255
#define MAX_ERRORS 16
#define MSG_CODE MAX_ERRORS * 2
#define MSG_PAYLOAD MSG_SIZE - MSG_CODE

int main() {
  srand(time(NULL)); // use current time as seed for random generator

  unsigned char message[MSG_PAYLOAD] = {
      0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0A, 0x0B,
      0x0C, 0x0D, 0x0E, 0x0F, 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17,
      0x18, 0x19, 0x1A, 0x1B, 0x1C, 0x1D, 0x1E, 0x1F, 0x20, 0x21, 0x22, 0x23,
      0x24, 0x25, 0x26, 0x27, 0x28, 0x29, 0x2A, 0x2B, 0x2C, 0x2D, 0x2E, 0x2F,
      0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x3A, 0x3B,
      0x3C, 0x3D, 0x3E, 0x3F, 0x40, 0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47,
      0x48, 0x49, 0x4A, 0x4B, 0x4C, 0x4D, 0x4E, 0x4F, 0x50, 0x51, 0x52, 0x53,
      0x54, 0x55, 0x56, 0x57, 0x58, 0x59, 0x5A, 0x5B, 0x5C, 0x5D, 0x5E, 0x5F,
      0x60, 0x61, 0x62, 0x63, 0x64, 0x65, 0x66, 0x67, 0x68, 0x69, 0x6A, 0x6B,
      0x6C, 0x6D, 0x6E, 0x6F, 0x70, 0x71, 0x72, 0x73, 0x74, 0x75, 0x76, 0x77,
      0x78, 0x79, 0x7A, 0x7B, 0x7C, 0x7D, 0x7E, 0x7F, 0x80, 0x81, 0x82, 0x83,
      0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x8D, 0x8E, 0x8F,
      0x90, 0x91, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98, 0x99, 0x9A, 0x9B,
      0x9C, 0x9D, 0x9E, 0x9F, 0xA0, 0xA1, 0xA2, 0xA3, 0xA4, 0xA5, 0xA6, 0xA7,
      0xA8, 0xA9, 0xAA, 0xAB, 0xAC, 0xAD, 0xAE, 0xAF, 0xB0, 0xB1, 0xB2, 0xB3,
      0xB4, 0xB5, 0xB6, 0xB7, 0xB8, 0xB9, 0xBA, 0xBB, 0xBC, 0xBD, 0xBE, 0xBF,
      0xC0, 0xC1, 0xC2, 0xC3, 0xC4, 0xC5, 0xC6, 0xC7, 0xC8, 0xC9, 0xCA, 0xCB,
      0xCC, 0xCD, 0xCE, 0xCF, 0xD0, 0xD1, 0xD2, 0xD3, 0xD4, 0xD5, 0xD6, 0xD7,
      0xD8, 0xD9, 0xDA, 0xDB, 0xDC, 0xDD, 0xDE};
  /*
  for (size_t i = 0; i < MSG_PAYLOAD; ++i) {
    message[i] = rand() % 256;
  }
  */

  const int fcr = 128 - MAX_ERRORS;

  struct rs *rs = (struct rs *)init_rs_char(8, 0x187, fcr, 11, MSG_CODE, 0);

  printf("CCSDS(255,223) internals:\n");
  print_details_ccsds_239(rs);

  unsigned char parity[MSG_CODE];
  encode_rs_ccsds(message, parity, 0);

  printf("Message:");
  for (size_t i = 0; i < MSG_PAYLOAD; i++) {
    if ((i % 16) == 0)
      printf("\n\t");
    printf("0x%02x,", message[i]);
  };

  printf("\n\n");

  printf("Parity");
  for (size_t i = 0; i < MSG_CODE; i++) {
    if ((i % 16) == 0)
      printf("\n\t");
    printf("0x%02x,", parity[i]);
  }

  printf("\n");
  return 0;
};

[Jazzinghen]

I read a bit more around (especially the wiki books entry about ReedSolomon) and tested out a bit more.
Turns out that even if I use another library, the result is the same:

import galois

rs_field = galois.GF(2**8, irreducible_poly='x^8 + x^7 + x^2 + x + 1', primitive_element=_GEN_ALPHA)
reed_solomon = galois.ReedSolomon(_CHUNK_SIZE, _CHUNK_PAYLOAD, c=_FCR, field=rs_field)

galois_out = reed_solomon.encode(msg_data)

So I went back to the definition book to read more, given I had a bit more knowledge. This is the part that threw me off:

The code generator polynomial shall be:

$$g(x) = \prod_{j=128–E}^{127+E} (x – \alpha^{11j}) = \sum_{i=0}^{2E} G_ix^i$$

over $GF(2^8)$, where $F(\alpha) = 0$

The fact that 11 is used to multiply alpha, instead of being alpha, alarmed so I went back to check and, while I thought that there were some optimization and didn't put much thought about it, the fact that libfec is using an exp lookup table of the powers of 2 reduced using the prime polynomial made me try a thing.

I made libfec print the generator polynomial and fed it to the manual functions in reedsolo and the output at that point is the same:

(log, exp, field_char) = crs.init_tables(0x187)
gen = bytearray([0x01,0x5b,0x7f,0x56,0x10,0x1e,0x0d,0xeb,0x61,0xa5,0x08,0x2a,0x36,0x56,0xab,0x20,0x71,0x20,0xab,0x56,0x36,0x2a,0x08,0xa5,0x61,0xeb,0x0d,0x1e,0x10,0x56,0x7f,0x5b,0x01])
conv_message = crs.rs_encode_msg(conv_data,32,fcr=112,gen=gen)

This allows me to move forward, but if I can use this without using the "low-level" functionalities it would be great!

Let me know if I can give you more data to work on.

[Jazzinghen]

Ouch.
Sadly even in decoding the change in definition makes it impossible to use. I needed to patch rs_calc_syndromes like this to make it work (i.e. add a *11 to the exponential power):

def rs_calc_syndromes(msg, nsym, fcr=0, generator=2):
    return [0] + [crs.gf_poly_eval(msg, crs.gf_pow(generator, (i+fcr)*11)) for i in range(nsym)]

This makes the function return all zeroes in case of a correct message.

[Jazzinghen]

Oh, wait a second, I did it. I should have used $\alpha^11$ as generator, not 11. This means generating the tables for the specified prime (i.e. 0x187), then get exp_table[11].
Man, I feel pretty dumb.

[lrq3000]

Dear @Jazzinghen , thank you so much for investigating this issue and providing so much details, this certainly helps me a LOT! :D

I think I can implement a fix and will do so shortly, but bear in mind this is a workaround for a specific implementation that doesn't follow the original model. Indeed, there is no reason there is such an exponent here, it's not a parameter described anywhere (I mean it is in the book you cited - great finding! But it's not a common parameter and it's redundant with the j exponent when generating - it seems to me it's kinda like a salt or a shift to make their implementation unique).

[Jazzinghen]

No problem, I am sorry I flooded you with messages, I was very excited to have understood what was wrong. It surely is a very specific detail here (albeit for a "famous" implementation), so I am not sure if it calls for an API change or just a mention in the readme :D
Thank you for keeping up with me on this