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) = 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) = 3cdot 3^6 + 2cdot 3^5 + 123cdot 3^4 + 456cdot 3^3 + 191cdot 3^2 + 487cdot 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	S_n+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 x² + 173 x + 1	173 x + 1	412	1
3	925	576	329 x² + 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 = 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 e₁x³ and e₂x⁴ from the received polynomial r reproduces the original codeword s.

1 Euclidean decoder
- 1.1 Example
2 Decoder using discrete Fourier transform
3 Decoding beyond the error-correction bound
4 Soft-decoding
5 Matlab Example
- 5.1 Encoder
- 5.2 Decoder
6 See also
7 Source

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 .

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	R_i	A_i
-1	001 x⁴ + 000 x³ + 000 x² + 000 x + 000	000
0	925 x³ + 762 x² + 637 x + 732	001
1	683 x² + 676 x + 024	697 x + 396
2	673 x + 596	608 x² + 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 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 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, and 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 . This algorithm produces a list of codewords (it is a 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 ‘s confidence in the correctness of the symbol. The advent of LDPC and turbo codes, which employ iterated belief propagation decoding methods to achieve error-correction performance close to the , has spurred interest in applying soft-decision decoding to conventional algebraic codes. In 2003, Ralf Koetter and 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

Source

http://wikipedia.org/