标题：[译] 为程序员写的Reed-Solomon码解释 出处：Felix021 时间：Sun, 24 Mar 2013 13:48:03 +0000 作者：felix021 地址：https://www.felix021.com/blog/read.php?2116 内容： 原文: 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码的典型可见特征。 点击在新窗口中浏览此图片 https://www.felix021.com/blog/attachment.php?fid=466 原图 1.1 掩码 编码器会采用一个掩码过程，从而避免可能对扫描器产生混淆的特征，诸如像定位器模式的形状，或者是大片的空白区域。掩码逆转某些模块（白色变成黑色，黑色变成白色），保留其他模块不变。 参考下面的图示，红色的区域编码了格式信息，并使用了一个固定的掩码模式。数据区域（黑白部分）使用一个可自定义的掩码模式。当QR码被创建的时候，编码器会尝试不同的掩码，选择使得不期望出现特征出现最少的那个。被选择的掩码模式信息会被保存在格式信息中，使得解码器知道该用哪个。浅灰色的区域是固定模式，因此并不包含任何信息。此外，在定位器模式中，也会有timing pattern(译注：不懂)，包含可选的黑、白模块。 点击在新窗口中浏览此图片 https://www.felix021.com/blog/attachment.php?fid=467 原图 通过使用异或（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开始的编号，从左上角开始算起） 点击在新窗口中浏览此图片 https://www.felix021.com/blog/attachment.php?fid=468 原图 格式信息中剩下的10个bit是用于对格式信息本身的错误校验，会在后续章节中解释。 1.3 数据 下图展示了经过反掩码操作后的放大了的图样。不同的区域被标记出来了，包括数据区域的边缘。 点击在新窗口中浏览此图片 https://www.felix021.com/blog/attachment.php?fid=469 原图，这是个SVG格式的图像。 数据比特从右下角开始，沿着最右边的两列向上走“之”字形。前三个字节分别是 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+ 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 。 #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<>> 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)计算 （译注：这个要怎么翻译啊。。Syndrome是症状的意思，这里的确也是在计算收到的RS码码字的症状，从而判断是否接受到了错误的码字。） （译注：后来看了各种文档以后，得知一般都翻译为“伴随式”） 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码实现的，且使用了固定的生成多项式，所以并不是完全通用。 Generated by Bo-blog 2.1.0