python - Errata(borrones+errores) Berlekamp-Massey para la decodificación Reed-Solomon
math error-correction (2)
Después de leer muchos artículos de investigación y libros, el único lugar donde encontré la respuesta está en el libro ( legible en línea en Google Books , pero no disponible como PDF):
"Códigos algebraicos para la transmisión de datos", Blahut, Richard E., 2003, editorial universitaria de Cambridge.
Aquí hay algunos extractos de este libro, que detalla la descripción exacta (a excepción de la representación matricial / vectorizada de las operaciones polinómicas) del algoritmo Berlekamp-Massey que implementé:
Y aquí está el algoritmo de errata (errores y borraduras) Berlekamp-Massey para Reed-Solomon:
Como puede ver, contrariamente a la descripción habitual de que solo necesita inicializar Lambda, el polinomio del localizador de errores, con el valor del polinomio del localizador de borrados calculado previamente , también debe omitir las primeras iteraciones, donde v es la cantidad de borraduras Tenga en cuenta que no es equivalente a omitir las últimas v iteraciones: debe omitir las primeras v iteraciones , porque r (el contador de iteraciones, K en mi implementación) se usa no solo para contar iteraciones sino también para generar el factor de discrepancia Delta correcto.
Aquí está el código resultante con las modificaciones para admitir borrones, así como errores hasta v+2*e <= (nk) :
def _berlekamp_massey(self, s, k=None, erasures_loc=None, erasures_eval=None, erasures_count=0):
''''''Computes and returns the errata (errors+erasures) locator polynomial (sigma) and the
error evaluator polynomial (omega) at the same time.
If the erasures locator is specified, we will return an errors-and-erasures locator polynomial and an errors-and-erasures evaluator polynomial, else it will compute only errors. With erasures in addition to errors, it can simultaneously decode up to v+2e <= (n-k) where v is the number of erasures and e the number of errors.
Mathematically speaking, this is equivalent to a spectral analysis (see Blahut, "Algebraic Codes for Data Transmission", 2003, chapter 7.6 Decoding in Time Domain).
The parameter s is the syndrome polynomial (syndromes encoded in a
generator function) as returned by _syndromes.
Notes:
The error polynomial:
E(x) = E_0 + E_1 x + ... + E_(n-1) x^(n-1)
j_1, j_2, ..., j_s are the error positions. (There are at most s
errors)
Error location X_i is defined: X_i = α^(j_i)
that is, the power of α (alpha) corresponding to the error location
Error magnitude Y_i is defined: E_(j_i)
that is, the coefficient in the error polynomial at position j_i
Error locator polynomial:
sigma(z) = Product( 1 - X_i * z, i=1..s )
roots are the reciprocals of the error locations
( 1/X_1, 1/X_2, ...)
Error evaluator polynomial omega(z) is here computed at the same time as sigma, but it can also be constructed afterwards using the syndrome and sigma (see _find_error_evaluator() method).
It can be seen that the algorithm tries to iteratively solve for the error locator polynomial by
solving one equation after another and updating the error locator polynomial. If it turns out that it
cannot solve the equation at some step, then it computes the error and weights it by the last
non-zero discriminant found, and delays the weighted result to increase the polynomial degree
by 1. Ref: "Reed Solomon Decoder: TMS320C64x Implementation" by Jagadeesh Sankaran, December 2000, Application Report SPRA686
The best paper I found describing the BM algorithm for errata (errors-and-erasures) evaluator computation is in "Algebraic Codes for Data Transmission", Richard E. Blahut, 2003.
''''''
# For errors-and-erasures decoding, see: "Algebraic Codes for Data Transmission", Richard E. Blahut, 2003 and (but it''s less complete): Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
# also see: Blahut, Richard E. "A universal Reed-Solomon decoder." IBM Journal of Research and Development 28.2 (1984): 150-158. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2084&rep=rep1&type=pdf
# and another good alternative book with concrete programming examples: Jiang, Yuan. A practical guide to error-control coding using Matlab. Artech House, 2010.
n = self.n
if not k: k = self.k
# Initialize, depending on if we include erasures or not:
if erasures_loc:
sigma = [ Polynomial(erasures_loc.coefficients) ] # copy erasures_loc by creating a new Polynomial, so that we initialize the errata locator polynomial with the erasures locator polynomial.
B = [ Polynomial(erasures_loc.coefficients) ]
omega = [ Polynomial(erasures_eval.coefficients) ] # to compute omega (the evaluator polynomial) at the same time, we also need to initialize it with the partial erasures evaluator polynomial
A = [ Polynomial(erasures_eval.coefficients) ] # TODO: fix the initial value of the evaluator support polynomial, because currently the final omega is not correct (it contains higher order terms that should be removed by the end of BM)
else:
sigma = [ Polynomial([GF256int(1)]) ] # error locator polynomial. Also called Lambda in other notations.
B = [ Polynomial([GF256int(1)]) ] # this is the error locator support/secondary polynomial, which is a funky way to say that it''s just a temporary variable that will help us construct sigma, the error locator polynomial
omega = [ Polynomial([GF256int(1)]) ] # error evaluator polynomial. We don''t need to initialize it with erasures_loc, it will still work, because Delta is computed using sigma, which itself is correctly initialized with erasures if needed.
A = [ Polynomial([GF256int(0)]) ] # this is the error evaluator support/secondary polynomial, to help us construct omega
L = [ 0 ] # update flag: necessary variable to check when updating is necessary and to check bounds (to avoid wrongly eliminating the higher order terms). For more infos, see https://www.cs.duke.edu/courses/spring11/cps296.3/decoding_rs.pdf
M = [ 0 ] # optional variable to check bounds (so that we do not mistakenly overwrite the higher order terms). This is not necessary, it''s only an additional safe check. For more infos, see the presentation decoding_rs.pdf by Andrew Brown in the doc folder.
# Fix the syndrome shifting: when computing the syndrome, some implementations may prepend a 0 coefficient for the lowest degree term (the constant). This is a case of syndrome shifting, thus the syndrome will be bigger than the number of ecc symbols (I don''t know what purpose serves this shifting). If that''s the case, then we need to account for the syndrome shifting when we use the syndrome such as inside BM, by skipping those prepended coefficients.
# Another way to detect the shifting is to detect the 0 coefficients: by definition, a syndrome does not contain any 0 coefficient (except if there are no errors/erasures, in this case they are all 0). This however doesn''t work with the modified Forney syndrome (that we do not use in this lib but it may be implemented in the future), which set to 0 the coefficients corresponding to erasures, leaving only the coefficients corresponding to errors.
synd_shift = 0
if len(s) > (n-k): synd_shift = len(s) - (n-k)
# Polynomial constants:
ONE = Polynomial(z0=GF256int(1))
ZERO = Polynomial(z0=GF256int(0))
Z = Polynomial(z1=GF256int(1)) # used to shift polynomials, simply multiply your poly * Z to shift
# Precaching
s2 = ONE + s
# Iteratively compute the polynomials n-k-erasures_count times. The last ones will be correct (since the algorithm refines the error/errata locator polynomial iteratively depending on the discrepancy, which is kind of a difference-from-correctness measure).
for l in xrange(0, n-k-erasures_count): # skip the first erasures_count iterations because we already computed the partial errata locator polynomial (by initializing with the erasures locator polynomial)
K = erasures_count+l+synd_shift # skip the FIRST erasures_count iterations (not the last iterations, that''s very important!)
# Goal for each iteration: Compute sigma[l+1] and omega[l+1] such that
# (1 + s)*sigma[l] == omega[l] in mod z^(K)
# For this particular loop iteration, we have sigma[l] and omega[l],
# and are computing sigma[l+1] and omega[l+1]
# First find Delta, the non-zero coefficient of z^(K) in
# (1 + s) * sigma[l]
# Note that adding 1 to the syndrome s is not really necessary, you can do as well without.
# This delta is valid for l (this iteration) only
Delta = ( s2 * sigma[l] ).get_coefficient(K) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial).
# Make it a polynomial of degree 0, just for ease of computation with polynomials sigma and omega.
Delta = Polynomial(x0=Delta)
# Can now compute sigma[l+1] and omega[l+1] from
# sigma[l], omega[l], B[l], A[l], and Delta
sigma.append( sigma[l] - Delta * Z * B[l] )
omega.append( omega[l] - Delta * Z * A[l] )
# Now compute the next support polynomials B and A
# There are two ways to do this
# This is based on a messy case analysis on the degrees of the four polynomials sigma, omega, A and B in order to minimize the degrees of A and B. For more infos, see https://www.cs.duke.edu/courses/spring10/cps296.3/decoding_rs_scribe.pdf
# In fact it ensures that the degree of the final polynomials aren''t too large.
if Delta == ZERO or 2*L[l] > K+erasures_count /
or (2*L[l] == K+erasures_count and M[l] == 0):
#if Delta == ZERO or len(sigma[l+1]) <= len(sigma[l]): # another way to compute when to update, and it doesn''t require to maintain the update flag L
# Rule A
B.append( Z * B[l] )
A.append( Z * A[l] )
L.append( L[l] )
M.append( M[l] )
elif (Delta != ZERO and 2*L[l] < K+erasures_count) /
or (2*L[l] == K+erasures_count and M[l] != 0):
# elif Delta != ZERO and len(sigma[l+1]) > len(sigma[l]): # another way to compute when to update, and it doesn''t require to maintain the update flag L
# Rule B
B.append( sigma[l] // Delta )
A.append( omega[l] // Delta )
L.append( K - L[l] ) # the update flag L is tricky: in Blahut''s schema, it''s mandatory to use `L = K - L - erasures_count` (and indeed in a previous draft of this function, if you forgot to do `- erasures_count` it would lead to correcting only 2*(errors+erasures) <= (n-k) instead of 2*errors+erasures <= (n-k)), but in this latest draft, this will lead to a wrong decoding in some cases where it should correctly decode! Thus you should try with and without `- erasures_count` to update L on your own implementation and see which one works OK without producing wrong decoding failures.
M.append( 1 - M[l] )
else:
raise Exception("Code shouldn''t have gotten here")
# Hack to fix the simultaneous computation of omega, the errata evaluator polynomial: because A (the errata evaluator support polynomial) is not correctly initialized (I could not find any info in academic papers). So at the end, we get the correct errata evaluator polynomial omega + some higher order terms that should not be present, but since we know that sigma is always correct and the maximum degree should be the same as omega, we can fix omega by truncating too high order terms.
if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])
# Return the last result of the iterations (since BM compute iteratively, the last iteration being correct - it may already be before, but we''re not sure)
return sigma[-1], omega[-1]
def _find_erasures_locator(self, erasures_pos):
''''''Compute the erasures locator polynomial from the erasures positions (the positions must be relative to the x coefficient, eg: "hello worldxxxxxxxxx" is tampered to "h_ll_ worldxxxxxxxxx" with xxxxxxxxx being the ecc of length n-k=9, here the string positions are [1, 4], but the coefficients are reversed since the ecc characters are placed as the first coefficients of the polynomial, thus the coefficients of the erased characters are n-1 - [1, 4] = [18, 15] = erasures_loc to be specified as an argument.''''''
# See: http://ocw.usu.edu/Electrical_and_Computer_Engineering/Error_Control_Coding/lecture7.pdf and Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
erasures_loc = Polynomial([GF256int(1)]) # just to init because we will multiply, so it must be 1 so that the multiplication starts correctly without nulling any term
# erasures_loc is very simple to compute: erasures_loc = prod(1 - x*alpha[j]**i) for i in erasures_pos and where alpha is the alpha chosen to evaluate polynomials (here in this library it''s gf(3)). To generate c*x where c is a constant, we simply generate a Polynomial([c, 0]) where 0 is the constant and c is positionned to be the coefficient for x^1. See https://en.wikipedia.org/wiki/Forney_algorithm#Erasures
for i in erasures_pos:
erasures_loc = erasures_loc * (Polynomial([GF256int(1)]) - Polynomial([GF256int(self.generator)**i, 0]))
return erasures_loc
Nota : Sigma, Omega, A, B, L y M son listas de polinomios (por lo que guardamos toda la historia de todos los polinomios intermedios que calculamos en cada iteración). Por supuesto, esto puede optimizarse porque solo necesitamos Sigma[l] , Sigma[l-1] , Omega[l] , Omega[l-1] , A[l] , B[l] , L[l] y M[l] (así que solo Sigma y Omega necesitan mantener la iteración anterior en la memoria, las otras variables no son necesarias).
Nota 2 : el indicador de actualización L es complicado: en algunas implementaciones, hacer exactamente lo mismo en el esquema de Blahut dará lugar a errores incorrectos al decodificar. En mi implementación anterior, era obligatorio utilizar L = K - L - erasures_count para decodificar correctamente los errores y borrados hasta el límite Singleton, pero en mi última implementación, tuve que usar L = K - L (incluso cuando hay borraduras) para evitar fallas de descodificación erróneas. Debería probar ambos en su propia implementación y ver cuál no produce fallas de decodificación erróneas. Vea a continuación en los números para obtener más información.
El único problema con este algoritmo es que no describe cómo calcular simultáneamente Omega, el polinomio del evaluador de errores (el libro describe cómo inicializar Omega sólo para errores, pero no al descodificar errores y borrados). Intenté varias variaciones y lo anterior funciona, pero no del todo: al final, Omega incluirá términos de orden superior que deberían haberse cancelado. Probablemente Omega o A, el polinomio de soporte del evaluador de errores, no se inicializa con el buen valor.
Sin embargo, puede solucionarlo recortando el polinomio de Omega de los términos de orden demasiado altos (ya que debe tener el mismo grado que Lambda / Sigma):
if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])
O puede calcular totalmente Omega desde cero después de BM utilizando el localizador de erratas Lambda / Sigma, que siempre se calcula correctamente:
def _find_error_evaluator(self, synd, sigma, k=None):
''''''Compute the error (or erasures if you supply sigma=erasures locator polynomial) evaluator polynomial Omega from the syndrome and the error/erasures/errata locator Sigma. Omega is already computed at the same time as Sigma inside the Berlekamp-Massey implemented above, but in case you modify Sigma, you can recompute Omega afterwards using this method, or just ensure that Omega computed by BM is correct given Sigma (as long as syndrome and sigma are correct, omega will be correct).''''''
n = self.n
if not k: k = self.k
# Omega(x) = [ Synd(x) * Error_loc(x) ] mod x^(n-k+1) -- From Blahut, Algebraic codes for data transmission, 2003
return (synd * sigma) % Polynomial([GF256int(1)] + [GF256int(0)] * (n-k+1)) # Note that you should NOT do (1+Synd(x)) as can be seen in some books because this won''t work with all primitive generators.
Estoy buscando una mejor solución en la siguiente pregunta sobre CSTheory .
/ EDIT: describiré algunos de los problemas que he tenido y cómo solucionarlos:
- no olvide iniciar el polinomio del localizador de errores con el polinomio localizador de borrados (que puede calcular fácilmente a partir de las posiciones de los síndromes y borrados).
- si solo puede decodificar errores y borrados sin problemas, pero limitado a
2*errors + erasures <= (nk)/2, entonces olvidó omitir las primeras v iteraciones. - si puede decodificar tanto borrones-y-errores como hasta
2*(errors+erasures) <= (nk), entonces olvidó actualizar la asignación de L:L = i+1 - L - erasures_countdeL = i+1 - L - erasures_countlugar deL = i+1 - LPero esto puede hacer que su decodificador falle en algunos casos, dependiendo de cómo implemente su decodificador, vea el siguiente punto. - mi primer decodificador estaba limitado a un solo generador / polinomio principal / fcr, pero cuando lo actualicé para que fuera universal y agregué pruebas unitarias estrictas, el decodificador falló cuando no debería. Parece que el esquema anterior de Blahut es incorrecto sobre L (el indicador de actualización): debe actualizarse usando
L = K - Ly noL = K - L - erasures_count, porque esto provocará que el decodificador falle a veces incluso cuando estamos bajo el Singleton obligado (¡y así deberíamos decodificar correctamente!). Esto parece confirmarse por el hecho de que la computaciónL = K - Lno solo reparará esos problemas de decodificación, sino que también dará exactamente el mismo resultado que la forma alternativa de actualizar sin usar la bandera de actualización L (es decir, la condiciónif Delta == ZERO or len(sigma[l+1]) <= len(sigma[l]):. Pero esto es extraño: en mi implementación anterior,L = K - L - erasures_countera obligatorio para la decodificación de errores y borrados, pero ahora parece que produce errores incorrectos. Entonces debería probar con y sin su propia implementación y si uno u otro produce fallas incorrectas para usted. - tenga en cuenta que la condición
2*L[l] > Kcambia a2*L[l] > K+erasures_count. Es posible que no note ningún efecto secundario sin agregar la condición+erasures_countal principio, pero en algunos casos la decodificación fallará cuando no sea así. - si puede corregir solo un error o una eliminación, verifique que su condición sea
2*L[l] > K+erasures_county no2*L[l] >= K+erasures_count(observe el>lugar de>=). - si puede corregir
2*errors + erasures <= (nk-2)(justo debajo del límite, por ejemplo, si tiene 10 símbolos ecc, puede corregir solo 4 errores en lugar de 5 normalmente) luego verifique su síndrome y su ciclo en el interior el BM algo: si el síndrome comienza con un coeficiente 0 para el término constante x ^ 0 (que a veces se recomienda en los libros), entonces su síndrome se desplaza, y luego su ciclo dentro de BM debe comenzar en1y terminar enn-k+1lugar de0:(nk)si no se desplaza. - Si puede corregir cada símbolo, excepto el último (el último símbolo ecc), luego verifique sus rangos, particularmente en su búsqueda Chien: no debe evaluar el polinomio localizador de errores de alfa ^ 0 a alfa ^ 255, sino de alfa ^ 1 a alpha ^ 256.
Estoy tratando de implementar un codificador-decodificador Reed-Solomon en Python que soporte la decodificación de borrados y errores, y eso me está volviendo loco.
En la actualidad, la implementación solo admite la decodificación de errores o solo borrados, pero no ambos al mismo tiempo (incluso si está por debajo del límite teórico de 2 * errores + borrados <= (nk)).
De los documentos de Blahut ( aquí y aquí ), parece que solo necesitamos inicializar el polinomio localizador de errores con el polinomio localizador de borrado para calcular implícitamente el polinomio localizador de erratas dentro de Berlekamp-Massey.
Este enfoque funciona parcialmente para mí: cuando tengo 2 * errores + borrones <(nk) / 2 funciona, pero de hecho después de la depuración solo funciona porque BM calcula un polinomio localizador de errores que obtiene exactamente el mismo valor que el polinomio localizador de borrado (porque estamos por debajo del límite para la corrección de solo errores), y por lo tanto se trunca a través de campos de galois y terminamos con el valor correcto del polinomio localizador de borrado (al menos así es como lo entiendo, puedo estar equivocado).
Sin embargo, cuando vamos por encima (nk) / 2, por ejemplo si n = 20 yk = 11, entonces tenemos (nk) = 9 símbolos borrados que podemos corregir, si alimentamos en 5 borrados, entonces BM simplemente sale mal. Si alimentamos en 4 borrados + 1 error (todavía estamos muy por debajo del límite ya que tenemos 2 * errores + borrados = 2 + 4 = 6 <9), el BM sigue saliendo mal.
El algoritmo exacto de Berlekamp-Massey I implementado se puede encontrar en esta presentación (páginas 15-17), pero una descripción muy similar se puede encontrar aquí y aquí , y aquí adjunto una copia de la descripción matemática:
Ahora, tengo una reproducción casi exacta de este algoritmo matemático en un código de Python. Lo que me gustaría es ampliarlo para admitir borrones, lo cual intenté inicializando el localizador de errores sigma con el localizador de borrado:
def _berlekamp_massey(self, s, k=None, erasures_loc=None):
''''''Computes and returns the error locator polynomial (sigma) and the
error evaluator polynomial (omega).
If the erasures locator is specified, we will return an errors-and-erasures locator polynomial and an errors-and-erasures evaluator polynomial.
The parameter s is the syndrome polynomial (syndromes encoded in a
generator function) as returned by _syndromes. Don''t be confused with
the other s = (n-k)/2
Notes:
The error polynomial:
E(x) = E_0 + E_1 x + ... + E_(n-1) x^(n-1)
j_1, j_2, ..., j_s are the error positions. (There are at most s
errors)
Error location X_i is defined: X_i = a^(j_i)
that is, the power of a corresponding to the error location
Error magnitude Y_i is defined: E_(j_i)
that is, the coefficient in the error polynomial at position j_i
Error locator polynomial:
sigma(z) = Product( 1 - X_i * z, i=1..s )
roots are the reciprocals of the error locations
( 1/X_1, 1/X_2, ...)
Error evaluator polynomial omega(z) is here computed at the same time as sigma, but it can also be constructed afterwards using the syndrome and sigma (see _find_error_evaluator() method).
''''''
# For errors-and-erasures decoding, see: Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
# also see: Blahut, Richard E. "A universal Reed-Solomon decoder." IBM Journal of Research and Development 28.2 (1984): 150-158. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2084&rep=rep1&type=pdf
# or alternatively see the reference book by Blahut: Blahut, Richard E. Theory and practice of error control codes. Addison-Wesley, 1983.
# and another good alternative book with concrete programming examples: Jiang, Yuan. A practical guide to error-control coding using Matlab. Artech House, 2010.
n = self.n
if not k: k = self.k
# Initialize:
if erasures_loc:
sigma = [ Polynomial(erasures_loc.coefficients) ] # copy erasures_loc by creating a new Polynomial
B = [ Polynomial(erasures_loc.coefficients) ]
else:
sigma = [ Polynomial([GF256int(1)]) ] # error locator polynomial. Also called Lambda in other notations.
B = [ Polynomial([GF256int(1)]) ] # this is the error locator support/secondary polynomial, which is a funky way to say that it''s just a temporary variable that will help us construct sigma, the error locator polynomial
omega = [ Polynomial([GF256int(1)]) ] # error evaluator polynomial. We don''t need to initialize it with erasures_loc, it will still work, because Delta is computed using sigma, which itself is correctly initialized with erasures if needed.
A = [ Polynomial([GF256int(0)]) ] # this is the error evaluator support/secondary polynomial, to help us construct omega
L = [ 0 ] # necessary variable to check bounds (to avoid wrongly eliminating the higher order terms). For more infos, see https://www.cs.duke.edu/courses/spring11/cps296.3/decoding_rs.pdf
M = [ 0 ] # optional variable to check bounds (so that we do not mistakenly overwrite the higher order terms). This is not necessary, it''s only an additional safe check. For more infos, see the presentation decoding_rs.pdf by Andrew Brown in the doc folder.
# Polynomial constants:
ONE = Polynomial(z0=GF256int(1))
ZERO = Polynomial(z0=GF256int(0))
Z = Polynomial(z1=GF256int(1)) # used to shift polynomials, simply multiply your poly * Z to shift
s2 = ONE + s
# Iteratively compute the polynomials 2s times. The last ones will be
# correct
for l in xrange(0, n-k):
K = l+1
# Goal for each iteration: Compute sigma[K] and omega[K] such that
# (1 + s)*sigma[l] == omega[l] in mod z^(K)
# For this particular loop iteration, we have sigma[l] and omega[l],
# and are computing sigma[K] and omega[K]
# First find Delta, the non-zero coefficient of z^(K) in
# (1 + s) * sigma[l]
# This delta is valid for l (this iteration) only
Delta = ( s2 * sigma[l] ).get_coefficient(l+1) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial).
# Make it a polynomial of degree 0, just for ease of computation with polynomials sigma and omega.
Delta = Polynomial(x0=Delta)
# Can now compute sigma[K] and omega[K] from
# sigma[l], omega[l], B[l], A[l], and Delta
sigma.append( sigma[l] - Delta * Z * B[l] )
omega.append( omega[l] - Delta * Z * A[l] )
# Now compute the next B and A
# There are two ways to do this
# This is based on a messy case analysis on the degrees of the four polynomials sigma, omega, A and B in order to minimize the degrees of A and B. For more infos, see https://www.cs.duke.edu/courses/spring10/cps296.3/decoding_rs_scribe.pdf
# In fact it ensures that the degree of the final polynomials aren''t too large.
if Delta == ZERO or 2*L[l] > K /
or (2*L[l] == K and M[l] == 0):
# Rule A
B.append( Z * B[l] )
A.append( Z * A[l] )
L.append( L[l] )
M.append( M[l] )
elif (Delta != ZERO and 2*L[l] < K) /
or (2*L[l] == K and M[l] != 0):
# Rule B
B.append( sigma[l] // Delta )
A.append( omega[l] // Delta )
L.append( K - L[l] )
M.append( 1 - M[l] )
else:
raise Exception("Code shouldn''t have gotten here")
return sigma[-1], omega[-1]
Polynomial y GF256int son implementaciones genéricas de, respectivamente, polinomios y campos de galois de más de 2 ^ 8. Estas clases se prueban en unidades y son, normalmente, a prueba de errores. Lo mismo ocurre con el resto de los métodos de codificación / descodificación para Reed-Solomon, como Forney y Chien search. El código completo con un caso de prueba rápida para el problema que estoy hablando aquí se puede encontrar aquí: http://codepad.org/l2Qi0y8o
Aquí hay un ejemplo de salida:
Encoded message:
hello world�ꐙ�Ī`>
-------
Erasures decoding:
Erasure locator: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Syndrome: 149x^9 + 113x^8 + 29x^7 + 231x^6 + 210x^5 + 150x^4 + 192x^3 + 11x^2 + 41x
Sigma: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Symbols positions that were corrected: [19, 18, 17, 16, 15]
(''Decoded message: '', ''hello world'', ''/xce/xea/x90/x99/x8d/xc4/xaa`>'')
Correctly decoded: True
-------
Errors+Erasures decoding for the message with only erasures:
Erasure locator: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Syndrome: 149x^9 + 113x^8 + 29x^7 + 231x^6 + 210x^5 + 150x^4 + 192x^3 + 11x^2 + 41x
Sigma: 101x^10 + 139x^9 + 5x^8 + 14x^7 + 180x^6 + 148x^5 + 126x^4 + 135x^3 + 68x^2 + 155x + 1
Symbols positions that were corrected: [187, 141, 90, 19, 18, 17, 16, 15]
(''Decoded message: '', ''/xf4/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00./x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00P/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/x00/xe3/xe6/xffO> world'', ''/xce/xea/x90/x99/x8d/xc4/xaa`>'')
Correctly decoded: False
-------
Errors+Erasures decoding for the message with erasures and one error:
Erasure locator: 77x^4 + 96x^3 + 6x^2 + 206x + 1
Syndrome: 49x^9 + 107x^8 + x^7 + 109x^6 + 236x^5 + 15x^4 + 8x^3 + 133x^2 + 243x
Sigma: 38x^9 + 98x^8 + 239x^7 + 85x^6 + 32x^5 + 168x^4 + 92x^3 + 225x^2 + 22x + 1
Symbols positions that were corrected: [19, 18, 17, 16]
(''Decoded message: '', "/xda/xe1''/xccA world", ''/xce/xea/x90/x99/x8d/xc4/xaa`>'')
Correctly decoded: False
Aquí, la decodificación del borrado siempre es correcta, ya que no utiliza BM en absoluto para calcular el localizador de borrado. Normalmente, los otros dos casos de prueba deben emitir el mismo sigma, pero simplemente no lo hacen.
El hecho de que el problema proviene de BM es evidente aquí cuando se comparan los dos primeros casos de prueba: el síndrome y el localizador de borrado son los mismos, pero la sigma resultante es totalmente diferente (en la segunda prueba se usa BM, mientras que en el primer caso de prueba con borraduras solo BM no se llama).
Muchas gracias por cualquier ayuda o alguna idea sobre cómo puedo solucionar esto. Tenga en cuenta que sus respuestas pueden ser matemáticas o de código, pero explique qué ha ido mal con mi enfoque.
/ EDITAR: todavía no encontré cómo implementar correctamente un decodificador BM errata (ver mi respuesta a continuación). La recompensa se ofrece a cualquiera que pueda solucionar el problema (o al menos guiarme a la solución).
/ EDIT2: tonto, lo siento, solo volví a leer el esquema y encontré que me perdí el cambio en la tarea L = r - L - erasures_count ... He actualizado el código para arreglarlo y volví a aceptar mi respuesta.
He hecho referencia a su código python y lo reescribí por C++ .
Funciona, su información y código de muestra son realmente útiles.
Y encontré que las fallas incorrectas podrían ser causadas por el valor M
Según "Códigos algebraicos para transmisión de datos", el valor
Mno debe ser miembro delif-else.
No obtuve ninguna falla incorrecta después de que se eliminó la M (o simplemente no falla aún)
Muchas gracias por compartir tu conocimiento.
// calculate C
Ref<ModulusPoly> T = C;
// M just for shift x
ArrayRef<int> m_temp(2);
m_temp[0]=1;
m_poly = new ModulusPoly(field_, m_temp);
// C = C - d*B*x
ArrayRef<int> d_temp(1);
d_temp[0] = d;
Ref<ModulusPoly> d_poly (new ModulusPoly(field_, d_temp));
d_poly = d_poly->multiply(m_poly);
d_poly = d_poly->multiply(B);
C = C->subtract(d_poly);
if(2*L<=n+e_size-1 && d!=0)
{
// b = d^-1
ArrayRef<int> b_temp(1);
b_temp[0] = field_.inverse(d);
b_poly = new ModulusPoly(field_, b_temp);
L = n-L-e_size;
// B = B*b = B*d^-1
B = T->multiply(b_poly);
}
else
{
// B = B*x
B = B->multiply(m_poly);
}