线性分组码之实现RM编码

参考资料：
http://www.mcs.csueastbay.edu/~malek/Class/Reed-Muller.pdf
Journal of Computer Science 4 (10): 792-798, 2008

RM码全称Reed-Muller Codes，是一种非常神奇的编码，早在1972年就应用在了水手9号火星探测器，采用RM(1,5)对火星的黑白照片进行编码处理。

RM码也是线性分组码，但是与通常的线性分组码描述方法(n,k,d)不同，采用另外一种参数描述方法：(r,m)，二者对应关系为
$$(n,k,d)=(2^m,\sum _{i=0}^r{C}_m^i,2^{m-r})$$

可见RM码的最小距离给的还是很大的，提供给我们纠错多位的可能。之前所说的RM(1,5)对应的线性分组码就是(32,6,16) woc能纠正7位发现1位

生成矩阵

下面我们看一下如何去设计一这样神奇的编码。同所有线性分组码一样，它也有它的生成矩阵
首先当$r=0$时，生成矩阵被定义为
$$G(0,m)=\begin{array}{c}\underbrace{\left[\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right]}\\ 2^m\end{array}\text{\hspace{1em}(1)}$$

当$r\ge 1$时，我们通过迭代的方式得到生成矩阵
$$G(m,m)=\left[\begin{array}{ccc}& G(m-1,m)& \\ 0& 0\cdots & 1\end{array}\right]\text{\hspace{1em}(2)}$$

$$G(r,m+1)=\left[\begin{array}{cc}G(r,m)& G(r,m)\\ 0& G(r-1,m)\end{array}\right]\text{\hspace{1em}(3)}$$

---

Exam:G(1,3)
$$G(1,3)=\left[\begin{array}{cc}G(1,2)& G(1,2)\\ 0& G(0,2)\end{array}\right]$$
$$G(1,2)=\left[\begin{array}{cc}G(1,1)& G(1,1)\\ 0& G(0,1)\end{array}\right]$$

综上求得G(1,3)
$$\begin{gathered} \left[\begin{array}{cccc}G(0,1)& G(0,1)& G(0,1)& G(0,1)\\ 01& 01& 01& 01\\ 0& & G(0,2)& \end{array}\right] \\ =\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 0& 1& 0& 1& 0& 1& 0& 1\\ 0& 0& 1& 1& 0& 0& 1& 1\\ 0& 0& 0& 0& 1& 1& 1& 1\end{array}\right] \end{gathered}$$

---

此外还有另一种构造方式
零阶RM码生成矩阵仍参照(1)
定义RM(1,m)的生成矩阵为：设矩阵H的列序号与列向量相同，那么RM(1,m)的生成矩阵为
$$G=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 0& & & \\ ⋮& & H& \\ 0& & & \end{array}\right]$$

一般来说很多人看到这个定义都会比较懵逼，一方面是数值与向量相等的操作比较费解，另一方面是这个相等是真的不好解释，下面给一个实例说明一下，仍以RM(1,3)的生成矩阵为例
$$\left[\begin{array}{cccccccc}& a_1& a_2& a_3& a_4& a_5& a_6& a_7\\ 1& 1& 1& 1& 1& 1& 1& 1\\ 0& 1& 0& 1& 0& 1& 0& 1\\ 0& 0& 1& 1& 0& 0& 1& 1\\ 0& 0& 0& 0& 1& 1& 1& 1\end{array}\right]$$

序号不要求从1起，也可以从0起

接下来是如何从一阶RM码得到高阶RM码。将组成H的行向量依次记为
$$H=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_m\end{array}\right]$$

高阶RM码生成矩阵就是在低阶生成矩阵的基础上添加全部同阶的行向量。这个地方还是不容易说明，依旧举例说明：

定义点乘运算：设有同维行向量$a(a_1,a_2,\cdots a_n)$和$b(b_1,b_2,\cdots b_n)$，二者点乘的结果为
$$ab=(a_1{b}_1,a_2{b}_2,\cdots ,a_n{b}_n)$$

由RM(1,3)得到RM(2,3)为例，方便起见将RM(1,3)写为如下形式:将首列的0替换为行向量的标识符
$$\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ x_1& 1& 0& 1& 0& 1& 0& 1\\ x_2& 0& 1& 1& 0& 0& 1& 1\\ x_3& 0& 0& 0& 1& 1& 1& 1\end{array}\right]$$

$x_1,x_2,x_3$最高阶次均为1次，因为是1阶组合；而$x_1{x}_2$为最高阶次为2次，这里的乘法运算就是前文提及的点乘，所以RM(2,3)的生成矩阵为
$$\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ x_1& 1& 0& 1& 0& 1& 0& 1\\ x_2& 0& 1& 1& 0& 0& 1& 1\\ x_3& 0& 0& 0& 1& 1& 1& 1\\ x_1{x}_2& 0& 0& 1& 0& 0& 0& 1\\ x_1{x}_3& 0& 0& 0& 0& 1& 0& 1\\ x_2{x}_3& 0& 0& 0& 0& 0& 1& 1\end{array}\right]$$

由于m个元素r阶组合的个数为$C_m^r$，所以可知RM(r,m)的信息位为${\sum }_{i=0}^r{C}_m^i$

译码

一阶译码

首先介绍一下什么是Kronecker Product：

矩阵A与B的Kronecker Product定义如下
$$A⨂B=\left[\begin{array}{cccc}{a}_{11}B& a_{12}B& \cdots & a_{1n}B\\ a_{21}B& a_{12}B& \cdots & a_{2n}B\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}B& a_{m2}B& \cdots & a_{mn}B\end{array}\right]$$

之后是定义n阶hadamard 矩阵

n阶hadamard 矩阵为nxn方阵，仅由 1 和 -1 俩个元素构成，且具有如下性质
$$H_n{H}_n^T=n{I}_n$$

hadamard矩阵是可以构造的，但这不是本文的重点 而且我也不会

之后是如何利用hadamard矩阵对RM码进行译码
首先利用2阶hadamard矩阵和kronecker Product定义如下矩阵的生成
$$H_m^j=I_{2^{m-j}}\otimes H_2\otimes I_{2^{j-1}}$$

已$H_3^3$为例
$$H_3^3=H\otimes I_4=\left[\begin{array}{cccccccc}1& 0& 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 1& 0& 0\\ 0& 0& 1& 0& 0& 0& 1& 0\\ 0& 0& 0& 1& 0& 0& 0& 1\\ 1& 0& 0& 0& -1& 0& 0& 0\\ 0& 1& 0& 0& 0& -1& 0& 0\\ 0& 0& 1& 0& 0& 0& -1& 0\\ 0& 0& 0& 1& 0& 0& 0& -1\end{array}\right]$$

现假设采用RM(1,m)编码，收到码字$r(r_0,r_1,\cdots r_{2^m-1})$，按如下步骤进行译码

1. 将接收码字中0替换为-1
   $$w_0=2w-[111\cdots 1]$$

2. 计算如下伴随式
   $$\begin{array}{rl}{w}_1& =w_0{H}_m^1\\ w_2& =w_1{H}_m^2\\ w_3& =w_2{H}_m^3\\ ⋮& \\ w_m& =w_{m-1}{H}_m^m\end{array}$$

3. 在行向量$w_m(w_{m,1}{w}_{m,2}\cdots w_{m,2^m-1})$中找位置$i$，使得$\forall pos\ne i$，都有$|w_{m,i}|>|w_{m,pos}|$

4. 根据$i$查询下表

   pos i	code v(i)
   0	000…0
   1	100…0
   2	010…0
   3	110…0
   …	…
   2m-1	111…1

   如果$w_{m,i}<0$，那么认为译码结果为$[0v(i)]$；如果$w_{m,i}>0$，则认为译码结果为$[1v(i)]$

5. 将求得的译码结果与生成矩阵$G(1,m)$相乘得到码字$w_{snd}$与接收码字$w$比较，核对结果

解码示例程序

完整一阶编译码程序-> 德莉莎世界第一可爱
下面给出matlab演示程序作为示例：

hadamard(n): 生成n阶hadamard矩阵
kron(A,B):计算矩阵A B的kronecker Product

close all
clear all
clc

H=hadamard(2);
w=[1 0 1 0 1 0 1 1];
%replace
w0=zeros(1,8);
for i=1:8
    w0(i)=2*w(i)-1;
end
%calculate w1 w2 w3
H13=kron(eye(4),H);
w1=w0*H13
H23=kron(kron(eye(2),H),eye(2));
w2=w1*H23;
H33=kron(H,eye(4));
w3=w2*H33
%get max element and its position
pos=1;
for i=2:8
    if w3(i)>w3(pos)
        pos=i;
    end
end
%decode
code=zeros(1,4);
if w3(pos)>0
    code(1)=1;
else
    code(1)=0;
end
pos=pos-1;
for i=2:4
    code(i)=mod(pos,2);
    pos=fix(pos/2);
end
%display
code