Reed–Solomon error correction
Reed–So 4}} (t 3 x2 + 2 x + 1}}, then the codeword is calculated as follows.
Errors in transmission might cause this to be received instead.
The syndromes are calculated by evaluating r at powers of α.
To correct the errors, first use the Berlekamp–Massey algorithm to calculate the error locator polynomial.
|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.
Subtracting e1x3 and e2x4 from the received polynomial r reproduces the original codeword s.
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:
The key equation is:
For t = 6 and e = 3:
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.
Using the same data as the Berlekamp Massey example above:
|-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:
Use through 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 to , the error locator polynomial, and these formulas
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 and its extensions.
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.
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
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 on BitcoinWiki
- BCH code
- Chien search
- Berlekamp–Massey algorithm
- Forward error correction
- Berlekamp–Welch algorithm
- Folded Reed–Solomon code