Reed–Solomon error correction

Reed–So 4}} (t 3 x² + 2 x + 1}}, then the codeword is calculated as follows.

$s_{r}(x)=p(x)\,x^{t}\mod g(x)=547x^{3}+738x^{2}+442x+455$

$s(x)=p(x)\,x^{t}-s_{r}(x)=3x^{6}+2x^{5}+1x^{4}+382x^{3}+191x^{2}+487x+474$

Errors in transmission might cause this to be received instead.

$r(x)=s(x)+e(x)=3x^{6}+2x^{5}+123x^{4}+456x^{3}+191x^{2}+487x+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 x₁ = 757 = 3⁻³ and x₂ = 562 = 3⁻⁴, corresponding to the error locations. To calculate the error values, apply the Forney algorithm.

$\Omega (x)=S(x)\Lambda (x)\mod x^{4}=546x+732\,$

$\Lambda '(x)=658x+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 e₁x³ and e₂x⁴ 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

A_i(x) S(x) + B_i(x) x^t = R_i(x)

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

A_i(x) = Λ(x)

B_i(x) = -Q(x)

R_i(x) = Ω(x).

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

R₋₁ = x^t

R₀ = S(x)

A₋₁ = 0

A₀ = 1

i = 0

while degree of R_i >= t/2

i = i + 1

Q = R_i-2 / R_i-1

R_i = R_i-2 - Q R_i-1

A_i = A_i-2 - Q A_i-1

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

Λ(x) = A_i / A_i(0)

Ω(x) = R_i / A_i(0)

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

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) = A₂ / 544 = 329 x² + 821 x + 001

Ω(x) = R₂ / 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\leq j\leq 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

% Add leading zeros
function xpad = pad(x,k)
 len = length(x);
 if (len<k)
 xpad = [zeros(1, k-len) x];
 end
end

See also

• BCH code
• Chien search
• Berlekamp–Massey algorithm
• Forward error correction
• Berlekamp–Welch algorithm
• Folded Reed–Solomon code

Source

http://wikipedia.org/