参考资料：
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{matrix} \underbrace{ \begin{bmatrix} 1&1&\cdots&1 \end{bmatrix} } \\ 2^m \end{matrix} \tag{1}$

当 $r\ge 1$时，我们通过迭代的方式得到生成矩阵
$G(m,m)= \begin{bmatrix} &G(m-1,m) \\ 0&0\cdots&1 \end{bmatrix} \tag{2}$

$G(r,m+1)= \begin{bmatrix} G(r,m) & G(r,m) \\ 0 & G(r-1,m) \end{bmatrix} \tag{3}$

Exam:G(1,3)
$G(1,3)= \begin{bmatrix} G(1,2) & G(1,2) \\ 0 & G(0,2) \end{bmatrix}$
$G(1,2)= \begin{bmatrix} G(1,1) & G(1,1) \\ 0 & G(0,1) \end{bmatrix}$

综上求得G(1,3)
$\begin{bmatrix} G(0,1) & G(0,1) & G(0,1) & G(0,1) \\ 01 & 01 & 01 & 01 \\ 0 && G(0,2) \end{bmatrix} \\ \quad \\ =\begin{bmatrix} 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{bmatrix}$

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

一般来说很多人看到这个定义都会比较懵逼，一方面是数值与向量相等的操作比较费解，另一方面是这个相等是真的不好解释，下面给一个实例说明一下，仍以RM(1,3)的生成矩阵为例
$\begin{bmatrix} & 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{bmatrix}$

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

接下来是如何从一阶RM码得到高阶RM码。将组成H的行向量依次记为
$H= \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_m \\ \end{bmatrix}$

高阶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替换为行向量的标识符
$\begin{bmatrix} 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{bmatrix}$

$x_1,x_2,x_3$最高阶次均为1次，因为是1阶组合；而 $x_1x_2$为最高阶次为2次，这里的乘法运算就是前文提及的点乘，所以RM(2,3)的生成矩阵为
$\begin{bmatrix} 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_1x_2 &0 & 0 & 1 & 0 & 0 & 0 & 1 \\ x_1x_3 &0 & 0 & 0 & 0 & 1 & 0 & 1 \\ x_2x_3 &0 & 0 & 0 & 0 & 0 & 1 & 1 \\ \end{bmatrix}$

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

译码

一阶译码

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

矩阵A与B的Kronecker Product定义如下
$A \bigotimes B = \begin{bmatrix} a_{11} B & a_{12} B &\cdots & a_{1n} B \\ a_{21} B & a_{12} B &\cdots & a_{2n} B \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} B & a_{m2} B & \cdots & a_{mn} B \end{bmatrix}$

之后是定义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 = \begin{bmatrix} 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{bmatrix}$

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

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

2. 计算如下伴随式
   $\begin{aligned} w_1 &= w_0 H^1_m \\ w_2 &= w_1 H^2_m \\ w_3 &= w_2 H^3_m \\ \vdots \\ w_m &= w_{m-1} H^m_m \end{aligned}$

3. 在行向量 $w_m (w_{m,1} \quad 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$，那么认为译码结果为 $[0\quad v(i)]$；如果 $w_{m,i}>0$，则认为译码结果为 $[1 \quad v(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