Encoding/Decoding CCSDS (255,223) messages with reedsolomon about reedsolomon HOT 10 OPEN

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).

from reedsolomon.

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.

from reedsolomon.

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

from reedsolomon.

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

from reedsolomon.

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

rs_test.txt

from reedsolomon.

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;
};

from reedsolomon.

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.

from reedsolomon.

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.

from reedsolomon.

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.

from reedsolomon.

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

from reedsolomon.