码迷,mamicode.com
首页 > 其他好文 > 详细

写给编码人员的里德-所罗门编码

时间:2018-06-08 15:58:12      阅读:761      评论:0      收藏:0      [点我收藏+]

标签:class   include   基本   Python程序员   pre   交换   its   sage   就会   

写给编码人员的里德-所罗门编码

先用谷歌翻译Reed–Solomon codes for coders,然后自己修正了下。by Dan Giant

 

  纠错码是纠正错误的信号处理技术。它们现在无处不在,例如通信(手机,互联网),数据存储和档案(硬盘,光盘CD / DVD /蓝光,档案磁带),仓库管理(条形码)和广告(QR码)等。

  里德 - 所罗门纠错码是一种特定类型的纠错码。它虽然是最古老的纠错码之一,但是仍然被广泛使用,因为它的定义非常清晰,并且现在可以在公共领域获得几种高效的算法。

  通常,纠错码是隐藏的,大多数用户甚至不知道它们,也不知道它们何时被使用。然而,它们是一些应用程序的关键组件,例如通信或数据存储。事实上,每隔几天随机丢失一次数据的硬盘将毫无用处,而且只有在无云天气的情况下才能打电话的手机很少会被使用。使用纠错码可以将损坏的消息恢复为完整的原始消息。

  条码和QR码是研究的有趣应用,因为它们具有在视觉上显示纠错码的特殊性,使好奇的用户容易访问这些码。

  在本文中,我们将尝试从程序员而不是数学家的角度介绍里德 - 所罗门码的原理,这意味着我们将更多地关注实践而不是理论,尽管我们也将解释理论,但只有直觉和实施的必要知识。将提供该领域内值得注意的参考资料,以便有兴趣的读者可以随意深入研究数学理论。我们将提供来自流行的QR码条码系统以及工作代码样本的真实示例。我们选择使用Python作为示例(主要是因为它看起来非常类似于伪代码),但我们将尝试向那些不熟悉它的人解释任何非显而易见的功能。所涉及的数学是高级的,因为它通常不是在大学水平以下教授的,但对于掌握高中代??数的人应该是可以理解的。

  我们首先将轻松介绍纠错码原理背后的直觉,然后在第二部分中我们将介绍QR码的结构设计,换句话说,信息如何存储在QR码中以及如何读取和生成QR码,以及第三部分,我们将通过实施Reed-Solomon解码器来研究纠错码,并快速介绍更大的BCH码族,以便可靠地读取损坏的QR码。

  请注意,好奇的读者可以在附录讨论页面找到扩展信息。 

 

纠错原理

  在详细说明代码之前,了解纠错背后的直觉可能很有用。事实上,尽管纠错码在数学上似乎令人生畏,但大多数数学运算都是高中级的(除Galois Fields外,其实对于任何程序员来说都很容易和通用:它只是对整数进行模运算数)。然而,纠错码背后数学独创性的复杂性掩盖了相当直观的目标和机制。

    纠错码可能看起来像一个挺难的数学概念,但实际上它们是基于一个巧妙的数学实现的直观理念:让我们使数据结构化,以便我们可以在数据损坏时,通过“修复”结构,“猜测”损坏的数据是什么。在数学上,我们使用伽罗瓦域的多项式来实现这个结构。

  让我们举一个更实用的比喻:假设您想要将消息传递给其他人,但这些消息可能会一路被破坏。纠错码的主要观点是,我们可以使用一小组仔细选择的单词,即“缩小字典”,而不是使用整个字典词典,以便每个字词与其他字词不同。这样,当我们收到一条消息时,我们只需要在缩小的字典内查找1)检测哪些字被损坏(因为它们不在我们缩小的字典中); 2)通过查找我们词典中最相似的单词来纠正被破坏的单词。

  我们举一个简单的例子:我们有一个缩小的字典,只有三个4个字母的单词:this,that和corn。假设我们收到一个损坏的单词:co **,其中*是数据丢失。由于我们的字典中只有3个单词,因此我们可以轻松地将收到的单词与我们的字典进行比较,以找到最接近的单词。在这种情况下,它是corn。因此,丢失的字母是rn。

  现在让我们架设收到的单词是“th**”。现在的问题是我们的字典中有两个词与收到的词匹配:this和that。在这种情况下,我们不能确定它是哪一个,因此我们无法解码。这意味着我们的字典不是很好,我们应该用另一个有更多不同的单词来替换,例如dash,以最大化每个单词之间的差异。这种差异,或者更确切地说,是我们字典中任意两个单词之间的最小不同字母数,被称为我们字典的最大汉明距离。确保字典中的任何两个单词在同一位置只共享最少数量的字母称为最大可分离性

  大多数纠错码都使用同样的原则:我们只生成一个简化字典,其中只包含具有最大可分离性的单词(我们将在第三部分详细介绍如何做),然后我们只用这个简化字典的单词进行通信。 伽罗瓦域提供的是结构(即,简化字典的基础),而里德 - 所罗门是一种自动创建合适结构的方法(为数据集量身定制最大可分离性的简化字典),并提供自动方法检测和纠正错误(即在简化字典中查找)。更确切地说,伽罗瓦域是结构(由于它们的循环本质,模以整数),里德 - 所罗门是基于伽罗瓦域的编解码器(编码器/解码器)。

  如果一个单词在通信中被破坏了,这没什么大不了的,因为我们可以通过查看字典并找到最接近的单词来很容易地修复它,这可能是正确的(但是如果输入消息太严重损坏,则有可能选择错误,但概率非常小)。而且,我们的单词越长,它们就越可分离,因为更多的字符可以被破坏而不受任何影响。

  生成最大可分词的字典的最简单方法是让词比实际长。

  我们再来看看我们的例子:

   t h i s
   t h a t
   c o r n

 

  添加一组独特的字符,以便在任何新增位置都没有重复的字符,并添加一个单词来帮助解释:

   t h i s a b c d
   t h a t b c d e
   c o r n c d e f
   t h i n d e f g

 

  请注意,此字典中的每个单词与其他单词至少有5个字符不同,因此距离为5。这允许在已知位置最多有4个错误,这些错误称为丢失,或者在未知位置中的2个错误需要纠正。

  假设发生4次错误:

   t * * * a b * d

 

  然后可以在字典中搜索4个未丢失的字符,并找到与这4个字符匹配的唯一条目,因为距离为5。

  假设在这些模式之一中有两个错误发生:

  t h i x a x c d
  x h i x a b c d
  t h x x a b c d

 

  这里的问题是错误的位置是未知的。8个中取6个字符的可能子集有28个,因此使用有这6个字符的28个子集中的每个子集进行搜索,因为距离为5,将只有一个匹配6个字符的子集(假设发生2个或更少的错误)。

  通过这些示例,你可以看到数据冗余在恢复丢失信息中的优势:冗余字符可帮助你恢复原始数据。前面的例子显示了粗略的错误纠正方案是如何工作的。Reed-Solomon的核心思想是相似的,将冗余数据附加到基于伽罗华域数学的消息上。原始纠错解码器与上述错误示例类似,搜索与有效消息相对应的接收消息的子集,并选择匹配最多的消息作为纠正的消息。对于较大的消息来说这是不实际的,因此开发了数学算法来在合理的时间内执行纠错。

 

QR码结构 

  本节介绍QR码的结构,QR码是如何将数据存储在QR码中的。本节中的信息故意不完整。这里仅介绍小型21×21尺寸符号(也称为版本1)的最常用功能,但请参阅附录以获取更多信息。

  下图是一个用作示例QR码图。 它由黑色和白色的正方形组成,在条形码中被称为模组。 角落上的三个方形定位器图案是QR码的的图像特征。

技术分享图片

技术分享图片

掩码

   掩码处理用于避免码图中出现可能会使扫码器混淆的特征,例如看起来像定位器的图案,和大面积空白区域等有误导性的图形。 掩码反转某些模块(白色变成黑色,黑色变成白色),同时不影响其他模块。

  在下图中,红色区域编码格式信息并使用固定的掩码图形。 数据区域(黑白)被可变模板掩码处理。 创建QR码时,编码器会尝试多种不同的掩码,并选择使结果中不希望出现的特征最小化的掩码图形。 然后在格式信息中指明所选掩码图形,以便解码器知道使用哪一个。浅灰色区域是不编码任何信息的固定图形。除了明显的定位器图形外,还有定时图形,其中包含交替的明暗模块。

技术分享图片

  掩码变换很容易通过异或操作(在许多编程语言中用插入符号 ^ 表示)应用(或删除)。格式信息的掩码去除如下所示。逆时针读取左上角的定位器图案,我们得到如下的位序列。白色模块表示0,黑色模块表示1。

 

Input       101101101001011
Mask      ^ 101010000010010
Output      000111101011001

 

 

格式化信息

  格式信息有两个相同的副本,因此即使该码元已损坏也仍可被解码。 第二个副本分成两部分,放在另外两个定位器的周围,同时也以逆时针方向读取(在左下角朝上,然后从左上角朝右上角)。

  格式化信息的前两位给出了用于消息数据的纠错级别。 这个尺寸的QR码包含26个字节的信息。 其中一些用于存储消息,一些用于纠错,如下表所示。左栏是级别的名字。 

Error Correction LevelLevel IndicatorError Correction BytesMessage Data Bytes
L 01 7 19
M 00 10 16
Q 11 13 13
H 10 17 9

  接下来的3个位的格式信息选择要在数据区域中使用的掩膜图案。这些图案如下图所示,其中包括数学公式指示模块是否为黑色(i和j分别为行号和列号,并以左上角的0开头)。

 技术分享图片

  格式信息的其余10个位用于纠正格式本身的错误。这将在后面的章节中说明。

 

消息数据

  这是一个更大的图案,显示“未掩码”的QR码。码图的不同区域被指示出来,包括消息数据字节的边界。

技术分享图片

  从右下角开始读取数据位,并以Z字形向上移动两个右列。前三个字节是01000000  11010010  01110101。接下来的两列是向下读取的,所以下一个字节是01000111。一旦到达底部,之后的两列将被向上读取。 继续以上下方式一直到符号的左侧(必要时跳过时序图案)。这是以16进制表示的完整信息。

  消息数据字节: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

 

解码

最后一步是将消息字节解码为可读的内容。 前4个位指示消息如何编码。QR码使用几种不同的编码方案,因此可以高效地存储不同类型的消息。这些总结在下表中。在模式指示符之后是长度字段,它指示存储了多少个字符。长度字段的大小取决于特定的编码。

 

Mode NameMode IndicatorLength BitsData Bits
Numeric 0001 10 10 bits per 3 digits
Alphanumeric 0010 9 11 bits per 2 characters
Byte 0100 8 8 bits per character
Kanji 1000 8 13 bits per character

   (上述长度字段仅适用于较小的QR码。)

  我们的示例消息从0100开始,表示每个字符有8个位。接下来的8个位是长度字段,00001101,十进制13。长度之后是消息的实际字符。前两个是00100111和01010100(撇号和T的ASCII码)。有兴趣的读者可能想要为自己解码其余的消息。

  数据位的最后是一个4个位的模式指示符。它可以与第一个不同,以允许在同一个QR码内混合不同的编码。当没有更多数据要存储时,给出特殊的消息结束码0000。(请注意,如果数据字节不足,标准允许省略消息结束码。)

  此时,我们知道如何解码或读取整个QR码。然而,在现实条件下,QR码很少是完整的:通常,它可以通过手机的相机扫描,在光线条件较差的情况下,或者在部分刮擦导致QR代码划痕的表面上,或在沾污的表面上,等等

  为了使我们的QR码解码器“可靠”,我们必需能够“纠正”错误。本文的下一部分将介绍如何通过构建BCH解码器,更具体地说是Reed-Solomon解码器,来纠正错误。

 

BCH码

  在本节中,我们将介绍一类通用的纠错码:BCH码,它是现代里德-所罗门码的父族,及其基本检测和纠正机制。

  格式化信息用BCH码编码,它能检测和纠正一定数量的比特错误。BCH码是里德-所罗门码的一般化(现代里德-所罗门码是BCH码)。在QR码的情况下,用于格式信息的BCH码比用于消息数据的里德-所罗门码简单得多,因此从用于格式信息的BCH码开始是有道理的。

 BCH错误检测

  检查编码信息的过程与长除法相似,但使用的是亦或而不是减法。 格式码在被代码的所谓生成器“除”时余数应为零。QR格式码使用生成器10100110111。下面的示例代码(000111101011001)中的格式信息说明了此过程。

 

                        00011
10100110111 ) 000111101011001
               ^ 10100110111 
                 010100110111
                ^ 10100110111
                  00000000000

  这是一个实现这个计算的Python函数。

 

1 def qr_check_format(fmt):
2    g = 0x537 # = 0b10100110111 in python 2.6+
3    for i in range(4,-1,-1):
4       if fmt & (1 << (i+10)):
5          fmt ^= g << i
6    return fmt

  Python说明:非Python程序员可能对范围函数不清楚。它产生一个从4到0的数字列表。在C派生语言中,for循环可以写为(i = 4; i> = 0; i--); 在Pascal派生语言中,for i := 4 downto 0。

  Python注2:&运算符按位执行,while<<是左移位。这与C类语言一致。

  该函数也可用于编码5个位的格式信息。

 

encoded_format = (format<<10) ^ qr_check_format(format<<10)

  读者可能会发现这是一个有趣的练习,将这个函数概括为不同的数字。例如,较大的QR码使用发生器1111100100101,包含6个位的版本信息和12个位的纠错位。

  在数学形式表示,这些二进制数被描述为系数是整数模2的多项式。数的每一位是一项的系数。例如:

10100110111 = 1 x10 + 0 x9 + 1 x8 + 0 x7 + 0 x6 + 1 x5 + 1 x4 + 0 x3 + 1 x2 + 1 x + 1 = x10 + x8 + x5 + x4 + x2 + x + 1

  如果由qr_check_format产生的余数不为零,则QR码被损坏或误读。下一步是确定哪种格式码最有可能是预期的格式码(即,在我们的简化字典中查找)。

 

BCH纠错

  尽管存在用于解码BCH码的复杂算法,但在这种情况下它们可能是大材小用。 由于只有32种可能的格式码,因此简单地尝试每种格式码,并选择与有问题的码不同的比特数(称为汉明距离)最小的更加容易。这种查找最接近码的方法称为穷举搜索,并且可能只是因为我们的码字非常少(码字是个有效消息,并且此处只有32个,除此外其他所有的二进制数字都不正确)。

  (请注意,里德-所罗门码也是基于这个原理,但由于可能的码字数量太大,我们无法进行穷举搜索,这就是为什么要设计巧妙但复杂的算法的原因,例如Berlekamp-Massey算法。)

 1 def hamming_weight(x):
 2    weight = 0
 3    while x > 0:
 4       weight += x & 1
 5       x >>= 1
 6    return weight
 7 
 8 def qr_decode_format(fmt):
 9    best_fmt = -1
10    best_dist = 15
11    for test_fmt in range(0,32):
12       test_code = (test_fmt<<10) ^ qr_check_format(test_fmt<<10)
13       test_dist = hamming_weight(fmt ^ test_code)
14       if test_dist < best_dist:
15          best_dist = test_dist
16          best_fmt = test_fmt
17       elif test_dist == best_dist:
18          best_fmt = -1
19    return best_fmt

 

  如果格式代码不能被明确解码,函数qr_decode_format返回-1。当两个或多个格式码与输入具有相同的距离时会发生这种情况。

  要在Python中运行此代码,请首先启动IDLE,Python的集成开发环境。您应该看到版本信息和交互式输入提示>>>。打开一个新窗口,将函数qr_check_format,hamming_weight和qr_decode_format复制到里面,并另存为qr.py。返回到提示符,并键入下面的>>>后面的行。

 

>>> from qr import *
>>> qr_decode_format(int("000111101011001",2))  # no errors
3
>>> qr_decode_format(int("111111101011001",2))  # 3 bit-errors
3
>>> qr_decode_format(int("111011101011001",2))  # 4 bit-errors
-1

  你也可以通过在命令提示符处键入python来启动Python。

  在接下来的部分中,我们将研究有限域算法和Reed-Solomon码,它是BCH码的一个子类。基本思想(即,使用具有最大可分离性的有限字词字典)是相同的,但是由于我们将编码更长的字(256字节而不是2字节),因此要用到更多的码元(在所有8个位上编码,因此有256个不同的可能值 ),我们不能使用这种天真而详尽的方法,因为它会花费太多时间:我们需要使用更聪明的算法,有限域数学将帮助我们做到这一点,通过给我们一个结构。

 

有限域算术

数学域的介绍

  在讨论用于消息的Reed-Solomon码之前,先介绍更多些数学是有用的。

  我们希望为8位字节定义加法,减法,乘法和除法,并总是产生8位字节,以避免任何溢出。天真地说,我们可能会尝试使用这些操作的正常定义,然后用256取模以防止溢出。这正是我们将要做的,也就是所谓的伽罗瓦域2 ^ 8。你可以很容易地想象为什么它适用于所有事情,除了除法:7/5是多少?

  下面简单介绍伽罗瓦域:有限域是一组数字,一个域需要有六个属性:封闭性,结合律,交换律,分配律,单位元和可逆。更简单地说,使用域允许研究该域的数字之间的关系,并将结果应用于遵循相同属性的任何其他域。例如,实数集?是一个域。换句话说,数学的域研究一个数集的构造。

  然而,整数?不是一个域,因为正如我们上面所说的,并非所有的除法都被定义(例如7/5),这违反了乘法逆性质(x如7 * x = 5不存在)。解决这个问题的一个简单方法是使用素数来取模,比如2:这样,我们保证存在7 * x = 5,因为我们只是环绕一下。?对2取模被称为伽罗瓦域,任何可被2整除的数都是伽罗瓦域(因为我们需要使用素数取模),所以256,8位码元的值,可以减少到2 ^ 8 ,因此我们说我们使用2 ^ 8的伽罗瓦域,或GF(2 ^ 8)。 在口语中,2是该域的特征,8是指数,256是该域的基数。关于有限域的更多信息可以在这里找到

  在这里我们将定义通常用于整数运算的数学运算,但适用于GF(2 ^ 8),它基本上是按照常规运算进行的,但是对2 ^ 8取模。

  考虑GF(2)和GF(2 ^ 8)之间关系的另一种方法,是将GF(2 ^ 8)看做表示8个二进制系数的多项式。例如,在GF(2 ^ 8)中,170等于10101010 = 1 * x ^ 7 + 0 * x ^ 6 + 1 * x ^ 5 + 0 * x ^ 4 + 1 * x ^ 3 + 0 * x ^ 2 + 1 * x + 0 = x ^ 7 + x ^ 5 + x ^ 3 + x。两种表示方式都是等价的,只是在第一种情况下,170是十进制的表示,而另一种情况是二进制的表示形式,可以将其视为按照惯例表示多项式(仅在GF(2 ^ p)中使用,如解释这里)。后者通常是学术书籍和硬件实现中使用的表示形式(因为逻辑门和寄存器在二进制级别工作)。对于软件实现来说,优先使用十进制表示,因为这种方式更清晰,更贴近数学的表示(这是我们将在本教程中使用的表示法,除了一些示例将使用二进制表示之外)。

  在任何时候,尽量不要将表示单个GF(2 ^ p)码元的多项式(每个系数是一个位/布尔值:0或1),与表示一组GF(2 ^ p)码元列表的多项式(在这种情况下,多项式等价于消息 + RS码,每个系数是介于0和2 ^ p之间的值,并表示消息 + RS码的一个字符)混淆。我们将首先描述单个码元的操作,然后描述码元列表上的多项式操作。

 

加法和减法

  在基于2的伽罗瓦域中,加法和减法都被替换为异或。这是合乎逻辑的:加模2与XOR完全相同,减模2与加模2完全相同。这是可能的,因为加法和减法在这个伽罗瓦域是无进位的。

  将我们的8位值看作系数为mod 2的多项式:

   0101 + 0110 = 0101 - 0110 = 0101 XOR 0110 = 0011

  以相同的方式(以两个GF(2 ^ 8)整数的二进制表示):

(x2 + 1) + (x2 + x) = 2 x2 + x + 1 = 0 x2 + x + 1 = x + 1

  由于(a ^ a)= 0,每个数都是自身取反,所以(x - y)与(x + y)相同。

  请注意,在书中,你将找到加法和减法来定义对GF域整数的一些数学运算,但实际上,只需XOR(只要你是在伽罗瓦域中,在其他域中不是这样)。

  这是等效的Python代码:

1 def gf_add(x, y):
2     return x ^ y
3 
4 def gf_sub(x, y):
5     return x ^ y # in binary galois field, subtraction is just the same as addition (since we mod 2)

 

乘法

  乘法同样基于多项式乘法。 简单地将输入写成多项式,然后用分配法将它们相乘。举例来说,10001001倍的00101010计算如下。

(x7 + x3 + 1) (x5 + x3 + x) = x7 (x5 + x3 + x) + x3 (x5 + x3 + x) + 1 (x5 + x3 + x)
x12 + x10 + 2 x8 + x6 + x5 + x4 + x3 + x
x12 + x10 + x6 + x5 + x4 + x3 + x

  同样的结果可以通过修改版本的标准中学乘法过程来获得,其中我们用异或代替加法。

       10001001
*      00101010
      10001001
^   10001001
^ 10001001
  1010001111010

  注意:这里的XOR乘法是没有进位的!如果你进行了进位,你会得到错误的结果1011001111010,而不是正确的结果1010001111010。

  这是一个在单一GF(2 ^ 8)整数域上实现该多项式乘法的Python函数。

  注意:这个函数(和下面的一些其他函数)使用了很多的位运算符,如>>和<<,因为它们在实现我们想做的事情上更快,更简洁。 这些操作符在大多数编程语言中都可用,它们并不仅限于Python,你可以在此处获得有关它们的更多信息

1 def cl_mul(x,y):
2     ‘‘‘Bitwise carry-less multiplication on integers‘‘‘
3     z = 0
4     i = 0
5     while (y>>i) > 0:
6         if y & (1<<i):
7             z ^= x<<i
8         i += 1
9     return z

  当然,结果不再是8位字节长度(本例中为13比特长),所以我们需要在完成之前再执行一个步骤。使用前面描述的长除法过程,将结果对数100011101(这个数字的选择在代码下面解释)取模。 在这种情况下,这被称为“模归约”,因为基本上我们所做的就是使用模数来做除法并保留余数。这在我们的例子中产生了最终答案11000011。

  1010001111010
^ 100011101
  0010110101010
  ^ 100011101
    00111011110
    ^ 100011101
      011000011

下面是用模归约来完成整个伽罗瓦域乘法的Python代码:
 1 def gf_mult_noLUT(x, y, prim=0):
 2     ‘‘‘Multiplication in Galois Fields without using a precomputed look-up table (and thus it‘s slower)
 3     by using the standard carry-less multiplication + modular reduction using an irreducible prime polynomial‘‘‘
 4 
 5     ### Define bitwise carry-less operations as inner functions ###
 6     def cl_mult(x,y):
 7         ‘‘‘Bitwise carry-less multiplication on integers‘‘‘
 8         z = 0
 9         i = 0
10         while (y>>i) > 0:
11             if y & (1<<i):
12                 z ^= x<<i
13             i += 1
14         return z
15 
16     def bit_length(n):
17         ‘‘‘Compute the position of the most significant bit (1) of an integer. Equivalent to int.bit_length()‘‘‘
18         bits = 0
19         while n >> bits: bits += 1
20         return bits
21 
22     def cl_div(dividend, divisor=None):
23         ‘‘‘Bitwise carry-less long division on integers and returns the remainder‘‘‘
24         # Compute the position of the most significant bit for each integers
25         dl1 = bit_length(dividend)
26         dl2 = bit_length(divisor)
27         # If the dividend is smaller than the divisor, just exit
28         if dl1 < dl2:
29             return dividend
30         # Else, align the most significant 1 of the divisor to the most significant 1 of the dividend (by shifting the divisor)
31         for i in range(dl1-dl2,-1,-1):
32             # Check that the dividend is divisible (useless for the first iteration but important for the next ones)
33             if dividend & (1 << i+dl2-1):
34                 # If divisible, then shift the divisor to align the most significant bits and XOR (carry-less subtraction)
35                 dividend ^= divisor << i
36         return dividend
37     
38     ### Main GF multiplication routine ###
39 
40     # Multiply the gf numbers
41     result = cl_mult(x,y)
42     # Then do a modular reduction (ie, remainder from the division) with an irreducible primitive polynomial so that it stays inside GF bounds
43     if prim > 0:
44         result = cl_div(result, prim)
45 
46     return result

结果:

>>> a = 0b10001001
>>> b = 0b00101010
>>> print bin(gf_mult_noLUT(a, b, 0)) # multiplication only
0b1010001111010
>>> print bin(gf_mult_noLUT(a, b, 0x11d)) # multiplication + modular reduction
0b11000011

  为什么模数取100011101(十六进制:0x11d)? 这里的数学有点复杂,但简而言之,100011101表示一个8次多项式,它是“不可约”的(意味着它不能被表示为两个较小多项式的乘积)。这个数字被称为本原多项式或不可约多项式,或质数多项式(我们将在本教程的其余部分主要使用后一个名称)。这是进行除法的必要条件,是为了保持伽罗华域的局限性,而不出现重复值。我们也可以选择其他的数值,但它们都基本相同,并且100011101(0x11d)是Reed-Solomon代码的共同质数多项式。如果您想知道如何生成这些质数,请参阅附录

  质数多项式的其他信息(可跳过):质数多项式是什么?它等同于一个质数,但是是在伽罗华域。请记住,伽罗华域使用2的倍数作为生成器。当然,标准算术中质数不能是2的倍数,但在伽罗华域中是可能的。为什么我们需要一个质数多项式? 因为如果要保持在域的界内,所以我们需要计算伽罗华域上的任何值的模。我们为什么不用伽罗华域的大小来模数? 因为我们最终会得到大量重复值,并且我们希望得到域中尽可能多的唯一值,所以当使用质数多项式进行模或XOR运算时,一个数总是有唯一的映射。

  感兴趣的读者请注意:作为使用聪明算法可以实现的一个例子,下面是以更简洁快速的方式实现GF域数值乘法的另一种方法,俄罗斯农民乘法算法:

 1 def gf_mult_noLUT(x, y, prim=0, field_charac_full=256, carryless=True):
 2     ‘‘‘Galois Field integer multiplication using Russian Peasant Multiplication algorithm (faster than the standard multiplication + modular reduction).
 3     If prim is 0 and carryless=False, then the function produces the result for a standard integers multiplication (no carry-less arithmetics nor modular reduction).‘‘‘
 4     r = 0
 5     while y: # while y is above 0
 6         if y & 1: r = r ^ x if carryless else r + x # y is odd, then add the corresponding x to r (the sum of all x‘s corresponding to odd y‘s will give the final product). Note that since we‘re in GF(2), the addition is in fact an XOR (very important because in GF(2) the multiplication and additions are carry-less, thus it changes the result!).
 7         y = y >> 1 # equivalent to y // 2
 8         x = x << 1 # equivalent to x*2
 9         if prim > 0 and x & field_charac_full: x = x ^ prim # GF modulo: if x >= 256 then apply modular reduction using the primitive polynomial (we just subtract, but since the primitive number can be above 256 then we directly XOR).
10 
11     return r

  请注意,使用带参数prim = 0和carryless = False的最后一个函数将返回标准整数乘法的结果(因此你可以看到非进位和进位加法之间的差异及其对乘法的影响)。

 

与对数的乘法

  上述过程并不是实现伽罗瓦域乘法最快的方法。乘以两个数最多需要八次迭代的乘法循环,然后是多达八次迭代的除法循环。但是,通过使用查找表,我们可以不用循环实现相乘。一种解决方案是在内存中构建整个乘法表,但这需要庞大的64k表。下面介绍的解决方案更加紧凑。

  首先,很容易注意到乘以2 = 00000010(按照惯例,这个数字被称为α或数字生成器):简单地左移一个位,然后如果需要的话与模数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

 

  如果这个表格以同样的方式继续下去,那么α的幂不会重复,直到α255= 00000001。因此,除零以外的每个元素都等于α的某个幂。我们定义的元素α被称为伽罗瓦域的基本元素或生成元素。

  这个结果提出了实现乘法的另一种方式:通过叠加α的指数。

  10001001 * 00101010 =α74*α142=α74+ 142 =α216= 11000011

  问题是,我们如何找到对应于10001001的α的倍数? 这被称为离散对数问题,没有有效的通用解决方案。但是,由于该域中只有256个元素,我们可以轻松构建一个对数表。至此,逆对数(指数)的相应表格也是有用的。关于这个技巧的更多数学信息可以在这里找到

 1 gf_exp = [0] * 512 # Create list of 512 elements. In Python 2.6+, consider using bytearray
 2 gf_log = [0] * 256
 3 
 4 def init_tables(prim=0x11d):
 5     ‘‘‘Precompute the logarithm and anti-log tables for faster computation later, using the provided primitive polynomial.‘‘‘
 6     # prim is the primitive (binary) polynomial. Since it‘s a polynomial in the binary sense,
 7     # it‘s only in fact a single galois field value between 0 and 255, and not a list of gf values.
 8     global gf_exp, gf_log
 9     gf_exp = [0] * 512 # anti-log (exponential) table
10     gf_log = [0] * 256 # log table
11     # For each possible value in the galois field 2^8, we will pre-compute the logarithm and anti-logarithm (exponential) of this value
12     x = 1
13     for i in range(0, 255):
14         gf_exp[i] = x # compute anti-log for this value and store it in a table
15         gf_log[x] = i # compute log at the same time
16         x = gf_mult_noLUT(x, 2, prim)
17 
18         # If you use only generator==2 or a power of 2, you can use the following which is faster than gf_mult_noLUT():
19         #x <<= 1 # multiply by 2 (change 1 by another number y to multiply by a power of 2^y)
20         #if x & 0x100: # similar to x >= 256, but a lot faster (because 0x100 == 256)
21             #x ^= prim # substract the primary polynomial to the current value (instead of 255, so that we get a unique set made of coprime numbers), this is the core of the tables generation
22 
23     # Optimization: double the size of the anti-log table so that we don‘t need to mod 255 to
24     # stay inside the bounds (because we will mainly use this table for the multiplication of two GF numbers, no more).
25     for i in range(255, 512):
26         gf_exp[i] = gf_exp[i - 255]
27     return [gf_log, gf_exp]

  Python注释:范围运算符的上限是唯一的,所以gf_exp [255]没有被上面的代码重复设置。

  gf_exp表大于其实际需要以简化乘法函数。这样,我们不必检查以确保gf_log [x] + gf_log [y]在表格大小内。

1 def gf_mul(x,y):
2     if x==0 or y==0:
3         return 0
4     return gf_exp[gf_log[x] + gf_log[y]] # should be gf_exp[(gf_log[x]+gf_log[y])%255] if gf_exp wasn‘t oversized

 

除法

  对数表方法的另一个优点是它允许我们使用对数的差来定义除法。在下面的代码中,添加了255以确保差不是负数。

1 def gf_div(x,y):
2     if y==0:
3         raise ZeroDivisionError()
4     if x==0:
5         return 0
6     return gf_exp[(gf_log[x] + 255 - gf_log[y]) % 255]

Python说明:raise语句抛出异常并中止gf_div函数的执行。

  有了这个除法的定义,对于任何x和任何非零y,gf_div(gf_mul(x,y),y)== x。

  高级程序员读者可能会发现写封装伽罗瓦域算法的类很有趣。可以通过运算符重载用熟悉的运算符 * 和 / 来代替对gf_mul和gf_div的调用,但是这会混淆正在执行的运算类型。某些细节可能的方式,将使类更广泛的实用推广。 例如,阿兹特克码使用五个不同的伽罗瓦域,元素大小为4到12位。

 

幂与逆

对数表方法在进行幂和逆运算时将再次简化和加速我们的计算:

1 def gf_pow(x, power):
2     return gf_exp[(gf_log[x] * power) % 255]
3 
4 def gf_inverse(x):
5     return gf_exp[255 - gf_log[x]] # gf_inverse(x) == gf_div(1, x)

 

多项式

  在进入Reed-Solomon码之前,我们需要定义几个系数为伽罗瓦域元素的多项式的运算。这是一个潜在的混淆根源,因为元素本身被描述为多项式;我的建议不是想太多。令人困惑的是x仍然被用作占位符。这个x与前面提到的x没有任何关系,所以把它们混淆起来。

  之前用于伽罗瓦域元素的二进制表示法在这一点上开始变得不方便,所以我将改为十六进制表示。

00000001 x4 + 00001111 x3 + 00110110 x2 + 01111000 x + 01000000 = 01 x4 + 0f x3 + 36 x2 + 78 x + 40

  在Python中,多项式由x的幂的递减顺序的一列数字表示,所以上述多项式变为[0x01, 0x0f, 0x36, 0x78, 0x40]。(也可以用相反的顺序;两种选择都有其优缺点。)

  第一个函数将多项式乘以标量。

1 def gf_poly_scale(p,x):
2     r = [0] * len(p)
3     for i in range(0, len(p)):
4         r[i] = gf_mul(p[i], x)
5     return r

  Python程序员注意:这个函数不是用“pythonic”风格编写的。它可以被优雅地表示为一个递推式构造列表,但是我将限制于更易于翻译为其他编程语言的语言特性上。

  函数 “adds” 对两个多项式做加法(使用异或)。

1 def gf_poly_add(p,q):
2     r = [0] * max(len(p),len(q))
3     for i in range(0,len(p)):
4         r[i+len(r)-len(p)] = p[i]
5     for i in range(0,len(q)):
6         r[i+len(r)-len(q)] ^= q[i]
7     return r

  下一个函数将两个多项式相乘。

 1 def gf_poly_mul(p,q):
 2     ‘‘‘Multiply two polynomials, inside Galois Field‘‘‘
 3     # Pre-allocate the result array
 4     r = [0] * (len(p)+len(q)-1)
 5     # Compute the polynomial multiplication (just like the outer product of two vectors,
 6     # we multiply each coefficients of p with all coefficients of q)
 7     for j in range(0, len(q)):
 8         for i in range(0, len(p)):
 9             r[i+j] ^= gf_mul(p[i], q[j]) # equivalent to: r[i + j] = gf_add(r[i+j], gf_mul(p[i], q[j]))
10                                                          # -- you can see it‘s your usual polynomial multiplication
11     return r

  最后,我们需要一个函数来计算多项式在特定值x的值,得到标量结果。使用霍纳方法(Horner‘s method)以避免明确计算x的幂。霍纳的方法通过连续地分解各项,以便我们总是以迭代方式处理x ^ 1,以避免计算更高阶的项:

01 x4 + 0f x3 + 36 x2 + 78 x + 40 = (((01 x + 0fx + 36x + 78x + 40

1 def gf_poly_eval(poly, x):
2     ‘‘‘Evaluates a polynomial in GF(2^p) given the value for x. This is based on Horner‘s scheme for maximum efficiency.‘‘‘
3     y = poly[0]
4     for i in range(1, len(poly)):
5         y = gf_mul(y, x) ^ poly[i]
6     return y

  还有一个我们需要的多项式操作:多项式除法。这比多项式的其他运算更复杂,所以我们将在下一章与里德-所罗门编码一起研究它。

 

里德-所罗门编码

  现在预备知识已经完成,我们已经准备好开始研究里德-所罗门编码。

洞察编码理论

  但首先,为什么我们必须学习有限域和多项式?因为这是里德-所罗门等纠错遍码的主要思想:不是将消息看作一系列(ASCII)数字,而是将其视为一个遵循明确定义的有限域运算法则的多项式。换句话说,通过使用多项式和有限域算术表示数据,我们为数据添加了一个结构。消息的值仍然是相同的,但是这个概念结构现在允许我们使用定义良好的数学规则在消息上对字符值进行操作。这种结构,我们总是知道,因为它在外部并且独立于数据,这使我们能够修复损坏的消息。

  因此,即使在你的代码实现中,你可能选择非显式表示多项式和有限域算术,但这些概念对于纠错编码产生作用是必不可少的,你会发现这些概念是隐含的(即隐式地)实现。

  现在我们将把这些想法付诸实践!

 

RS编码

编码大纲

  像BCH码一样,里德-所罗门码通过用表示消息的多项式除以不可约生成多项式来编码,余数的就是RS码,我们将其附加到原始消息。

  为什么?我们之前说过,BCH码和其他大多数纠错码背后的原理是使用一个单词差异很大的简化字典,从而使单词之间的距离最大化,并且更长的单词具有更大的距离:这里的原理相同,因为我们用附加码元(余数)来加长原始信息,这增加了距离;其次是因为余数几乎是唯一的(得益于精心设计的不可约生成多项式),所以它可以被聪明的算法用来推导原始消息。

  总结一下,与加密类似:生成多项式是编码字典,多项式除法是用来将我们的消息通过字典(生成多项式)转换成RS码的运算。

 

异常管理

  为了管理我们无法更正消息的错误和情况,我们将通过引发异常来显示有意义的错误消息。 我们将制作我们自己的自定义异常,以便用户可以轻松捕捉并管理它们:

1 class ReedSolomonError(Exception):
2     pass

  要通过抛出自定义异常来显示错误,可以简单地执行以下操作:

1 raise ReedSolomonError("Some error message")

  通过使用try / except块,可以轻松捕获此异常以进行管理:

1 try:
2     raise ReedSolomonError("Some error message")
3 except ReedSolomonError, e:
4     pass # do something here

 

RS生成多项式

  里德-所罗门码使用与BCH码类似的生成多项式(不要与生成器数alpha混淆)。生成器是因子(x - αn)的乘积,从QR码的n = 0开始。 例如:

  g4(x) = (x - α0) (x - α1) (x - α2) (x - α3) = 01 x4 + 0f x3 + 36 x2 + 78 x + 40

  这是一个计算给定数量的纠错符号的生成多项式的函数。

1 def rs_generator_poly(nsym):
2     ‘‘‘Generate an irreducible generator polynomial (necessary to encode a message into Reed-Solomon)‘‘‘
3     g = [1]
4     for i in range(0, nsym):
5         g = gf_poly_mul(g, [1, gf_pow(2, i)])
6     return g

  这个函数有点低效,因为它为g分配了连续的更大的数组。尽管这在实践中不太可能是一个性能问题,但是读者可能会发现重写优化器使得g只被分配一次,或者可以计算一次并记住g,因为它对于给定的nsym是固定的, 你可以重用g。

 

多项式除法

  存在几种用于多项式分割的算法,在高中经常教导的最简单的算法是长分割。 这个例子显示了消息12 34 56的计算。

                             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

  注意:起初看起来很复杂,但当你用基数为10的整数进行多项式长整除时,它的工作方式与通常的方式完全相同,只不过它是在16(十六进制)中。 如果你在理解这个例子时遇到了一些麻烦,你可以尝试将数字转换为10进制,你很快就会看到结果是当然的。

  将其余数与消息相连,所以编码后的消息是12 34 56 37 e6 78 d9。

  但是,长除法很慢,因为它需要进行很多次递归迭代才结束。可以设计更有效的策略,例如使用合成除法(也称为霍纳方法[Horner‘s method],可以在可汗学院上找到一个很好的教学视频)。这是一个实现GF(2 ^ p)多项式的扩展合成除法的函数(扩展是因为除数是一个多项式而不是单项):

 1 def gf_poly_div(dividend, divisor):
 2     ‘‘‘Fast polynomial division by using Extended Synthetic Division and optimized for GF(2^p) computations
 3     (doesn‘t work with standard polynomials outside of this galois field, see the Wikipedia article for generic algorithm).‘‘‘
 4     # CAUTION: this function expects polynomials to follow the opposite convention at decoding:
 5     # the terms must go from the biggest to lowest degree (while most other functions here expect
 6     # a list from lowest to biggest degree). eg: 1 + 2x + 5x^2 = [5, 2, 1], NOT [1, 2, 5]
 7 
 8     msg_out = list(dividend) # Copy the dividend
 9     #normalizer = divisor[0] # precomputing for performance
10     for i in range(0, len(dividend) - (len(divisor)-1)):
11         #msg_out[i] /= normalizer # for general polynomial division (when polynomials are non-monic), the usual way of using
12                                   # synthetic division is to divide the divisor g(x) with its leading coefficient, but not needed here.
13         coef = msg_out[i] # precaching
14         if coef != 0: # log(0) is undefined, so we need to avoid that case explicitly (and it‘s also a good optimization).
15             for j in range(1, len(divisor)): # in synthetic division, we always skip the first coefficient of the divisior,
16                                               # because it‘s only used to normalize the dividend coefficient
17                 if divisor[j] != 0: # log(0) is undefined
18                     msg_out[i + j] ^= gf_mul(divisor[j], coef) # equivalent to the more mathematically correct
19                                                                # (but xoring directly is faster): msg_out[i + j] += -divisor[j] * coef
20 
21     # The resulting msg_out contains both the quotient and the remainder, the remainder being the size of the divisor
22     # (the remainder has necessarily the same degree as the divisor -- not length but degree == length-1 -- since it‘s
23     # what we couldn‘t divide from the dividend), so we compute the index where this separation is, and return the quotient and remainder.
24     separator = -(len(divisor)-1)
25     return msg_out[:separator], msg_out[separator:] # return quotient, remainder.

 

编码主函数

  现在,下面是如何编码一条消息以获取其RS码:

 1 def rs_encode_msg(msg_in, nsym):
 2     ‘‘‘Reed-Solomon main encoding function‘‘‘
 3     gen = rs_generator_poly(nsym)
 4 
 5     # Pad the message, then divide it by the irreducible generator polynomial
 6     _, remainder = gf_poly_div(msg_in + [0] * (len(gen)-1), gen)
 7     # The remainder is our RS code! Just append it to our original message to get our full codeword (this represents a polynomial of max 256 terms)
 8     msg_out = msg_in + remainder
 9     # Return the codeword
10     return msg_out

  很简单,不是吗? 编码实际上是里德-所罗门最简单的部分,并且它始终是用相同的方法(多项式除法)。解码是里德-所罗门的难题,你会发现很多不同的算法,这取决于你的需求,但我们稍后会谈到。

  这个函数非常快速,但是由于编码非常关键,下面是一个增强的编码函数,它内联 [inline] 多项式合成除法,这是里德-所罗门软件库中最常见的形式:

 1 def rs_encode_msg(msg_in, nsym):
 2     ‘‘‘Reed-Solomon main encoding function, using polynomial division (algorithm Extended Synthetic Division)‘‘‘
 3     if (len(msg_in) + nsym) > 255: raise ValueError("Message is too long (%i when max is 255)" % (len(msg_in)+nsym))
 4     gen = rs_generator_poly(nsym)
 5     # Init msg_out with the values inside msg_in and pad with len(gen)-1 bytes (which is the number of ecc symbols).
 6     msg_out = [0] * (len(msg_in) + len(gen)-1)
 7     # Initializing the Synthetic Division with the dividend (= input message polynomial)
 8     msg_out[:len(msg_in)] = msg_in
 9 
10     # Synthetic division main loop
11     for i in range(len(msg_in)):
12         # Note that it‘s msg_out here, not msg_in. Thus, we reuse the updated value at each iteration
13         # (this is how Synthetic Division works: instead of storing in a temporary register the intermediate values,
14         # we directly commit them to the output).
15         coef = msg_out[i]
16 
17         # log(0) is undefined, so we need to manually check for this case. There‘s no need to check
18         # the divisor here because we know it can‘t be 0 since we generated it.
19         if coef != 0:
20             # in synthetic division, we always skip the first coefficient of the divisior, because it‘s only used to normalize the dividend coefficient (which is here useless since the divisor, the generator polynomial, is always monic)
21             for j in range(1, len(gen)):
22                 msg_out[i+j] ^= gf_mul(gen[j], coef) # equivalent to msg_out[i+j] += gf_mul(gen[j], coef)
23 
24     # At this point, the Extended Synthetic Divison is done, msg_out contains the quotient in msg_out[:len(msg_in)]
25     # and the remainder in msg_out[len(msg_in):]. Here for RS encoding, we don‘t need the quotient but only the remainder
26     # (which represents the RS code), so we can just overwrite the quotient with the input message, so that we get
27     # our complete codeword composed of the message + code.
28     msg_out[:len(msg_in)] = msg_in
29 
30     return msg_out

  该算法速度更快,但对实际应用中仍然非常慢,特别是在Python中。有一些方法可以使用各种技巧来优化速度,例如通过内联(该操作替代而gf_mul以避免函数调用);通过预先计算(gen和coef的对数,甚至通过生成一个乘法表 - 但似乎后者在Python中不能很好地工作);通过使用静态类型的构造(将gf_log和gf_exp赋值给 array.array(‘i‘,[...]);通过使用内存视图(就像更改所有列表为bytearrays);通过使用PyPy运行;或者通过将算法转换成Cython或C扩展[1]。

  此示例显示了前面介绍的示例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版本说明:print函数的语法已更改,本示例使用Python 3.0+版本。在之前版本的Python(特别是Python 2.x)中,用print hex(msg [i])(包括最终逗号)和range by xrange替换打印行。

 

RS解码

解码大纲

  里德-所罗门解码是从潜在损坏的消息和RS码中返回纠正的消息的过程。换句话说,解码是使用先前计算的RS码修复信息的过程。

  虽然里德-所罗门只有一种编码消息的方法,但却有很多不同的解码方法,因此存在许多不同的解码算法。

  但是,我们通常可以通过5个步骤概述解码过程[2] [3]

  1. 计算症状多项式。 这使我们能够使用Berlekamp-Massey(或其他算法)分析哪些字符错误,并且还可以快速检查输入消息是否完全被破坏。

  2. 计算丢失/错误定位多项式(利用证据)。 这是由Berlekamp-Massey算法计算出来的,它是一个能够准确告诉我们哪些字符被破坏的检测器。

  3. 计算丢失/错误评估多项式(通过症状和丢失/错误定位多项式)。对字符被篡改的程度的评估是必要的(有助于计算幅度)。

  4. 计算丢失/错误幅度多项式(从上述所有3个多项式):该多项式也可以称为错误多项式,因为实际上它准确地存储了需要从接收到的消息中减去的值,以获得原始正确的消息( 即用正确的丢失字符值)。 换句话说,在此刻,我们提取噪声并将其存储在这个多项式中,我们只需要从输入消息中去除噪声来修复它。

  5. 简单地通过从输入消息中减去幅度多项式来修复输入消息

  我们将在下面描述这五个步骤中的每一个。

  此外,解码器还可以根据它们可以修复的错误类型进行分类:丢失(我们知道被破坏的字符的位置,但不是大小),错误(我们忽略位置和大小),或者同时发生了错误和丢失。 我们将描述如何支持所有这些类型。

 

症状的计算

  对里德-所罗门消息进行解码需要几个步骤。第一步是计算消息的“症状”。将消息视为多项式,并在α0,α1,α2,...,αn处对其进行估值。由于这些是生成多项式的零,如果扫描的消息未损坏,结果应为零(可用于在纠正了损坏的消息并完全修复后之后检查消息是否损坏)。 如果不是,则症状包含确定应该进行校正的所有必要信息。编写函数来计算症状很简单。

 1 def rs_calc_syndromes(msg, nsym):
 2     ‘‘‘Given the received codeword msg and the number of error correcting symbols (nsym), computes the syndromes polynomial.
 3     Mathematically, it‘s essentially equivalent to a Fourrier Transform (Chien search being the inverse).
 4     ‘‘‘
 5     # Note the "[0] +" : we add a 0 coefficient for the lowest degree (the constant). This effectively shifts the syndrome, and will shift every computations depending on the syndromes (such as the errors locator polynomial, errors evaluator polynomial, etc. but not the errors positions).
 6     # This is not necessary, you can adapt subsequent computations to start from 0 instead of skipping the first iteration (ie, the often seen range(1, n-k+1)),
 7     synd = [0] * nsym
 8     for i in range(0, nsym):
 9         synd[i] = gf_poly_eval(msg, gf_pow(2,i))
10     return [0] + synd # pad with one 0 for mathematical precision (else we can end up with weird calculations sometimes)

  继续这个例子,我们看到原始码字没有任何损坏的症状确实是零。在消息或其RS代码中引入至少一个字符的损坏会产生非零的症状。

>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0] # not corrupted message = all 0 syndromes
>>> 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] # when corrupted, the syndromes will be non zero

  以下是自动执行此检查的代码:

1 def rs_check(msg, nsym):
2     ‘‘‘Returns true if the message + ecc has no error of false otherwise (may not always catch a wrong decoding or a wrong message, particularly if there are too many errors -- above the Singleton bound --, but it usually does)‘‘‘
3     return ( max(rs_calc_syndromes(msg, nsym)) == 0 )

 

丢失数据的恢复

  如果错误的位置已知,纠正代码中的错误是最简单的。这被称为丢失纠正。可以对添加到代码中的每个纠错符号校正一个丢失符号(即字符)。如果错误位置未知,每个符号错误需要两个EC符号(这样您可以纠正两次以内的错误)。这使得丢失校正在实践中很有用,如果被扫描的QR码的一部分被覆盖或被物理撕掉。 但扫描仪可能很难确定是否发生了这种情况,所以不是所有的QR码扫描仪都可以执行丢失校正。

  现在我们已经有了症状,我们需要计算定位多项式。 这很容易:

 1 def rs_find_errata_locator(e_pos):
 2     ‘‘‘Compute the erasures/errors/errata locator polynomial from the erasures/errors/errata positions
 3        (the positions must be relative to the x coefficient, eg: "hello worldxxxxxxxxx" is tampered to "h_ll_ worldxxxxxxxxx"
 4        with xxxxxxxxx being the ecc of length n-k=9, here the string positions are [1, 4], but the coefficients are reversed
 5        since the ecc characters are placed as the first coefficients of the polynomial, thus the coefficients of the
 6        erased characters are n-1 - [1, 4] = [18, 15] = erasures_loc to be specified as an argument.‘‘‘
 7 
 8     e_loc = [1] # just to init because we will multiply, so it must be 1 so that the multiplication starts correctly without nulling any term
 9     # erasures_loc = product(1 - x*alpha**i) for i in erasures_pos and where alpha is the alpha chosen to evaluate polynomials.
10     for i in e_pos:
11         e_loc = gf_poly_mul( e_loc, gf_poly_add([1], [gf_pow(2, i), 0]) )
12     return e_loc

  接下来,从定位多项式中计算丢失/错误评估多项式很容易,它只是先做一个多项式乘法,然后再做一个多项式除法(你可以将其替换为列表分部计算,因为这是我们最终需要的效果):

 1 def rs_find_error_evaluator(synd, err_loc, nsym):
 2     ‘‘‘Compute the error (or erasures if you supply sigma=erasures locator polynomial, or errata) evaluator polynomial Omega
 3        from the syndrome and the error/erasures/errata locator Sigma.‘‘‘
 4 
 5     # Omega(x) = [ Synd(x) * Error_loc(x) ] mod x^(n-k+1)
 6     _, remainder = gf_poly_div( gf_poly_mul(synd, err_loc), ([1] + [0]*(nsym+1)) ) # first multiply syndromes * errata_locator, then do a
 7                                                                                    # polynomial division to truncate the polynomial to the
 8                                                                                    # required length
 9 
10     # Faster way that is equivalent
11     #remainder = gf_poly_mul(synd, err_loc) # first multiply the syndromes with the errata locator polynomial
12     #remainder = remainder[len(remainder)-(nsym+1):] # then slice the list to truncate it (which represents the polynomial), which
13                                                                           # is equivalent to dividing by a polynomial of the length we want
14 
15     return remainder

最后,Forney算法用于计算校正值(也称为误差幅度多项式)。 它在下面的函数中实现。

 1 def rs_correct_errata(msg_in, synd, err_pos): # err_pos is a list of the positions of the errors/erasures/errata
 2     ‘‘‘Forney algorithm, computes the values (error magnitude) to correct the input message.‘‘‘
 3     # calculate errata locator polynomial to correct both errors and erasures (by combining the errors positions given by the error locator polynomial found by BM with the erasures positions given by caller)
 4     coef_pos = [len(msg_in) - 1 - p for p in err_pos] # need to convert the positions to coefficients degrees for the errata locator algo to work (eg: instead of [0, 1, 2] it will become [len(msg)-1, len(msg)-2, len(msg) -3])
 5     err_loc = rs_find_errata_locator(coef_pos)
 6     # calculate errata evaluator polynomial (often called Omega or Gamma in academic papers)
 7     err_eval = rs_find_error_evaluator(synd[::-1], err_loc, len(err_loc)-1)[::-1]
 8 
 9     # Second part of Chien search to get the error location polynomial X from the error positions in err_pos (the roots of the error locator polynomial, ie, where it evaluates to 0)
10     X = [] # will store the position of the errors
11     for i in range(0, len(coef_pos)):
12         l = 255 - coef_pos[i]
13         X.append( gf_pow(2, -l) )
14 
15     # Forney algorithm: compute the magnitudes
16     E = [0] * (len(msg_in)) # will store the values that need to be corrected (substracted) to the message containing errors. This is sometimes called the error magnitude polynomial.
17     Xlength = len(X)
18     for i, Xi in enumerate(X):
19 
20         Xi_inv = gf_inverse(Xi)
21 
22         # Compute the formal derivative of the error locator polynomial (see Blahut, Algebraic codes for data transmission, pp 196-197).
23         # the formal derivative of the errata locator is used as the denominator of the Forney Algorithm, which simply says that the ith error value is given by error_evaluator(gf_inverse(Xi)) / error_locator_derivative(gf_inverse(Xi)). See Blahut, Algebraic codes for data transmission, pp 196-197.
24         err_loc_prime_tmp = []
25         for j in range(0, Xlength):
26             if j != i:
27                 err_loc_prime_tmp.append( gf_sub(1, gf_mul(Xi_inv, X[j])) )
28         # compute the product, which is the denominator of the Forney algorithm (errata locator derivative)
29         err_loc_prime = 1
30         for coef in err_loc_prime_tmp:
31             err_loc_prime = gf_mul(err_loc_prime, coef)
32         # equivalent to: err_loc_prime = functools.reduce(gf_mul, err_loc_prime_tmp, 1)
33 
34         # Compute y (evaluation of the errata evaluator polynomial)
35         # This is a more faithful translation of the theoretical equation contrary to the old forney method. Here it is an exact reproduction:
36         # Yl = omega(Xl.inverse()) / prod(1 - Xj*Xl.inverse()) for j in len(X)
37         y = gf_poly_eval(err_eval[::-1], Xi_inv) # numerator of the Forney algorithm (errata evaluator evaluated)
38         y = gf_mul(gf_pow(Xi, 1), y)
39         
40         # Compute the magnitude
41         magnitude = gf_div(y, err_loc_prime) # magnitude value of the error, calculated by the Forney algorithm (an equation in fact): dividing the errata evaluator with the errata locator derivative gives us the errata magnitude (ie, value to repair) the ith symbol
42         E[err_pos[i]] = magnitude # store the magnitude for this error into the magnitude polynomial
43 
44     # Apply the correction of values to get our message corrected! (note that the ecc bytes also gets corrected!)
45     # (this isn‘t the Forney algorithm, we just apply the result of decoding here)
46     msg_in = gf_poly_add(msg_in, E) # equivalent to Ci = Ri - Ei where Ci is the correct message, Ri the received (senseword) message, and Ei the errata magnitudes (minus is replaced by XOR since it‘s equivalent in GF(2^p)). So in fact here we substract from the received message the errors magnitude, which logically corrects the value to what it should be.
47     return msg_in

数学注释:误差值表达式的分母是误差定位多项式q的形式导数。这通过用n cn xn-1替换每个项cn xn的通常程序来计算。由于我们在特征二的领域工作,当n是奇数时,cn等于cn,当n是偶数时,等于0。因此,我们可以简单地移除偶数系数(导致多项式qprime)并评估qprime(x2)。

Python注释:该函数使用[:: - 1]来反转列表中元素的顺序。这是必要的,因为这些函数并不都使用相同的顺序约定(即,一些首先使用最少的项目,另一些首先使用最大的项目)。它还使用列表理解,这只是一种简洁的方式来编写for循环,其中项目被附加到列表中,但Python解释器可以优化这一点而不是循环。

继续这个例子,这里我们使用rs_correct_errata来恢复消息的第一个字节。

>>> msg[0] = 0
>>> synd = rs_calc_syndromes(msg, 10)
>>> msg = rs_correct_errata(msg, synd, [0]) # [0] is the list of the erasures locations, here it‘s the first character, at position 0
>>> print(hex(msg[0]))
0x40

 

错误数据的纠正

  在错误位置未知的情况下(我们通常称之为错误,与位置已知的丢失相对应),我们将使用与处理丢失数据相同的步骤,但现在我们需要额外的步骤来查找位置。Berlekamp-Massey算法用于计算错误定位多项式,我们稍后可以使用它来确定错误位置:

 1 def rs_find_error_locator(synd, nsym, erase_loc=None, erase_count=0):
 2     ‘‘‘Find error/errata locator and evaluator polynomials with Berlekamp-Massey algorithm‘‘‘
 3     # The idea is that BM will iteratively estimate the error locator polynomial.
 4     # To do this, it will compute a Discrepancy term called Delta, which will tell us if the error locator polynomial needs an update or not
 5     # (hence why it‘s called discrepancy: it tells us when we are getting off board from the correct value).
 6 
 7     # Init the polynomials
 8     if erase_loc: # if the erasure locator polynomial is supplied, we init with its value, so that we include erasures in the final locator polynomial
 9         err_loc = list(erase_loc)
10         old_loc = list(erase_loc)
11     else:
12         err_loc = [1] # This is the main variable we want to fill, also called Sigma in other notations or more formally the errors/errata locator polynomial.
13         old_loc = [1] # BM is an iterative algorithm, and we need the errata locator polynomial of the previous iteration in order to update other necessary variables.
14     #L = 0 # update flag variable, not needed here because we use an alternative equivalent way of checking if update is needed (but using the flag could potentially be faster depending on if using length(list) is taking linear time in your language, here in Python it‘s constant so it‘s as fast.
15 
16     # 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.
17     # 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, which set to 0 the coefficients corresponding to erasures, leaving only the coefficients corresponding to errors.
18     synd_shift = 0
19     if len(synd) > nsym: synd_shift = len(synd) - nsym
20 
21     for i in range(0, nsym-erase_count): # generally: nsym-erase_count == len(synd), except when you input a partial erase_loc and using the full syndrome instead of the Forney syndrome, in which case nsym-erase_count is more correct (len(synd) will fail badly with IndexError).
22         if erase_loc: # if an erasures locator polynomial was provided to init the errors locator polynomial, then we must skip the FIRST erase_count iterations (not the last iterations, this is very important!)
23             K = erase_count+i+synd_shift
24         else: # if erasures locator is not provided, then either there‘s no erasures to account or we use the Forney syndromes, so we don‘t need to use erase_count nor erase_loc (the erasures have been trimmed out of the Forney syndromes).
25             K = i+synd_shift
26 
27         # Compute the discrepancy Delta
28         # Here is the close-to-the-books operation to compute the discrepancy Delta: it‘s a simple polynomial multiplication of error locator with the syndromes, and then we get the Kth element.
29         #delta = gf_poly_mul(err_loc[::-1], synd)[K] # theoretically it should be gf_poly_add(synd[::-1], [1])[::-1] instead of just synd, but it seems it‘s not absolutely necessary to correctly decode.
30         # But this can be optimized: since we only need the Kth element, we don‘t need to compute the polynomial multiplication for any other element but the Kth. Thus to optimize, we compute the polymul only at the item we need, skipping the rest (avoiding a nested loop, thus we are linear time instead of quadratic).
31         # This optimization is actually described in several figures of the book "Algebraic codes for data transmission", Blahut, Richard E., 2003, Cambridge university press.
32         delta = synd[K]
33         for j in range(1, len(err_loc)):
34             delta ^= gf_mul(err_loc[-(j+1)], synd[K - j]) # delta is also called discrepancy. Here we do a partial polynomial multiplication (ie, we compute the polynomial multiplication only for the term of degree K). Should be equivalent to brownanrs.polynomial.mul_at().
35         #print "delta", K, delta, list(gf_poly_mul(err_loc[::-1], synd)) # debugline
36 
37         # Shift polynomials to compute the next degree
38         old_loc = old_loc + [0]
39 
40         # Iteratively estimate the errata locator and evaluator polynomials
41         if delta != 0: # Update only if there‘s a discrepancy
42             if len(old_loc) > len(err_loc): # Rule B (rule A is implicitly defined because rule A just says that we skip any modification for this iteration)
43             #if 2*L <= K+erase_count: # equivalent to len(old_loc) > len(err_loc), as long as L is correctly computed
44                 # Computing errata locator polynomial Sigma
45                 new_loc = gf_poly_scale(old_loc, delta)
46                 old_loc = gf_poly_scale(err_loc, gf_inverse(delta)) # effectively we are doing err_loc * 1/delta = err_loc // delta
47                 err_loc = new_loc
48                 # Update the update flag
49                 #L = K - L # the update flag L is tricky: in Blahut‘s schema, it‘s mandatory to use `L = K - L - erase_count` (and indeed in a previous draft of this function, if you forgot to do `- erase_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 `- erase_count` to update L on your own implementation and see which one works OK without producing wrong decoding failures.
50 
51             # Update with the discrepancy
52             err_loc = gf_poly_add(err_loc, gf_poly_scale(old_loc, delta))
53 
54     # Check if the result is correct, that there‘s not too many errors to correct
55     while len(err_loc) and err_loc[0] == 0: del err_loc[0] # drop leading 0s, else errs will not be of the correct size
56     errs = len(err_loc) - 1
57     if (errs-erase_count) * 2 + erase_count > nsym:
58         raise ReedSolomonError("Too many errors to correct")    # too many errors to correct
59 
60     return err_loc

  然后,使用错误定位器多项式,我们简单地使用称为试验替代的蛮力方法来找出该多项式的零点,该零点标识错误位置(即,需要校正的字符的索引)。 存在更高效的Chien搜索算法,该算法避免了在每个迭代步骤重新计算整个估值,但该算法仅作为练习给读者。

 1 def rs_find_errors(err_loc, nmess): # nmess is len(msg_in)
 2     ‘‘‘Find the roots (ie, where evaluation = zero) of error polynomial by brute-force trial, this is a sort of Chien‘s search
 3     (but less efficient, Chien‘s search is a way to evaluate the polynomial such that each evaluation only takes constant time).‘‘‘
 4     errs = len(err_loc) - 1
 5     err_pos = []
 6     for i in range(nmess): # normally we should try all 2^8 possible values, but here we optimize to just check the interesting symbols
 7         if gf_poly_eval(err_loc, gf_pow(2, i)) == 0: # It‘s a 0? Bingo, it‘s a root of the error locator polynomial,
 8                                                                                    # in other terms this is the location of an error
 9             err_pos.append(nmess - 1 - i)
10     # Sanity check: the number of errors/errata positions found should be exactly the same as the length of the errata locator polynomial
11     if len(err_pos) != errs:
12         # couldn‘t find error locations
13         raise ReedSolomonError("Too many (or few) errors found by Chien Search for the errata locator polynomial!")
14     return err_pos

数学注释:当错误定位多项式是线性的(err_poly具有长度2)时,可以很容易地解决,而不诉诸强力方法。上面介绍的功能没有利用这个事实,但感兴趣的读者可能希望实施更有效的解决方案。类似地,当误差定位器是二次的时,它可以通过使用二次公式的推广来解决。一个更有雄心的读者可能会希望实现这个程序。

  以下是消息中三个错误得到纠正的示例:

>>> print(hex(msg[10]))
0x96
>>> msg[0] = 6
>>> msg[10] = 7
>>> msg[20] = 8
>>> synd = rs_calc_syndromes(msg, 10)
>>> err_loc = rs_find_error_locator(synd, nsym)
>>> pos = rs_find_errors(err_loc[::-1], len(msg)) # find the errors locations
>>> print(pos)
[20, 10, 0]
>>> msg = rs_correct_errata(msg, synd, pos)
>>> print(hex(msg[10]))
0x96

 

错误和丢失数据的纠正

  里德-所罗门解码器可以同时解码丢失和错误,受限于2 * e + v <=(n - k)(称为单例束缚),其中e是错误的数量, v是丢失的数量,(n - k)是RS码字符的数量(在代码中称为nsym)。基本上,这意味着对于每一个数据丢失,只需要一个RS码字符来修复;而对于每个错误,需要两个RS码字符(因为除了估值/幅度以外,还需要找到的错误位置以纠正)。这种解码器被称为错误和丢失解码器,或勘误解码器。

  为了纠正错误和丢失,我们必须防止数据丢失干扰错误定位过程。可以通过计算Forney症状来完成,如下所示。

 1 def rs_forney_syndromes(synd, pos, nmess):
 2     # Compute Forney syndromes, which computes a modified syndromes to compute only errors (erasures are trimmed out). Do not confuse this with Forney algorithm, which allows to correct the message based on the location of errors.
 3     erase_pos_reversed = [nmess-1-p for p in pos] # prepare the coefficient degree positions (instead of the erasures positions)
 4 
 5     # Optimized method, all operations are inlined
 6     fsynd = list(synd[1:])      # make a copy and trim the first coefficient which is always 0 by definition
 7     for i in range(0, len(pos)):
 8         x = gf_pow(2, erase_pos_reversed[i])
 9         for j in range(0, len(fsynd) - 1):
10             fsynd[j] = gf_mul(fsynd[j], x) ^ fsynd[j + 1]
11 
12     # Equivalent, theoretical way of computing the modified Forney syndromes: fsynd = (erase_loc * synd) % x^(n-k)
13     # See Shao, H. M., Truong, T. K., Deutsch, L. J., & Reed, I. S. (1986, April). A single chip VLSI Reed-Solomon decoder. In Acoustics, Speech, and Signal Processing, IEEE International Conference on ICASSP‘86. (Vol. 11, pp. 2151-2154). IEEE.ISO 690
14     #erase_loc = rs_find_errata_locator(erase_pos_reversed, generator=generator) # computing the erasures locator polynomial
15     #fsynd = gf_poly_mul(erase_loc[::-1], synd[1:]) # then multiply with the syndrome to get the untrimmed forney syndrome
16     #fsynd = fsynd[len(pos):] # then trim the first erase_pos coefficients which are useless. Seems to be not necessary, but this reduces the computation time later in BM (thus it‘s an optimization).
17 
18     return fsynd

  然后可以使用Forney症状代替错误定位过程中的常规症状。

  下面的函数rs_correct_msg将完整的过程结合在一起。 通过在消息msg_in(完整接收消息:原始消息+ ecc)中提供erase_pos(消息索引位置列表)来指示丢失。

 1 def rs_correct_msg(msg_in, nsym, erase_pos=None):
 2     ‘‘‘Reed-Solomon main decoding function‘‘‘
 3     if len(msg_in) > 255: # can‘t decode, message is too big
 4         raise ValueError("Message is too long (%i when max is 255)" % len(msg_in))
 5 
 6     msg_out = list(msg_in)     # copy of message
 7     # erasures: set them to null bytes for easier decoding (but this is not necessary, they will be corrected anyway, but debugging will be easier with null bytes because the error locator polynomial values will only depend on the errors locations, not their values)
 8     if erase_pos is None:
 9         erase_pos = []
10     else:
11         for e_pos in erase_pos:
12             msg_out[e_pos] = 0
13     # check if there are too many erasures to correct (beyond the Singleton bound)
14     if len(erase_pos) > nsym: raise ReedSolomonError("Too many erasures to correct")
15     # prepare the syndrome polynomial using only errors (ie: errors = characters that were either replaced by null byte
16     # or changed to another character, but we don‘t know their positions)
17     synd = rs_calc_syndromes(msg_out, nsym)
18     # check if there‘s any error/erasure in the input codeword. If not (all syndromes coefficients are 0), then just return the message as-is.
19     if max(synd) == 0:
20         return msg_out[:-nsym], msg_out[-nsym:]  # no errors
21 
22     # compute the Forney syndromes, which hide the erasures from the original syndrome (so that BM will just have to deal with errors, not erasures)
23     fsynd = rs_forney_syndromes(synd, erase_pos, len(msg_out))
24     # compute the error locator polynomial using Berlekamp-Massey
25     err_loc = rs_find_error_locator(fsynd, nsym, erase_count=len(erase_pos))
26     # locate the message errors using Chien search (or brute-force search)
27     err_pos = rs_find_errors(err_loc[::-1] , len(msg_out))
28     if err_pos is None:
29         raise ReedSolomonError("Could not locate error")    # error location failed
30 
31     # Find errors values and apply them to correct the message
32     # compute errata evaluator and errata magnitude polynomials, then correct errors and erasures
33     msg_out = rs_correct_errata(msg_out, synd, (erase_pos + err_pos)) # note that we here use the original syndrome, not the forney syndrome
34                                                                                                                                   # (because we will correct both errors and erasures, so we need the full syndrome)
35     # check if the final message is fully repaired
36     synd = rs_calc_syndromes(msg_out, nsym)
37     if max(synd) > 0:
38         raise ReedSolomonError("Could not correct message")     # message could not be repaired
39     # return the successfully decoded message
40     return msg_out[:-nsym], msg_out[-nsym:] # also return the corrected ecc block so that the user can check()

Python注释:列表erase_pos和err_pos使用+运算符连接。

  这是完全的错误与丢失纠正Reed-Solomon解码器所需的最后一部分。如果您想深入研究勘误(错误和擦除纠正)解码器的内部工作原理,您可以阅读由Richard E. Blahut撰写的优秀著作《数据传输的代数代码【Algebraic Codes for Data Transmission】》(2003)。

  数学注释:在一些软件实现中,特别是那些使用为线性代数和矩阵运算优化的语言的人,你会发现算法(编码,Berlekamp-Massey等)看起来完全不同,并且使用傅里叶变换。这是因为这是完全等同的:当用谱估计术语表述时,解码Reed-Solomn由傅立叶变换(症状计算)组成,随后是频谱分析(Berlekamp-Massey或Euclidian算法),随后是傅立叶逆变换变换(Chien搜索)。有关更多信息,请参阅Blahut书籍[4]。事实上,如果您正在使用针对线性代数进行优化的编程语言,或者您想要使用快速线性代数库,那么使用傅立叶变换可能是一个非常好的主意,因为它现在非常快速(特别是快速傅里叶变换或数字理论变换【Number Theoretic Transform】[5])。

 

用一个例子来总结

  以下是如何使用刚刚完成的函数,以及如何解码错误和丢失的示例:

 1 # Configuration of the parameters and input message
 2 prim = 0x11d
 3 n = 20 # set the size you want, it must be > k, the remaining n-k symbols will be the ECC code (more is better)
 4 k = 11 # k = len(message)
 5 message = "hello world" # input message
 6 
 7 # Initializing the log/antilog tables
 8 init_tables(prim)
 9 
10 # Encoding the input message
11 mesecc = rs_encode_msg([ord(x) for x in message], n-k)
12 print("Original: %s" % mesecc)
13 
14 # Tampering 6 characters of the message (over 9 ecc symbols, so we are above the Singleton Bound)
15 mesecc[0] = 0
16 mesecc[1] = 2
17 mesecc[2] = 2
18 mesecc[3] = 2
19 mesecc[4] = 2
20 mesecc[5] = 2
21 print("Corrupted: %s" % mesecc)
22 
23 # Decoding/repairing the corrupted message, by providing the locations of a few erasures, we get below the Singleton Bound
24 # Remember that the Singleton Bound is: 2*e+v <= (n-k)
25 corrected_message, corrected_ecc = rs_correct_msg(mesecc, n-k, erase_pos=[0, 1, 2])
26 print("Repaired: %s" % (corrected_message+corrected_ecc))
27 print(‘‘.join([chr(x) for x in corrected_message]))

  应该输出以下内容:

Original:  [104, 101, 108, 108, 111,  32, 119, 111, 114, 108, 100, 145, 124, 96, 105, 94, 31, 179, 149, 163]
Corrupted: [  0,   2,   2,   2,   2,   2, 119, 111, 114, 108, 100, 145, 124, 96, 105, 94, 31, 179, 149, 163]
Repaired:  [104, 101, 108, 108, 111,  32, 119, 111, 114, 108, 100, 145, 124, 96, 105, 94, 31, 179, 149, 163]
hello world

 

结论与深入探讨

  里德-所罗门码的基本原理已在本文中提出。已经包含了针对特定实现(QR码使用通用里德-所罗门编解码器纠正误读)的Python代码。这里介绍的代码是非常通用的,可用于QR码以外的任何需要纠正错误/丢失的目的,如文件保护,网络连接等。这些想法的许多变化和改进都是可能的,因为编码理论是非常丰富的学习领域。

  如果你的代码仅用于你自己的数据(例如,你希望能够生成并阅读自己的QR码),那很好,但如果你打算使用他人提供的数据(例如,你想要读取和解码其他应用的QR码),那么这个解码器可能还不够,因为有一些隐藏的参数在这里被简化了(即:发生器/alpha数和第一个连续的根)。如果你想解码由其他库生成的里德-所罗门码,则需要使用通用里德-所罗门编解码器,这将允许指定自己的参数,甚至超出字段2 ^ 8。

  在补充资源页面上,你会发现这里介绍的代码的扩展通用版本,可以使用它来解码几乎任何Reed-Solomon代码,还有一个生成素数多项式列表的函数和一个检测未知RS码的参数。请注意,无论使用哪个参数,修复功能应始终保持一致:对数/反对数表和生成多项式的生成值不会更改Reed-Solomon代码的结构,因此始终可以获得相同的功能,无论参数如何。

  可能注意到的一个直接问题是,我们只能编码最多256个字符的消息。这种限制可以通过几种方式来规避,其中最常见的有三种:

  1. 使用更高的伽罗瓦域(例如216,其将允许65536个字符或232, 264, 2128等)。这里的问题在于,对于大的多项式来说,编码和解码里德-所罗门所需的多项式计算变得非常昂贵(大多数算法是在二次时间内,最有效的是n log n ,如数理论变换[5])。

  2. 通过“组块”,这意味着您只需通过256个字符块编码你的大数据流。

  3. 使用包含分组大小的变体算法,如Cauchy Reed-Solomon(见下文)。

  如果你想进一步研究,有很多关于里德-所罗门码的书籍和科学文章,一个很好的起点是着名的作者理查德 布拉胡特(Richard Blahut)。 此外,还有很多不同的方法可以对里德-所罗门码进行编码和解码,因此会发现许多不同的算法,特别是解码(Berlekamp-Massey,Berlekamp-Welch,欧几里德算法等)。

  如果你正在寻找更高的性能,在文献中找到这里介绍的算法的几种变体,例如Cauchy-Reed-Solomon。编程实现在你的Reed-Solomon编解码器的性能中也扮演着重要角色,这可能会导致1000倍的速度差异。更多相关信息,请阅读其他资源的“优化性能”部分

  即使近期最优的前向纠错算法由于其速度很快而如今已经风靡一时(如LDPC码,Turbo码等),Reed-Solomon是一种最佳的FEC,这意味着它可以达到已知的理论极限作为Singleton的界限。实际上,这意味着RS可以同时纠正多达2 * e + v <=(n-k)的错误和丢失,其中e是错误数,v是丢失个数,k是消息大小,n是消息+代码大小和(n-k)最小距离。这并不是说近乎最佳的FEC是无用的:它们比Reed-Solomon所能达到的速度难以想象,并且它们可能受到悬崖效应的影响较小(这意味着它们仍然可以部分解码消息的部分,即使存在纠正所有错误的错误太多了,与RS相反,这肯定会失败,甚至通过对错误信息进行解码而无声无息地进行解码[6]),但它们肯定无法纠正与Reed-Solomon一样多的错误。选择接近最优的FEC和最优的FEC主要是速度问题。

  最近,Reed-Solomon的研究领域自发现w:List_decoding (不要与软解码混淆)之后恢复了一些活力,它允许解码/修复比理论最优限制更多的符号。其核心思想是,而不是标准的Reed-Solomon只做一个独特的解码(意味着它总是导致一个单一的解决方案,如果它不能,因为它超过理论极限,解码器将返回一个错误或错误的结果),带有列表解码的Reed-Solomon仍然会尝试超出极限解码并获得几个可能的结果,但通过巧妙地检查不同的结果,通常可以区分只有一个可能是正确的多项式。

  一些列表解码算法已经可用,它允许修复高达n - sqrt(n*k)[7] 而不是 2*e+v <= (n-k),以及其他列表解码算法(更高效或解码更多符号)目前正在研究中。

 

写给编码人员的里德-所罗门编码

标签:class   include   基本   Python程序员   pre   交换   its   sage   就会   

原文地址:https://www.cnblogs.com/swordc007/p/9151205.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!