Reed–Solomon error correction

This is the approved revision of this page, as well as being the most recent.

Reed–So 4}} (t 3 x2 + 2 x + 1}}, then the codeword is calculated as follows. $s_r(x) = p(x) \, x^t \mod g(x) = 547 x^3 + 738 x^2 + 442 x + 455$ $s(x) = p(x) \, x^t - s_r(x) = 3 x^6 + 2 x^5 + 1 x^4 + 382 x^3 + 191 x^2 + 487 x + 474$

Errors in transmission might cause this to be received instead. $r(x) = s(x) + e(x) = 3 x^6 + 2 x^5 + 123 x^4 + 456 x^3 + 191 x^2 + 487 x + 474$

The syndromes are calculated by evaluating r at powers of α. $S_1 = r(3^1) = 3\cdot 3^6 + 2\cdot 3^5 + 123\cdot 3^4 + 456\cdot 3^3 + 191\cdot 3^2 + 487\cdot 3 + 474 = 732$ $S_2 = r(3^2) = 637,\;S_3 = r(3^3) = 762,\;S_4 = r(3^4) = 925$

To correct the errors, first use the Berlekamp–Massey algorithm to calculate the error locator polynomial.

n Sn+1 d C B b m
0 732 732 197 x + 1 1 732 1
1 637 846 173 x + 1 1 732 2
2 762 412 634 x2 + 173 x + 1 173 x + 1 412 1
3 925 576 329 x2 + 821 x + 1 173 x + 1 412 2

The final value of C is the error locator polynomial, Λ(x). The zeros can be found by trial substitution. They are x1 = 757 = 3−3 and x2 = 562 = 3−4, corresponding to the error locations. To calculate the error values, apply the Forney algorithm. $\Omega(x) = S(x) \Lambda(x) \mod x^4 = 546 x + 732\,$ $\Lambda'(x) = 658 x + 821\,$ $e_1 = -\Omega(x_1)/\Lambda'(x_1) = -649/54 = 280 \times 843 = 74\,$ $e_2 = -\Omega(x_2)/\Lambda'(x_2) = 122\,$

Subtracting e1x3 and e2x4 from the received polynomial r reproduces the original codeword s.

Euclidean decoder

Another iterative method for calculating both the error locator polynomial and the error value polynomial is based on Sugiyama's adaptation of the Extended Euclidean algorithm .

Define S(x), Λ(x), and Ω(x) for t syndromes and e errors: $S(x) = S_{t} x^{t-1} + S_{t-1} x^{t-2} + \cdots + S_2 x + S_1$ $\Lambda(x) = \Lambda_{e} x^{e} + \Lambda_{e-1} x^{e-1} + \cdots + \Lambda_{1} x + 1$ $\Omega(x) = \Omega_{e} x^{e} + \Omega_{e-1} x^{e-1} + \cdots + \Omega_{1} x + \Omega_{0}$

The key equation is: $\Lambda(x) S(x) = Q(x) x^{t} + \Omega(x)$

For t = 6 and e = 3: $\begin{bmatrix} \Lambda_3 S_6 & x^8 \\ \Lambda_2 S_6 + \Lambda_3 S_5 & x^7 \\ \Lambda_1 S_6 + \Lambda_2 S_5 + \Lambda_3 S_4 & x^6 \\ S_6 + \Lambda_1 S_5 + \Lambda_2 S_4 + \Lambda_3 S_3 & x^5 \\ S_5 + \Lambda_1 S_4 + \Lambda_2 S_3 + \Lambda_3 S_2 & x^4 \\ S_4 + \Lambda_1 S_3 + \Lambda_2 S_2 + \Lambda_3 S_1 & x^3 \\ S_3 + \Lambda_1 S_2 + \Lambda_2 S_1 & x^2 \\ S_2 + \Lambda_1 S_1 & x \\ S_1 & \\ \end{bmatrix} = \begin{bmatrix} Q_2 x^8 \\ Q_1 x^7 \\ Q_0 x^6 \\ 0 \\ 0 \\ 0 \\ \Omega_2 x^2 \\ \Omega_1 x \\ \Omega_0 \\ \end{bmatrix}$

The middle terms are zero due to the relationship between Λ and syndromes.

The extended Euclidean algorithm can find a series of polynomials of the form

Ai(x) S(x) + Bi(x) xt = Ri(x)

where the degree of R decreases as i increases. Once the degree of Ri(x) < t/2, then

Ai(x) = Λ(x)

Bi(x) = -Q(x)

Ri(x) = Ω(x).

B(x) and Q(x) don't need to be saved, so the algorithm becomes:

R−1 = xt
R0 = S(x)
A−1 = 0
A0 = 1
i = 0
while degree of Ri >= t/2
i = i + 1
Q = Ri-2 / Ri-1
Ri = Ri-2 - Q Ri-1
Ai = Ai-2 - Q Ai-1

to set low order term of Λ(x) to 1, divide Λ(x) and Ω(x) by Ai(0):

Λ(x) = Ai / Ai(0)
Ω(x) = Ri / Ai(0)

Ai(0) is the constant (low order) term of Ai.

Example

Using the same data as the Berlekamp Massey example above:

i Ri Ai
-1 001 x4 + 000 x3 + 000 x2 + 000 x + 000 000
0 925 x3 + 762 x2 + 637 x + 732 001
1 683 x2 + 676 x + 024 697 x + 396
2 673 x + 596 608 x2 + 704 x + 544
Λ(x) = A2 / 544 = 329 x2 + 821 x + 001
Ω(x) = R2 / 544 = 546 x + 732

Decoder using discrete Fourier transform

A discrete Fourier transform can be used for decoding. To avoid conflict with syndrome names, let c(x) = s(x) the encoded codeword. r(x) and e(x) are the same as above. Define C(x), E(x), and R(x) as the discrete Fourier transforms of c(x), e(x), and r(x). Since r(x) = c(x) + e(x), and since a discrete Fourier transform is a linear operator, R(x) = C(x) + E(x).

Transform r(x) to R(x) using discrete Fourier transform. Since the calculation for a discrete Fourier transform is the same as the calculation for syndromes, t coefficients of R(x) and E(x) are the same as the syndromes: $R_j = E_j = S_j = r(\alpha^j)$ $for \ 1 \le j \le t$

Use $R_1$ through $R_t$ as syndromes (they're the same) and generate the error locator polynomial using the methods from any of the above decoders.

Let v = number of errors. Generate E(x) using the known coefficients $E_1$ to $E_t$, the error locator polynomial, and these formulas $E_0 = - \frac{1}{\sigma_v}(E_{v} + \sigma_1 E_{v-1} + \cdots + \sigma_{v-1} E_{1})$ $E_j = -(\sigma_1 E_{j-1} + \sigma_2 E_{j-2} + \cdots + \sigma_v E_{j-v})$ $for \ t < j < n$

Then calculate C(x) = R(x) - E(x) and take the inverse transform of C(x) to produce c(x).

Decoding beyond the error-correction bound

The Singleton bound states that the minimum distance d of a linear block code of size (n,k) is upper-bounded by n − k + 1. The distance d was usually understood to limit the error-correction capability to ⌊d/2⌋. The Reed–Solomon code achieves this bound with equality, and can thus correct up to ⌊(n − k + 1)/2⌋ errors. However, this error-correction bound is not exact.

In 1999, Madhu Sudan and Venkatesan Guruswami at MIT published "Improved Decoding of Reed–Solomon and Algebraic-Geometry Codes" introducing an algorithm that allowed for the correction of errors beyond half the minimum distance of the code. It applies to Reed–Solomon codes and more generally to algebraic geometric codes. This algorithm produces a list of codewords (it is a list-decoding algorithm) and is based on interpolation and factorization of polynomials over $GF(2^m)$ and its extensions.

Soft-decoding

The algebraic decoding methods described above are hard-decision methods, which means that for every symbol a hard decision is made about its value. For example, a decoder could associate with each symbol an additional value corresponding to the channel demodulator's confidence in the correctness of the symbol. The advent of LDPC and turbo codes, which employ iterated soft-decision belief propagation decoding methods to achieve error-correction performance close to the theoretical limit, has spurred interest in applying soft-decision decoding to conventional algebraic codes. In 2003, Ralf Koetter and Alexander Vardy presented a polynomial-time soft-decision algebraic list-decoding algorithm for Reed–Solomon codes, which was based upon the work by Sudan and Guruswami. In 2016, Steven J. Franke and Joseph H. Taylor published a novel soft-decision decoder.

Matlab Example

Encoder

Here we present a simple Matlab implementation for an encoder.

function [ encoded ] = rsEncoder( msg, m, prim_poly, n, k )
%RSENCODER Encode message with the Reed-Solomon algorithm
% m is the number of bits per symbol
% prim_poly: Primitive polynomial p(x). Ie for DM is 301
% k is the size of the message
% n is the total size (k+redundant)
% Example: msg = uint8('Test')
% enc_msg = rsEncoder(msg, 8, 301, 12, numel(msg));

% Get the alpha
alpha = gf(2, m, prim_poly);

% Get the Reed-Solomon generating polynomial g(x)
g_x = genpoly(k, n, alpha);

% Multiply the information by X^(n-k), or just pad with zeros at the end to
% get space to add the redundant information
msg_padded = gf([msg zeros(1, n-k)], m, prim_poly);

% Get the remainder of the division of the extended message by the
% Reed-Solomon generating polynomial g(x)
[~, reminder] = deconv(msg_padded, g_x);

% Now return the message with the redundant information
encoded = msg_padded - reminder;

end

% Find the Reed-Solomon generating polynomial g(x), by the way this is the
% same as the rsgenpoly function on matlab
function g = genpoly(k, n, alpha)
g = 1;
% A multiplication on the galois field is just a convolution
for k = mod(1 : n-k, n)
g = conv(g, [1 alpha .^ (k)]);
end
end

Decoder

Now the decoding part:

function [ decoded, error_pos, error_mag, g, S ] = rsDecoder( encoded, m, prim_poly, n, k )
%RSDECODER Decode a Reed-Solomon encoded message
% Example:
% [dec, ~, ~, ~, ~] = rsDecoder(enc_msg, 8, 301, 12, numel(msg))
max_errors = floor((n-k)/2);
orig_vals = encoded.x;
% Initialize the error vector
errors = zeros(1, n);
g = [];
S = [];

% Get the alpha
alpha = gf(2, m, prim_poly);

% Find the syndromes (Check if dividing the message by the generator
% polynomial the result is zero)
Synd = polyval(encoded, alpha .^ (1:n-k));
Syndromes = trim(Synd);

% If all syndromes are zeros (perfectly divisible) there are no errors
if isempty(Syndromes.x)
decoded = orig_vals(1:k);
error_pos = [];
error_mag = [];
g = [];
S = Synd;
return;
end

% Prepare for the euclidean algorithm (Used to find the error locating
% polynomials)
r0 = [1, zeros(1, 2*max_errors)]; r0 = gf(r0, m, prim_poly); r0 = trim(r0);
size_r0 = length(r0);
r1 = Syndromes;
f0 = gf([zeros(1, size_r0-1) 1], m, prim_poly);
f1 = gf(zeros(1, size_r0), m, prim_poly);
g0 = f1; g1 = f0;

% Do the euclidean algorithm on the polynomials r0(x) and Syndromes(x) in
% order to find the error locating polynomial
while true
% Do a long division
[quotient, remainder] = deconv(r0, r1);
% Add some zeros
quotient = pad(quotient, length(g1));

% Find quotient*g1 and pad
c = conv(quotient, g1);
c = trim(c);
c = pad(c, length(g0));

% Update g as g0-quotient*g1
g = g0 - c;

% Check if the degree of remainder(x) is less than max_errors
if all(remainder(1:end - max_errors) == 0)
break;
end

% Update r0, r1, g0, g1 and remove leading zeros
r0 = trim(r1); r1 = trim(remainder);
g0 = g1; g1 = g;
end

% Remove leading zeros
g = trim(g);

% Find the zeros of the error polynomial on this galois field
evalPoly = polyval(g, alpha .^ (n-1 : -1 : 0));
error_pos = gf(find(evalPoly == 0), m);

% If no error position is found we return the received work, because
% basically is nothing that we could do and we return the received message
if isempty(error_pos)
decoded = orig_vals(1:k);
error_mag = [];
return;
end

% Prepare a linear system to solve the error polynomial and find the error
% magnitudes
size_error = length(error_pos);
Syndrome_Vals = Syndromes.x;
b(:, 1) = Syndrome_Vals(1:size_error);
for idx = 1 : size_error
e = alpha .^ (idx*(n-error_pos.x));
err = e.x;
er(idx, :) = err;
end

% Solve the linear system
error_mag = (gf(er, m, prim_poly) \ gf(b, m, prim_poly))';
% Put the error magnitude on the error vector
errors(error_pos.x) = error_mag.x;
% Bring this vector to the galois field
errors_gf = gf(errors, m, prim_poly);

% Now to fix the errors just add with the encoded code
decoded_gf = encoded(1:k) + errors_gf(1:k);
decoded = decoded_gf.x;

end

% Remove leading zeros from galois array
function gt = trim(g)
gx = g.x;
gt = gf(gx(find(gx, 1) : end), g.m, g.prim_poly);
end