Files
SDRPlusPlus/falcon9_decoder/src/libcorrect/reed-solomon/decode.c
2021-03-20 21:53:44 +01:00

509 lines
24 KiB
C

#include "correct/reed-solomon/encode.h"
// calculate all syndromes of the received polynomial at the roots of the generator
// because we're evaluating at the roots of the generator, and because the transmitted
// polynomial was made to be a product of the generator, we know that the transmitted
// polynomial is 0 at these roots
// any nonzero syndromes we find here are the values of the error polynomial evaluated
// at these roots, so these values give us a window into the error polynomial. if
// these syndromes are all zero, then we can conclude the error polynomial is also
// zero. if they're nonzero, then we know our message received an error in transit.
// returns true if syndromes are all zero
static bool reed_solomon_find_syndromes(field_t field, polynomial_t msgpoly, field_logarithm_t **generator_root_exp,
field_element_t *syndromes, size_t min_distance) {
bool all_zero = true;
memset(syndromes, 0, min_distance * sizeof(field_element_t));
for (unsigned int i = 0; i < min_distance; i++) {
// profiling reveals that this function takes about 50% of the cpu time of
// decoding. so, in order to speed it up a little, we precompute and save
// the successive powers of the roots of the generator, which are
// located in generator_root_exp
field_element_t eval = polynomial_eval_lut(field, msgpoly, generator_root_exp[i]);
if (eval) {
all_zero = false;
}
syndromes[i] = eval;
}
return all_zero;
}
// Berlekamp-Massey algorithm to find LFSR that describes syndromes
// returns number of errors and writes the error locator polynomial to rs->error_locator
static unsigned int reed_solomon_find_error_locator(correct_reed_solomon *rs, size_t num_erasures) {
unsigned int numerrors = 0;
memset(rs->error_locator.coeff, 0, (rs->min_distance + 1) * sizeof(field_element_t));
// initialize to f(x) = 1
rs->error_locator.coeff[0] = 1;
rs->error_locator.order = 0;
memcpy(rs->last_error_locator.coeff, rs->error_locator.coeff, (rs->min_distance + 1) * sizeof(field_element_t));
rs->last_error_locator.order = rs->error_locator.order;
field_element_t discrepancy;
field_element_t last_discrepancy = 1;
unsigned int delay_length = 1;
for (unsigned int i = rs->error_locator.order; i < rs->min_distance - num_erasures; i++) {
discrepancy = rs->syndromes[i];
for (unsigned int j = 1; j <= numerrors; j++) {
discrepancy = field_add(rs->field, discrepancy,
field_mul(rs->field, rs->error_locator.coeff[j], rs->syndromes[i - j]));
}
if (!discrepancy) {
// our existing LFSR describes the new syndrome as well
// leave it as-is but update the number of delay elements
// so that if a discrepancy occurs later we can eliminate it
delay_length++;
continue;
}
if (2 * numerrors <= i) {
// there's a discrepancy, but we still have room for more taps
// lengthen LFSR by one tap and set weight to eliminate discrepancy
// shift the last locator by the delay length, multiply by discrepancy,
// and divide by the last discrepancy
// we move down because we're shifting up, and this prevents overwriting
for (int j = rs->last_error_locator.order; j >= 0; j--) {
// the bounds here will be ok since we have a headroom of numerrors
rs->last_error_locator.coeff[j + delay_length] = field_div(
rs->field, field_mul(rs->field, rs->last_error_locator.coeff[j], discrepancy), last_discrepancy);
}
for (int j = delay_length - 1; j >= 0; j--) {
rs->last_error_locator.coeff[j] = 0;
}
// locator = locator - last_locator
// we will also update last_locator to be locator before this loop takes place
field_element_t temp;
for (int j = 0; j <= (rs->last_error_locator.order + delay_length); j++) {
temp = rs->error_locator.coeff[j];
rs->error_locator.coeff[j] =
field_add(rs->field, rs->error_locator.coeff[j], rs->last_error_locator.coeff[j]);
rs->last_error_locator.coeff[j] = temp;
}
unsigned int temp_order = rs->error_locator.order;
rs->error_locator.order = rs->last_error_locator.order + delay_length;
rs->last_error_locator.order = temp_order;
// now last_locator is locator before we started,
// and locator is (locator - (discrepancy/last_discrepancy) * x^(delay_length) * last_locator)
numerrors = i + 1 - numerrors;
last_discrepancy = discrepancy;
delay_length = 1;
continue;
}
// no more taps
// unlike the previous case, we are preserving last locator,
// but we'll update locator as before
// we're basically flattening the two loops from the previous case because
// we no longer need to update last_locator
for (int j = rs->last_error_locator.order; j >= 0; j--) {
rs->error_locator.coeff[j + delay_length] =
field_add(rs->field, rs->error_locator.coeff[j + delay_length],
field_div(rs->field, field_mul(rs->field, rs->last_error_locator.coeff[j], discrepancy),
last_discrepancy));
}
rs->error_locator.order = (rs->last_error_locator.order + delay_length > rs->error_locator.order)
? rs->last_error_locator.order + delay_length
: rs->error_locator.order;
delay_length++;
}
return rs->error_locator.order;
}
// find the roots of the error locator polynomial
// Chien search
bool reed_solomon_factorize_error_locator(field_t field, unsigned int num_skip, polynomial_t locator_log, field_element_t *roots,
field_logarithm_t **element_exp) {
// normally it'd be tricky to find all the roots
// but, the finite field is awfully finite...
// just brute force search across every field element
unsigned int root = num_skip;
memset(roots + num_skip, 0, (locator_log.order) * sizeof(field_element_t));
for (field_operation_t i = 0; i < 256; i++) {
// we make two optimizations here to help this search go faster
// a) we have precomputed the first successive powers of every single element
// in the field. we need at most n powers, where n is the largest possible
// degree of the error locator
// b) we have precomputed the error locator polynomial in log form, which
// helps reduce some lookups that would be done here
if (!polynomial_eval_log_lut(field, locator_log, element_exp[i])) {
roots[root] = (field_element_t)i;
root++;
}
}
// this is where we find out if we are have too many errors to recover from
// berlekamp-massey may have built an error locator that has 0 discrepancy
// on the syndromes but doesn't have enough roots
return root == locator_log.order + num_skip;
}
// use error locator and syndromes to find the error evaluator polynomial
void reed_solomon_find_error_evaluator(field_t field, polynomial_t locator, polynomial_t syndromes,
polynomial_t error_evaluator) {
// the error evaluator, omega(x), is S(x)*Lamba(x) mod x^(2t)
// where S(x) is a polynomial constructed from the syndromes
// S(1) + S(2)*x + ... + S(2t)*x(2t - 1)
// and Lambda(x) is the error locator
// the modulo is implicit here -- we have limited the max length of error_evaluator,
// which polynomial_mul will interpret to mean that it should not compute
// powers larger than that, which is the same as performing mod x^(2t)
polynomial_mul(field, locator, syndromes, error_evaluator);
}
// use error locator, error roots and syndromes to find the error values
// that is, the elements in the finite field which can be added to the received
// polynomial at the locations of the error roots in order to produce the
// transmitted polynomial
// forney algorithm
void reed_solomon_find_error_values(correct_reed_solomon *rs) {
// error value e(j) = -(X(j)^(1-c) * omega(X(j)^-1))/(lambda'(X(j)^-1))
// where X(j)^-1 is a root of the error locator, omega(X) is the error evaluator,
// lambda'(X) is the first formal derivative of the error locator,
// and c is the first consecutive root of the generator used in encoding
// first find omega(X), the error evaluator
// we generate S(x), the polynomial constructed from the roots of the syndromes
// this is *not* the polynomial constructed by expanding the products of roots
// S(x) = S(1) + S(2)*x + ... + S(2t)*x(2t - 1)
polynomial_t syndrome_poly;
syndrome_poly.order = rs->min_distance - 1;
syndrome_poly.coeff = rs->syndromes;
memset(rs->error_evaluator.coeff, 0, (rs->error_evaluator.order + 1) * sizeof(field_element_t));
reed_solomon_find_error_evaluator(rs->field, rs->error_locator, syndrome_poly, rs->error_evaluator);
// now find lambda'(X)
rs->error_locator_derivative.order = rs->error_locator.order - 1;
polynomial_formal_derivative(rs->field, rs->error_locator, rs->error_locator_derivative);
// calculate each e(j)
for (unsigned int i = 0; i < rs->error_locator.order; i++) {
if (rs->error_roots[i] == 0) {
continue;
}
rs->error_vals[i] = field_mul(
rs->field, field_pow(rs->field, rs->error_roots[i], rs->first_consecutive_root - 1),
field_div(
rs->field, polynomial_eval_lut(rs->field, rs->error_evaluator, rs->element_exp[rs->error_roots[i]]),
polynomial_eval_lut(rs->field, rs->error_locator_derivative, rs->element_exp[rs->error_roots[i]])));
}
}
void reed_solomon_find_error_locations(field_t field, field_logarithm_t generator_root_gap,
field_element_t *error_roots, field_logarithm_t *error_locations,
unsigned int num_errors, unsigned int num_skip) {
for (unsigned int i = 0; i < num_errors; i++) {
// the error roots are the reciprocals of the error locations, so div 1 by them
// we do mod 255 here because the log table aliases at index 1
// the log of 1 is both 0 and 255 (alpha^255 = alpha^0 = 1)
// for most uses it makes sense to have log(1) = 255, but in this case
// we're interested in a byte index, and the 255th index is not even valid
// just wrap it back to 0
if (error_roots[i] == 0) {
continue;
}
field_operation_t loc = field_div(field, 1, error_roots[i]);
for (field_operation_t j = 0; j < 256; j++) {
if (field_pow(field, j, generator_root_gap) == loc) {
error_locations[i] = field.log[j];
break;
}
}
}
}
// erasure method -- take given locations and convert to roots
// this is the inverse of reed_solomon_find_error_locations
static void reed_solomon_find_error_roots_from_locations(field_t field, field_logarithm_t generator_root_gap,
const field_logarithm_t *error_locations,
field_element_t *error_roots, unsigned int num_errors) {
for (unsigned int i = 0; i < num_errors; i++) {
field_element_t loc = field_pow(field, field.exp[error_locations[i]], generator_root_gap);
// field_element_t loc = field.exp[error_locations[i]];
error_roots[i] = field_div(field, 1, loc);
// error_roots[i] = loc;
}
}
// erasure method -- given the roots of the error locator, create the polynomial
static polynomial_t reed_solomon_find_error_locator_from_roots(field_t field, unsigned int num_errors,
field_element_t *error_roots,
polynomial_t error_locator,
polynomial_t *scratch) {
// multiply out roots to build the error locator polynomial
return polynomial_init_from_roots(field, num_errors, error_roots, error_locator, scratch);
}
// erasure method
static void reed_solomon_find_modified_syndromes(correct_reed_solomon *rs, field_element_t *syndromes, polynomial_t error_locator, field_element_t *modified_syndromes) {
polynomial_t syndrome_poly;
syndrome_poly.order = rs->min_distance - 1;
syndrome_poly.coeff = syndromes;
polynomial_t modified_syndrome_poly;
modified_syndrome_poly.order = rs->min_distance - 1;
modified_syndrome_poly.coeff = modified_syndromes;
polynomial_mul(rs->field, error_locator, syndrome_poly, modified_syndrome_poly);
}
void correct_reed_solomon_decoder_create(correct_reed_solomon *rs) {
rs->has_init_decode = true;
rs->syndromes = calloc(rs->min_distance, sizeof(field_element_t));
rs->modified_syndromes = calloc(2 * rs->min_distance, sizeof(field_element_t));
rs->received_polynomial = polynomial_create(rs->block_length - 1);
rs->error_locator = polynomial_create(rs->min_distance);
rs->error_locator_log = polynomial_create(rs->min_distance);
rs->erasure_locator = polynomial_create(rs->min_distance);
rs->error_roots = calloc(2 * rs->min_distance, sizeof(field_element_t));
rs->error_vals = malloc(rs->min_distance * sizeof(field_element_t));
rs->error_locations = malloc(rs->min_distance * sizeof(field_logarithm_t));
rs->last_error_locator = polynomial_create(rs->min_distance);
rs->error_evaluator = polynomial_create(rs->min_distance - 1);
rs->error_locator_derivative = polynomial_create(rs->min_distance - 1);
// calculate and store the first block_length powers of every generator root
// we would have to do this work in order to calculate the syndromes
// if we save it, we can prevent the need to recalculate it on subsequent calls
// total memory usage is min_distance * block_length bytes e.g. 32 * 255 ~= 8k
rs->generator_root_exp = malloc(rs->min_distance * sizeof(field_logarithm_t *));
for (unsigned int i = 0; i < rs->min_distance; i++) {
rs->generator_root_exp[i] = malloc(rs->block_length * sizeof(field_logarithm_t));
polynomial_build_exp_lut(rs->field, rs->generator_roots[i], rs->block_length - 1, rs->generator_root_exp[i]);
}
// calculate and store the first min_distance powers of every element in the field
// we would have to do this for chien search anyway, and its size is only 256 * min_distance bytes
// for min_distance = 32 this is 8k of memory, a pittance for the speedup we receive in exchange
// we also get to reuse this work during error value calculation
rs->element_exp = malloc(256 * sizeof(field_logarithm_t *));
for (field_operation_t i = 0; i < 256; i++) {
rs->element_exp[i] = malloc(rs->min_distance * sizeof(field_logarithm_t));
polynomial_build_exp_lut(rs->field, i, rs->min_distance - 1, rs->element_exp[i]);
}
rs->init_from_roots_scratch[0] = polynomial_create(rs->min_distance);
rs->init_from_roots_scratch[1] = polynomial_create(rs->min_distance);
}
ssize_t correct_reed_solomon_decode(correct_reed_solomon *rs, const uint8_t *encoded, size_t encoded_length,
uint8_t *msg) {
if (encoded_length > rs->block_length) {
return -1;
}
// the message is the non-remainder part
size_t msg_length = encoded_length - rs->min_distance;
// if they handed us a nonfull block, we'll write in 0s
size_t pad_length = rs->block_length - encoded_length;
if (!rs->has_init_decode) {
// initialize rs for decoding
correct_reed_solomon_decoder_create(rs);
}
// we need to copy to our local buffer
// the buffer we're given has the coordinates in the wrong direction
// e.g. byte 0 corresponds to the 254th order coefficient
// so we're going to flip and then write padding
// the final copied buffer will look like
// | rem (rs->min_distance) | msg (msg_length) | pad (pad_length) |
for (unsigned int i = 0; i < encoded_length; i++) {
rs->received_polynomial.coeff[i] = encoded[encoded_length - (i + 1)];
}
// fill the pad_length with 0s
for (unsigned int i = 0; i < pad_length; i++) {
rs->received_polynomial.coeff[i + encoded_length] = 0;
}
bool all_zero = reed_solomon_find_syndromes(rs->field, rs->received_polynomial, rs->generator_root_exp,
rs->syndromes, rs->min_distance);
if (all_zero) {
// syndromes were all zero, so there was no error in the message
// copy to msg and we are done
for (unsigned int i = 0; i < msg_length; i++) {
msg[i] = rs->received_polynomial.coeff[encoded_length - (i + 1)];
}
return msg_length;
}
unsigned int order = reed_solomon_find_error_locator(rs, 0);
// XXX fix this vvvv
rs->error_locator.order = order;
for (unsigned int i = 0; i <= rs->error_locator.order; i++) {
// this is a little strange since the coeffs are logs, not elements
// also, we'll be storing log(0) = 0 for any 0 coeffs in the error locator
// that would seem bad but we'll just be using this in chien search, and we'll skip all 0 coeffs
// (you might point out that log(1) also = 0, which would seem to alias. however, that's ok,
// because log(1) = 255 as well, and in fact that's how it's represented in our log table)
rs->error_locator_log.coeff[i] = rs->field.log[rs->error_locator.coeff[i]];
}
rs->error_locator_log.order = rs->error_locator.order;
if (!reed_solomon_factorize_error_locator(rs->field, 0, rs->error_locator_log, rs->error_roots, rs->element_exp)) {
// roots couldn't be found, so there were too many errors to deal with
// RS has failed for this message
return -1;
}
reed_solomon_find_error_locations(rs->field, rs->generator_root_gap, rs->error_roots, rs->error_locations,
rs->error_locator.order, 0);
reed_solomon_find_error_values(rs);
for (unsigned int i = 0; i < rs->error_locator.order; i++) {
rs->received_polynomial.coeff[rs->error_locations[i]] =
field_sub(rs->field, rs->received_polynomial.coeff[rs->error_locations[i]], rs->error_vals[i]);
}
for (unsigned int i = 0; i < msg_length; i++) {
msg[i] = rs->received_polynomial.coeff[encoded_length - (i + 1)];
}
return msg_length;
}
ssize_t correct_reed_solomon_decode_with_erasures(correct_reed_solomon *rs, const uint8_t *encoded,
size_t encoded_length, const uint8_t *erasure_locations,
size_t erasure_length, uint8_t *msg) {
if (!erasure_length) {
return correct_reed_solomon_decode(rs, encoded, encoded_length, msg);
}
if (encoded_length > rs->block_length) {
return -1;
}
if (erasure_length > rs->min_distance) {
return -1;
}
// the message is the non-remainder part
size_t msg_length = encoded_length - rs->min_distance;
// if they handed us a nonfull block, we'll write in 0s
size_t pad_length = rs->block_length - encoded_length;
if (!rs->has_init_decode) {
// initialize rs for decoding
correct_reed_solomon_decoder_create(rs);
}
// we need to copy to our local buffer
// the buffer we're given has the coordinates in the wrong direction
// e.g. byte 0 corresponds to the 254th order coefficient
// so we're going to flip and then write padding
// the final copied buffer will look like
// | rem (rs->min_distance) | msg (msg_length) | pad (pad_length) |
for (unsigned int i = 0; i < encoded_length; i++) {
rs->received_polynomial.coeff[i] = encoded[encoded_length - (i + 1)];
}
// fill the pad_length with 0s
for (unsigned int i = 0; i < pad_length; i++) {
rs->received_polynomial.coeff[i + encoded_length] = 0;
}
for (unsigned int i = 0; i < erasure_length; i++) {
// remap the coordinates of the erasures
rs->error_locations[i] = rs->block_length - (erasure_locations[i] + pad_length + 1);
}
reed_solomon_find_error_roots_from_locations(rs->field, rs->generator_root_gap, rs->error_locations,
rs->error_roots, erasure_length);
rs->erasure_locator =
reed_solomon_find_error_locator_from_roots(rs->field, erasure_length, rs->error_roots, rs->erasure_locator, rs->init_from_roots_scratch);
bool all_zero = reed_solomon_find_syndromes(rs->field, rs->received_polynomial, rs->generator_root_exp,
rs->syndromes, rs->min_distance);
if (all_zero) {
// syndromes were all zero, so there was no error in the message
// copy to msg and we are done
for (unsigned int i = 0; i < msg_length; i++) {
msg[i] = rs->received_polynomial.coeff[encoded_length - (i + 1)];
}
return msg_length;
}
reed_solomon_find_modified_syndromes(rs, rs->syndromes, rs->erasure_locator, rs->modified_syndromes);
field_element_t *syndrome_copy = malloc(rs->min_distance * sizeof(field_element_t));
memcpy(syndrome_copy, rs->syndromes, rs->min_distance * sizeof(field_element_t));
for (unsigned int i = erasure_length; i < rs->min_distance; i++) {
rs->syndromes[i - erasure_length] = rs->modified_syndromes[i];
}
unsigned int order = reed_solomon_find_error_locator(rs, erasure_length);
// XXX fix this vvvv
rs->error_locator.order = order;
for (unsigned int i = 0; i <= rs->error_locator.order; i++) {
// this is a little strange since the coeffs are logs, not elements
// also, we'll be storing log(0) = 0 for any 0 coeffs in the error locator
// that would seem bad but we'll just be using this in chien search, and we'll skip all 0 coeffs
// (you might point out that log(1) also = 0, which would seem to alias. however, that's ok,
// because log(1) = 255 as well, and in fact that's how it's represented in our log table)
rs->error_locator_log.coeff[i] = rs->field.log[rs->error_locator.coeff[i]];
}
rs->error_locator_log.order = rs->error_locator.order;
/*
for (unsigned int i = 0; i < erasure_length; i++) {
rs->error_roots[i] = field_div(rs->field, 1, rs->error_roots[i]);
}
*/
if (!reed_solomon_factorize_error_locator(rs->field, erasure_length, rs->error_locator_log, rs->error_roots, rs->element_exp)) {
// roots couldn't be found, so there were too many errors to deal with
// RS has failed for this message
free(syndrome_copy);
return -1;
}
polynomial_t temp_poly = polynomial_create(rs->error_locator.order + erasure_length);
polynomial_mul(rs->field, rs->erasure_locator, rs->error_locator, temp_poly);
polynomial_t placeholder_poly = rs->error_locator;
rs->error_locator = temp_poly;
reed_solomon_find_error_locations(rs->field, rs->generator_root_gap, rs->error_roots, rs->error_locations,
rs->error_locator.order, erasure_length);
memcpy(rs->syndromes, syndrome_copy, rs->min_distance * sizeof(field_element_t));
reed_solomon_find_error_values(rs);
for (unsigned int i = 0; i < rs->error_locator.order; i++) {
rs->received_polynomial.coeff[rs->error_locations[i]] =
field_sub(rs->field, rs->received_polynomial.coeff[rs->error_locations[i]], rs->error_vals[i]);
}
rs->error_locator = placeholder_poly;
for (unsigned int i = 0; i < msg_length; i++) {
msg[i] = rs->received_polynomial.coeff[encoded_length - (i + 1)];
}
polynomial_destroy(temp_poly);
free(syndrome_copy);
return msg_length;
}