原文: Reed–Solomon codes for coders
参考: AN2407.pdf
WIKI: 里德-所罗门码
实现：Pypi ReedSolo
#译注：最近看到了RS码，发现还挺有意思的，找了一些资料学习了下，发现对于程序员来说，从这篇看起会比较容易。看完以后想着翻译一下试试，看看自己到底看懂了多少，于是就有了这篇。本文有部分错误，以及一些排版不对的地方，有兴趣的还是看原文更好:)
为程序员写的Reed-Solomon码解释
Reed-Solomon纠错码（以下简称RS码）广泛用于数据存储（如CD）和传输应用中。然而，在这些应用中，码字是藏在了电子设备里，所以无法一窥它们的模样以及它们是如何生效的。有些复杂的条形码设计也采用了RS码，能够暴露出所有的细节，对于想要获得这种技术如何生效的第一手技术的爱好者，这是一种很有趣的方式。
在这篇文章里，我是试图从程序员的视角（而不是数学家的视角）来介绍RS码的基本原理。我会用以当下流行的QR码作为例子来介绍。我选择了Python（主要是因为写出来的代码看起来整洁美观），但是我也会介绍一些不那么显而易见的Python特性，以便那些不熟悉Python的人也能看懂。里头涉及到的数学知识对读者有一定要求，并且一般是大学才教授的，但是应当能让对高中代数掌握较好的人看懂。
内容：
1 QR码结构
1.1 掩码
1.2 格式信息
1.3 数据
1.4 解码
2 BCH码
2.1 BCH错误检测
2.2 BCH纠错
3 有限域理论
3.1 乘法
3.2 基于对数的乘法
3.3 除法
3.4 多项式
4 RS码
4.1 RS生成多项式
4.2 RS编码
4.3 伴随式(Syndrome)计算
4.4 消除(erasure)纠正
4.5 错误(error)纠正
4.6 消除和错误纠正
1. QR码结构
这一节详细介绍QR码的结构。本节的信息不完整，这是有意为之，只介绍了一个小的21x21的QR码（也被称为version 1）的常见特征，否则需要用好几页表格才能详细定义更大的QR码。完整的标准参见ISO标准18004. 很不幸地，这个文档并不是可以开放获取的。期望了解完整细节的读者，可以去查看众多开源QR码的实现细节。
这是一个用于当例子的QR码。它由深色和浅色的方格组成，在条形码领域被称为“模块”（module）。在角落的3个方形定位器模式是QR码的典型可见特征。
1.1 掩码
编码器会采用一个掩码过程，从而避免可能对扫描器产生混淆的特征，诸如像定位器模式的形状，或者是大片的空白区域。掩码逆转某些模块（白色变成黑色，黑色变成白色），保留其他模块不变。
参考下面的图示，红色的区域编码了格式信息，并使用了一个固定的掩码模式。数据区域（黑白部分）使用一个可自定义的掩码模式。当QR码被创建的时候，编码器会尝试不同的掩码，选择使得不期望出现特征出现最少的那个。被选择的掩码模式信息会被保存在格式信息中，使得解码器知道该用哪个。浅灰色的区域是固定模式，因此并不包含任何信息。此外，在定位器模式中，也会有timing pattern(译注：不懂)，包含可选的黑、白模块。
通过使用异或（XOR，eXclusive-or，通常在变成语言中用 ^ 来表示），掩码过程可以很容易地被加载/移除。对格式信息的反掩码操作如下所示。逆时针读取左上角的定位器模式，我们能够得到下面的比特序列，白色表示0，黑色表示1。
Input 101101101001011
Mask ^ 101010000010010
Output 000111101011001
1.2 格式信息
格式信息有另一份可辨别的副本，因此即使其中一份被毁坏，也仍然有机会被识别。副本被分成两个部分，分别放在另外两个定位器的边上，同样也是逆时针方向阅读（沿着左下角定位器往上，然后是右上角定位器边缘从左往右）。
格式信息的前2位给出了用于数据的纠错级别。这个尺寸的QR码包含26字节信息，其中一些用于保存原数据，一些用于保存校验码，如下表所示。左边第一列只是给纠错级别起了个简单的名字。
纠错级别 级别指示器 纠错码字节数 原数据字节数
L 01 7 19
M 00 10 16
Q 11 13 13
H 10 17 9
格式信息中的接下来三个bit用于指定对数据区域使用的掩码模式。模式使用6x6放歌的方式定义，根据需要不断重复以覆盖整个区域。模式如下所示，包含了对应的数学公式指明（掩码中的）每个模块是否是黑色（i和j分别是行列从0开始的编号，从左上角开始算起）
格式信息中剩下的10个bit是用于对格式信息本身的错误校验，会在后续章节中解释。
1.3 数据
下图展示了经过反掩码操作后的放大了的图样。不同的区域被标记出来了，包括数据区域的边缘。
数据比特从右下角开始，沿着最右边的两列向上走“之”字形。前三个字节分别是 01000000 11010010 01110101（译注：可参见图中右下角的小箭头）。接下来两列从上向下读取，因此接下来的字节是 01000111 。当读取到底部后，再反过来从下往上读取接下来两列……按照这种方式一直读到最左边的列（如果有必要，跳过timing pattern）。下面是完整的十六进制表示的数据：
原始信息: 40 d2 75 47 76 17 32 06 27 26 96 c6 c6 96 70 ec
错误纠正: bc 2a 90 13 6b af ef fd 4b e0
1.4 解码
最后的步骤是将信息解码成可读格式。前4个bit指明了信息是如何编码的。QR码使用几种不同的编码方案，使得不同类型的信息可以更高效地被存储，可参见下表的总结。在方案指示器之后是“长度”字段，告诉解码器保存了多少个字符。字符的长度取决于指定的编码方案。
方案名称 模式指示器 长度字节数 数据字节数
数字 0001 10 10 bits per 3 digits
字母数字 0010 9 11 bits per 2 characters
字节 0100 8 8 bits per character
汉字 1000 8 13 bits per character
(长度字节数只对小的QR码有效）
我们的样例数据开头是 0100 ，表明接下来是每个字符8个bit。接下来8个bit是长度字段，00001101（10进制的13），表明有13个字符。之后才是数据字符。前两个是00100111和01010100（对于ASCII字符 ' 和 T）。有兴趣的读者可以自行解码后续的字符。（译注：用微信二维码扫扫就行了。。。）
在数据bit之后是另外一个4bit模式指示器，可以跟前一个不同，从而允许在一个QR码中混合多个编码方案。如果没有其他数据了，用 0000 来标记结尾（注意，标准允许忽略这个标记，如果存储空间不够的话）。
2 BCH码
格式信息是用BCH码编码的，能够纠正一定数量的bit错误。BCH码是RS码的普遍形式（所有的RS码都是BCH码）。在QR二维码中，BCH码只用于格式信息，比数据信息中用到的RS码要简单得多，因此我们可以先从这里开始。
2.1 BCH错误检测
用于检测编码信息的过程类似整数除法，但是使用异或来代替减法。格式串应该能够被称为“生成子”（generator）的码“整除”。QR码使用 10100110111 这个生成子。下面使用前述格式串 000111101011001 演示了这个过程：
00011
10100110111 ) 000111101011001
^ 10100110111
010100110111
^ 10100110111
00000000000
下面这个Python函数实现了这个过程
def qr_check_format(fmt):
g = 0x537 # = 0b10100110111 in python 2.6+,0b表示二进制
for i in range(4,-1,-1):
if fmt & (1 << (i+10)):
fmt ^= g << i
return fmt
#Python注记1：range函数可能对于非python程序猿来说不够明确。它产生一个数字序列 [4, 3, 2, 1, 0]。源于C的语言中，它类似于 for (i = 4; i >= 0; i--); 在pascal类语言中则类似于 for i := 4 downto 0 。
之所以为5位数组，是 000111101011001位数减10100110111位数
#Python注记2：&操作符是“按位与”，<<操作符是左移位，与C类语言一致。
#译注：这个函数的返回值是参数 fmt 除以生成子的余数
这个函数也可以被用于编码5比特的格式信息：
encoded_format = (format<<10) ^ qr_check_format(format<<10)
#译注：由于format左移了10个bit，因此这里用 ^ 和用 + 是等价的。实际上因为这里的 + 是按位加（其实也就是异或了），所以它也等同于 - ，这一点对于理解它很重要。如果记格式信息为F，生成子为G，(F<<10)/G的商为Q、余数为R (即F<<10 == Q*G + R)，则最终的编码信息 C = (F << 10) ^ ((F << 10) mod G) = (Q*G + R) - R = Q*G，从而C应当是能够被G整除的。如果收到的C不能被G整除，说明传输出错了。
读者也许会对修改此函数让它能除以不同数字感兴趣。例如，更大些的QR码包含6个比特的版本信息和12个错误校验码，并使用生成子 1111100100101 。
使用数学的规范格式，这些二进制数字可以用一个系数为"整数模2"的多项式来描述。数字中的每一个bit是多项式中对应一项的系数（译注：该项的幂等于该bit在数字中的位置）。例如：
10100110111 = 1 x^10 + 0 x^9 + 1 x^8 + 0 x^7 + 0 x^6 + 1 x^5 + 1 x^4 + 0 x^3 + 1 x^2 + 1 x^1 + 1
= x^10 + x^8 + x^5 + x^4 + x^2 + x^1 + 1
2.2 BCH纠错
如果qr_check_format(fmt)得到的余数不是0，那么这个码被损坏或者是读取错误了（译注：即使是0也不能100%保证就对了）。下一步是要找出哪一个格式码最可能是原数据。虽然对于BCH码有许多复杂的解码算法，但是在这里杀鸡用不上牛刀。因为总共只有32个格式串（译注：15个bit中前5个是原信息，后10个是根据原信息生成的），因此遍历找出所有码字中与fmt不同位数最小的那个会更简单（这个被称为汉明距离, hamming distance）
def hamming_weight(x): #不同bit的数量
weight = 0
while x > 0:
weight += x & 1
x >>= 1
return weight
def qr_decode_format(fmt):
best_fmt = -1
best_dist = 15
for test_fmt in range(0,32):
test_code = (test_fmt<<10) ^ qr_check_format(test_fmt<<10)
test_dist = hamming_weight(fmt ^ test_code)
if test_dist < best_dist:
best_dist = test_dist
best_fmt = test_fmt
elif test_dist == best_dist: #如多个码字与fmt距离相同，则都不选
best_fmt = -1
return best_fmt
如果发现fmt不能被唯一地解码，则函数qr_decode_format返回-1。这种情况是由于有多个码字与fmt具有相同的距离。
>>> print(qr_decode_format(int("000111101011001",2))) # no errors
3
>>> print(qr_decode_format(int("111111101011001",2))) # 3 bit-errors
3
>>> print(qr_decode_format(int("111011101011001",2))) # 4 bit-errors
-1
#Python注记： >>> 是python解释器的提示符
3 有限域理论
在讲解用于编码数据的RS码之前，还需要介绍一点数学。类似于乘法和除法，我们定义两个对应的操作，应用于8-bit字节（译注：这里应该指的是整数，下同），且其结果也是8-bit字节。许多类似的算术规则对于新的定义也仍然有效，例如，任意元素乘以1（幺元）结果不必，任意元素乘以0（零元）得0，不允许除以0。其他有用的数学属性（例如分布率）也仍然有效。（基于这些操作的元素集合）被称为一个有限域(finite field)，或者叫伽罗华域(Galois field)。
3.1 乘法
乘法是基于上面定义的多项式。以多项式的形式把输入写下来，然后像我们熟悉的那样，用乘法分配率来计算。我们用 10001001 x 00101010 来举例：
(x^7 + x^3 + 1) * (x^5 + x^3 + x^1)
= x^7 * (x^5 + x^3 + x^) + x^3 * (x^5 + x^3 + x^1) + 1 * (x^5 + x^3 + x^1)
= x^12 + x^10 + x^6 + x^5 + x^4 + x^3 + x^1
#译注：由于最后的加法是异或（模2加），因此系数为偶数的项都消去了。
相同的结果也可以通过竖式计算的一个修改版得到，只要把其中的加改成异或即可：
10001001
* 00101010
10001001
^ 10001001
^ 10001001
1010001111010
下面这个Python函数实现了这个多项式乘法。
def poly_mult(x,y):
z = 0
i = 0
while (y>>i) > 0:
if y & (1<<i):
z ^= x<<i
i += 1
return z
当然，由于结果不能用8-bit来表示了（在这个例子里结果是13个bit），所以我们在得到最终结果之前还需要一个步骤：模100011101，使用前面提到的除法过程。在本例中，最终的结果是 11000011。
#译注：100011101是有限域GF(256)的一个本原多项式（记为p）。简单地理解，GF(256)就是8bit整数0~255，记零元为0，幺元为1，对于α=00000010，α的i次方(1<=i<=254)%p都在F(256)中且各不相同、不是0、也不是1。特别地，α的255次方模p等于1；因而α是GF(256)的一个本原根(primitive root)。α和p并不是唯一的，只是通常都这么用。
3.2 基于对数的乘法
上述过程并不是最适合用来计算伽罗华域的乘法。对两个数进行乘法的过程需要8次循环，除法的过程也需要8次循环（当然，这2种循环可以合并以简化计算），然而我们实际上可以通过使用查表的方式来避免循环。一个科学的结果是在内存中建立完整的乘法表，但是那需要创建64k大小的表格（译注：256×256）。下面给出的结果要紧凑的多。
首先，注意到使用00000010来进行乘法是相当简单的（通常把它记为α），只需要向左移1位，然后与100011101进行异或操作（如果有必要）。下面是α的前几次幂：
α^0 = 00000001 α^4 = 00010000 α^8 = 00011101 α^12 = 11001101
α^1 = 00000010 α^5 = 00100000 α^9 = 00111010 α^13 = 10000111
α^2 = 00000100 α^6 = 01000000 α^10 = 01110100 α^14 = 00010011
α^3 = 00001000 α^7 = 10000000 α^11 = 11101000 α^15 = 00100110
#译注：α^7 * α超过了8-bit，需要与100011101异或得到α^8，依此类推。
如果这个表格依次类推，α的i次幂不会重复，直到α^255=00000001。由此，这个域中的每个元素，除了0之外，都是α的某次幂。α被称为伽罗华域的本原元素（primitive element）或者生成子(generator)。
仔细观察，可以发现另一种实现乘法的方法：把α的幂加起来。
10001001 * 00101010 = α^74 * α^142 = α^74 + 142 = α^216 = 11000011
问题是，我们怎么计算10001001是α的几次幂呢？这个问题被称为离散对数（discrete logarithm），目前没有已知的高效解法。然而，因为域中只有256个元素，我们可以很容易地建立起一个指数表；而同时，对应的对数表也也有用。
gf_exp = [1] * 512 #512个元素的列表. Python 2.6+可以考虑使用bytearray类型
gf_log = [0] * 256
x = 1
for i in range(1,255):
x <<= 1
if x & 0x100:
x ^= 0x11d
gf_exp[i] = x
gf_log[x] = i
for i in range(255,512):
gf_exp[i] = gf_exp[i-255]
表gf_exp的冗余（实际只需要256，但是重复了一遍）是为了简化乘法。这样我们就不用保证gf_log[x]和gf_log[y]会超出表的范围。
def gf_mul(x,y):
if x==0 or y==0:
return 0
return gf_exp[gf_log[x] + gf_log[y]]
3.3 除法
实现指数表的另一个好处是可以使用幂的差来定义除法。在下面的代码里，加上255是为了保证差不是负数。
def gf_div(x,y):
if y==0:
raise ZeroDivisionError()
if x==0:
return 0
return gf_exp[gf_log[x] + 255 - gf_log[y]]
#Python注记：raise语句抛出一个异常从而终止gf_div函数的执行。
使用这种方式来实现除法，对于任意x和非0元素y，gf_div(gf_mul(x, y), y) == x。
3.4 多项式
继续介绍RS码之前，我们需要定义一些伽罗华域多项式（系数属于伽罗华域）的操作。这可能会带来一点混淆，因为伽罗华域的元素本身也是多项式（译注：是系数仅为0、1的多项式）。我的建议是不要想太多。更混乱的是，x也被用来当作占位符（多项式里的未知数）。这个x和前面提到的多项式的x没有任何关系，因此别把它们混起来。
前面给出的伽罗华域元素的二进制表示发在这里用起来会很罗嗦，因此我换成十六进制。
00000001 x^4 + 00001111 x^3 + 00110110 x^2 + 01111000 x^1 + 01000000
= 01 x^4 + 0f x^3 + 36 x^2 + 78 x^1 + 40
在Python里，多项式可以按x的幂递减的list来表示，list中的元素是对于项的系数，从而上述多项式变成了[0x01, 0x0f, 0x36, 0x78, 0x40]。（实现中也可以用反过来的顺序，各有优劣。）
第一个函数将一个多项式和一个标量相乘：
def gf_poly_scale(p,x):
r = [0] * len(p)
for i in range(0, len(p)):
r[i] = gf_mul(p[i], x)
return r
#Python程序员注意：这个函数不是用Pythonic的方式实现的，它可以用列表推导式(list comprehension)更优雅地实现（译注：指的是可以用return [gf_mul(i, x) for i in p] 来实现）。但是我限制使用这些语言特性以便它可以更方便地被转换成其他语言。
这个函数将两个多项式"加"起来（照旧使用异或操作）
def gf_poly_add(p,q):
r = [0] * max(len(p),len(q))
for i in range(0,len(p)):
r[i+len(r)-len(p)] = p[i]
for i in range(0,len(q)):
r[i+len(r)-len(q)] ^= q[i]
return r
下一个函数将两个多项式乘起来：
def gf_poly_mul(p,q):
r = [0] * (len(p)+len(q)-1)
for j in range(0, len(q)):
for i in range(0, len(p)):
r[i+j] ^= gf_mul(p[i], q[j])
return r
最后，我们需要对一个多项式求职，对于某个x值，求出其标量结果。我们引入秦九韶算法(Horner Scheme, 霍纳算法，中国人很牛逼有没有!)来避免计算x的n次幂：
01 x^4 + 0f x^3 + 36 x^2 + 78 x^1 + 40
= (((01 x^ + 0f) x^ + 36) x^ + 78) x^ + 40
比如 用秦九韶算法求多项式p(x)=3x^-2x^+x+7在x=3处的值
解答（注意：这里没有在有限域里面进行，如果是有限域里面，乘和加法改成相应的函数即可）
p(x)=3x^+0x^-2x^+0x^+x+ =(3x^+0x^-2x^+0x+)x+ =((3x^+0x^-2x+)x+)x+ =(((3x^+0x-)x+)x+)x+ =((((3x+)x-)x+)x+)x+ v[]=3x+= v[]=v[]x-= v[]=v[]x+= v[]=v[]x+= v[]=v[]x+=
def gf_poly_eval(p,x):
y = p[0]
for i in range(1, len(p)):
y = gf_mul(y,x) ^ p[i]
return y
4 RS码
好了，背景知识都介绍完了，我们可以开始看看RS码了。 #译注吐槽：终于进入重点了。
4.1 RS生成多项式
RS码使用类似BCH码的生成多项式，将 (x - α^i) 乘起来。例如：
g(x) = (x - α^0) (x - α^1) (x - α^2) (x - α^3)
= 01 x^4 + 0f x^3 + 36 x^2 + 78 x^1 + 40
下面这个函数计算了对于指定n个校验符号的RS码需要的生成多项式：
def rs_generator_poly(nsym):
g = [1]
for i in range(0,nsym):
g = gf_poly_mul(g, [1, gf_exp[i]])
return g
这个函数不是很高效，因为它在每次循环中陆续分配了更大的数组给g；但是在实际应用中它并不是瓶颈，优化癖读者有兴趣的话可以重写它，让它只给g分配一次。
4.2 编码
类似BCH码，RS码也是用类似于整数除法的流程来编码。下面的例子展示了数据 12 34 56 （注记：这是16进制表示的）是如何被编码的：
12 da df
01 0f 36 78 40 ) 12 34 56 00 00 00 00
^ 12 ee 2b 23 f4
da 7d 23 f4 00
^ da a2 85 79 84
df a6 8d 84 00
^ df 91 6b fc d9
37 e6 78 d9
余数和原数据连起来得到编码后的数据 12 34 56 37 e6 78 d9。下面的这个函数实现了编码过程（译注：可以参考2.1的BCH编码过程 (fmt<<10) ^ qr_check_format(fmt << 10)，那里有个译注解释为什么这么编码）：
def rs_encode_msg(msg_in, nsym):
gen = rs_generator_poly(nsym)
msg_out = [0] * (len(msg_in) + nsym)
for i in range(0, len(msg_in)):
msg_out[i] = msg_in[i]
for i in range(0, len(msg_in)):
coef = msg_out[i]
if coef != 0:
for j in range(0, len(gen)):
msg_out[i+j] ^= gf_mul(gen[j], coef)
for i in range(0, len(msg_in)):
msg_out[i] = msg_in[i]
return msg_out
下面这个例子演示了编码函数如何对第一节的样例QR码中的数据进行编码。它计算出的结果中的第二行（译注：输出的第二行）与样例QR码中解码出来的纠错码相符
>>> msg_in = [ 0x40, 0xd2, 0x75, 0x47, 0x76, 0x17, 0x32, 0x06,
... 0x27, 0x26, 0x96, 0xc6, 0xc6, 0x96, 0x70, 0xec ]
>>> msg = rs_encode_msg(msg_in, 10)
>>> for i in range(0,len(msg)):
... print(hex(msg[i]), end=' ')
...
0x40 0xd2 0x75 0x47 0x76 0x17 0x32 0x6 0x27 0x26 0x96 0xc6 0xc6 0x96 0x70 0xec
0xbc 0x2a 0x90 0x13 0x6b 0xaf 0xef 0xfd 0x4b 0xe0
#Python版本注记：这是Python 3.0+的例子，因为print的语法变了。早版本的Python可以用这句来代替：“print hex(msg[i]), ” （包含末尾的逗号）
4.3 伴随式(Syndrome)计算
RS码的译码操作需要多个步骤。第一个步骤就是计算数据的伴随式。把数据当成一个多项式，使用x = α^0, α^1, α^2, ..., α^n对其求值（译注：得到n个伴随式的值）。因为这些x值使得生成多项式的值为0，因此如果读取到的数据没有损坏，结果应当是0。如果不是，这些伴随式里就包含了完成纠错所必需的信息。计算伴随式的实现很简单：
def rs_calc_syndromes(msg, nsym):
synd = [0] * nsym
for i in range(0, nsym):
synd[i] = gf_poly_eval(msg, gf_exp[i])
return synd
继续上面的例子，我们可以看到（编码后的）数据的伴随式的确都是0；而引入一个错误以后就会得到非零伴随式。
>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
>>> msg[0] = 0 # deliberately damage the message
>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)
[64, 192, 93, 231, 52, 92, 228, 49, 83, 245]
# === 译注开始 ===
收到的数据 r(x) 是由编码后的数据（RS码的码字） c(x) 和错误 e(x) 组合得到的：
r(x) = c(x) + e(x)
记v=nsym由于c(x)在α^0, α^1, α^2, ..., α^(v-1)都是生成多项式的根（去看生成多项式的定义就知道了，注意+实际上也是-），而c(x)能够表示为 q(x) * g(x) （不明白的话可以参看2.1节的译注），所以这里对c(x)求值的结果必然都是0。
如果对r(x)的求值都是0，说明e(x)也都是0，说明要么是没问题，要么就是错到根本发现不出来（正好是另一个码字）。如果r(x)不等于0，那么伴随式S正好就是所有e(x)的值：
S[0] = e(1), S[1] = e(α^1), S[2] = e(α^2), ..., S[v-1] = e(α^(v-1))
如果错误的码字个数不超过v/2，通过S是可以还原信息的。
# === 译注结束 ===
4.4 消除(erasure)纠正
如果错误的位置是已知的（译注：某些信道可以预知，比如BEC信道，wiki这么说的，我也不知道是啥；不过QR码是另一个例子），纠正它是最简单的。这被称为消除纠正。对于每一个添加的纠错码，都可以纠正一个消除（译注：这里应该是说，添加的n个纠错码能够保证纠正n个消除码字）。如果错误位置未知，那么对于一个错误，需要2个错误校验符号。这在实际应用中非常有用，比如QR码的某些位置被覆盖或被剪掉什么的。扫描器很难知道发生了什么，所以不是所有的QR码扫描器都能纠正消除。
Forney算法被用来纠正消除，实现如下：
def rs_correct_errata(msg, synd, pos):
#计算错误定位器多项式 calculate error locator polynomial
q = [1]
for i in range(0, len(pos)):
x = gf_exp[len(msg)-1-pos[i]]
q = gf_poly_mul(q, [x,1])
#计算错误求值多项式 calculate error evaluator polynomial
p = synd[0:len(pos)]
p.reverse()
p = gf_poly_mul(p, q)
p = p[len(p)-len(pos):len(p)]
#去掉偶数项 formal derivative of error locator eliminates even terms
q = q[len(q)&1:len(q):2]
#纠错 compute corrections
for i in range(0, len(pos)):
x = gf_exp[pos[i]+256-len(msg)]
y = gf_poly_eval(p, x)
z = gf_poly_eval(q, gf_mul(x,x))
msg[pos[i]] ^= gf_div(y, gf_mul(x,z))
#Python注记：这个函数使用数组分片来获取数组的某一部分。表达式 synd[0:len(pos)] 返回synd的前几个元素，而 p[len(p)-len(pos):len(p)] 返回p的最后几个元素。更复杂一点的这个表达式q[len(q)&1:len(q):2] 每隔两个元素提取一个，跳过第一个元素，如果q的长度是奇数（译注：就是提取p中的奇数项了...）
#译注：这个算法我没去看，有兴趣的读者可以看本文提到的另一篇讲解AN2407.pdf。
4.5 错误(error)纠正
更可能的情况是未知位置的错误，所以第一步是找出它们。Berlekamp-Massey算法（通常简称为BM算法）用来计算错误定位多项式。然后我们只需要计算这个（错误定位）多项式的零点（译注：应该就是指多项式的根），这标志了错误的位置。
def rs_find_errors(synd, nmess):
# find error locator polynomial with Berlekamp-Massey algorithm
err_poly = [1]
old_poly = [1]
for i in range(0,len(synd)):
old_poly.append(0)
delta = synd[i]
for j in range(1,len(err_poly)):
delta ^= gf_mul(err_poly[len(err_poly)-1-j], synd[i-j])
if delta != 0:
if len(old_poly) > len(err_poly):
new_poly = gf_poly_scale(old_poly, delta)
old_poly = gf_poly_scale(err_poly, gf_div(1,delta))
err_poly = new_poly
err_poly = gf_poly_add(err_poly, gf_poly_scale(old_poly, delta))
errs = len(err_poly)-1
if errs*2 > len(synd):
return None # too many errors to correct
# find zeros of error polynomial
err_pos = []
for i in range(0, nmess):
if gf_poly_eval(err_poly, gf_exp[255-i]) == 0:
err_pos.append(nmess-1-i)
if len(err_pos) != errs:
return None # couldn't find error locations
return err_pos
#译注：这个算法我也没细看，大体上，它用到了牛顿恒等式（根据4.3计算的伴随式）来生成错误定位多项式，然后求值得到错误的位置。更详细一点的信息可参考AN2407.pdf
下面是一个纠正了编码后数据中3个错误的例子：
>>> print(hex(msg[10]))
0x96
>>> msg[0] = 6
>>> msg[10] = 7
>>> msg[20] = 8
>>> synd = rs_calc_syndromes(msg, 10)
>>> pos = rs_find_errors(synd, len(msg))
>>> print(pos)
[20, 10, 0]
>>> rs_correct_errata(msg, synd, pos)
>>> print(hex(msg[10]))
0x96
4.6 消除和错误纠正
为了能够纠正错误和消除，我们需要让已知的消除不影响查找错误位置。这可以通过计算Forney syndrome来实现，如下所示：
def rs_forney_syndromes(synd, pos, nmess):
fsynd = list(synd) # make a copy
for i in range(0, len(pos)):
x = gf_exp[nmess-1-pos[i]]
for i in range(0,len(fsynd)-1):
fsynd[i] = gf_mul(fsynd[i], x) ^ fsynd[i+1]
fsynd.pop()
return fsynd
Forney syndrome可以用来替换常规错误位置查找中syndrome。
下面的这个rs_correct_errata函数给出了完整的过程，在msg_in被消除的位置中使用-1来表示。
def rs_correct_msg(msg_in, nsym):
msg_out = list(msg_in) # copy of message
# find erasures
erase_pos = []
for i in range(0, len(msg_out)):
if msg_out[i] < 0:
msg_out[i] = 0
erase_pos.append(i)
if len(erase_pos) > nsym:
return None # too many erasures to correct
synd = rs_calc_syndromes(msg_out, nsym)
if max(synd) == 0:
return msg_out # no errors
fsynd = rs_forney_syndromes(synd, erase_pos, len(msg_out))
err_pos = rs_find_errors(fsynd, len(msg_out))
if err_pos == None:
return None # error location failed
rs_correct_errata(msg_out, synd, erase_pos + err_pos)
synd = rs_calc_syndromes(msg_out, nsym)
if max(synd) > 0:
return None # message is still not right
return msg_out
#Python注记：erase_pos和err_pos这两个数组用 + 运算连接起来。
这是实现一个全功能的错误/消除纠正RS解码器所需的最后一个片段。
#译注："Bobmath"将上述代码整合以后提交到了pypi，包名叫reedsolo，可以在 https://pypi.python.org/pypi/reedsolo 找到。稍微要注意的一点是，上述的所有代码实际上是针对GF(256)上面的RS码实现的，且使用了固定的生成多项式，所以并不是完全通用。
--
转载请注明出自 https://www.felix021.com/blog/read.php?2116 ，如是转载文则注明原出处，谢谢:)
RSS订阅地址: http://www.felix021.com/blog/feed.php 。