发表更新大作业20 分钟读完 (大约2979个字)

极大似然技术编码译码

1 实验环境

语言：C++

2 实验内容

2.1 编程实现 $(m,n)$ 群码

2.1.1 输入输出

输入：$m,n$ 和一致性校验矩阵 $H$ 中的 $H_{m \times r}$ 部分，其中 $r=n-m$.

输出：所有的码字，按原码对应的十进制整数升序排列.

2.1.2 设计思路及核心算法

  1. 首先输入 $m,n$：

    1
    std::cin >> m >> n;

  2. 调用 inputH(m,n) 函数输入一致性校验矩阵 $H$：

    1
    BitMatrix H = inputH(m, n);

    inputH(m,n) 函数代码如下：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    uint64_t inputBit() {
        std::string line;
        std::cin >> line;
        uint64_t x = 0;
        int k = 0;
        for (size_t i = 0; i < line.size(); i++)
            if (line[i] == '1' || line[i] == '0')
                x |= (line[i] - '0') << k++;
        return x;
    }
    BitMatrix inputH(int m, int n) {
        BitMatrix H(n, n - m);
        for (int i = 0; i < m; i++) H.d[i] = inputBit();
        for (int i = m; i < n; i++) H.d[i] = 1 << (i - m);
        return H;
    }

    inputH 函数传入行列参数m,n，声明一个 BitMatrix 类，用来存放一致性校验矩阵的 $m \times r$部分.

  3. 调用 HtoH2(H) 函数生成 $m \times n$ 的编码矩阵：
    $$
    \left[
    \begin{matrix}
    I_m & H_{m\times r} \
    \end{matrix}
    \right]
    $$

    1
    2
    3
    4
    5
    6
    7
    8
    BitMatrix HtoH2(BitMatrix &H) {
        BitMatrix H2(H.n - H.m, H.n);
        for (int i = 0; i < H.n - H.m; i++) H2.d[i] = 1 << i;
        for (int i = 0; i < H.n - H.m; i++)
            for (int j = 0; j < H.m; j++)
                H2.d[i] |= (H.d[i] >> j & 1) << (j + H.n - H.m);
        return H2;
    }

  4. 定义一个 Bm 矩阵用于存放原码，矩阵每行从上到下按原码对应的十进制递增顺序排列，并进行运算：

    $$
    e_H(B^m)=B^m\times
    \left[
    \begin{matrix}
    I_m & H_{m\times r} \
    \end{matrix}
    \right]
    $$

    求出群码对应的矩阵并输出.

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    BitMatrix H2 = HtoH2(H), Bm(1 << m, m);
    for (int i = 0; i < (1 << m); i++)
        for (int j = 0; j < m; j++)
            Bm.d[i] |= (i >> j & 1) << (m - j - 1);
    BitMatrix N = Bm * H2;
    std::cout << "1. e(B^m):" << std::endl;
    for (int i = 0; i < (1 << m); i++) {
        std::cout << "  e(";
        printBit(Bm.d[i], m);
        std::cout << ") = ";
        printBit(N.d[i], n);
        std::cout << std::endl;
    }

2.2 构造陪集表

2.2.1 设计思路及核心算法

生成 coset table 的算法原理如下：

  1. 首先，定义了一个空的哈希表 cosets 用于存储 coset 的信息，以及一个空的哈希表 visited 用于标记已经访问过的向量；

  2. 创建一个向量 numbers，其中存储了所有可能的向量（共 2^n 个向量），并按照向量中 1 的个数进行排序，从少到多；

  3. 遍历 numbers 中的每个向量 i，对于每个向量，执行以下步骤：

    a. 创建一个空的 coset 对象 coset，用于存储属于同一个 coset 的向量；

    b. 对于矩阵 N 的每一行（表示每个 coset），将当前向量 i 与该行的向量 x 进行异或运算，得到新的向量 x' = N.d[j] ^ uint64_t(i)；

    c. 检查 visited 哈希表，如果在 visited 中找到了向量 x'，说明该向量已经属于某个 coset，跳过当前循环；

    d. 如果在 visited 中没有找到向量 x'，则将 x' 添加到当前 coset 中，并将 x' 添加到 visited 哈希表中，值为当前 coset 的索引；

    e. 将当前 coset 的陪集头 e 更新为 coset 中具有最少 1 的个数的向量；

    f. 如果当前 coset 的陪集头 e 在 cosets 哈希表中不存在，则将当前 coset 添加到 cosets 中，并将 e 对应的特征值（通过一致性校验矩阵 H 计算）添加到 syncdrome 哈希表中，值为当前 coset 的索引；

    g. 检查 cosetLeaders 的大小，如果已经达到了 $2^{n-m}$（即余类的数量），则停止循环.

  4. 最后，得到了所有的 coset 信息，并按照 coset 陪集头对应的十进制递增的顺序输出.

通过上述算法，可以生成所有的 coset，并找到每个 coset 的陪集头和特征值.

2.2.2 BitMatrix 类

该程序将译码时所需的比特矩阵进行了封装，能够实现乘法 (元素对应异或) 操作.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class BitMatrix {
public:
    int n, m;
    uint64_t d[100];
    BitMatrix(int n, int m): n(n), m(m) {
        memset(d, 0, sizeof(uint64_t) * n);
    }
    BitMatrix operator*(const BitMatrix& rhs) const {
        assert (m == rhs.n);
        BitMatrix res(n, rhs.m);
        for (int i = 0; i < n; i++)
            for (int k = 0; k < m; k++)
                if (d[i] >> k & 1)
                    res.d[i] ^= rhs.d[k];
        return res;
    }
    void print() const {
        for (int i = 0; i < n; i++) {
            printBit(d[i], m);
            std::cout << std::endl;
        }
    }
};

2.2.3 Coset 类

Coset 类由陪集、陪集头及添加操作组成，陪集中的元素使用一维向量进行存储，陪集头用整型变量单独存储，方便调用. 添加操作在每一次执行时更新一次陪集头.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class Coset {
public:
    std::vector<uint64_t> d; // coset
    uint64_t e; // coset leader
    Coset() : e(0) {}
    void add(uint64_t x) {
        if (d.empty())
            d.push_back(e = x);
        else {
            d.push_back(x);
            if (bitcount(x) < bitcount(e)) e = x;
        }
    }
};

2.3 极大似然译码

基于 2.2 构造陪集表 求出的左陪集表，用和函数 $e_H$ 相关联的极大似然译码函数 $d$ ，对 $x_t\in B^n$ 进行译码.

2.3.1 输入输出格式

  • 输入：首行是译码个数 $N$ ，整数格式；
  • 然后是 $N$ 行比特串 $x_t \in B^n$ ，字符串格式；
  • 输出：$N$行译码结果，格式为 $d(x_t) = b_1 b_2 … b_m$ . 如：d(001011)=001.

2.3.2 设计思路及核心算法

  1. 输入 $N$ ，然后输入待解码的 $N$ 个比特串. 输入后需要进行字符串转整数操作；

    1
    2
    3
    4
    5
    6
    7
    std::cin >> k;
    for (int i = 0; i < k; i++) {
        std::cerr << "Enter the word #" << i + 1 << ": ";
        uint64_t x = inputBit();
    /*
        ......................
    */

    inputBit() 函数如下：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    uint64_t inputBit() {
        std::string line;
        std::cin >> line;
        uint64_t x = 0;
        int k = 0;
        for (size_t i = 0; i < line.size(); i++)
            if (line[i] == '1' || line[i] == '0')
                    x |= (line[i] - '0') << k++;
        return x;
    }

  2. 确定输入的 $x_t$ 位于哪个 coset 中，可通过 visited[x] 得到所属陪集编号；

  3. 取出所属陪集的陪集头 $\epsilon$ ，与 $x_t$ 进行异或运算，得到的结果必然在群码对应的矩阵 $N$ 中找到，由此得到解码后结果.

    核心代码如下：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    std::vector<uint64_t > cosetLeaders;
        for (auto i : numbers) {
            Coset coset;
            bool notFound = true;
            for (int j = 0; j < N.n; j++) {
                uint64_t x = N.d[j] ^ uint64_t(i);
                if (visited.find(x) != visited.end()) {
                    notFound = false;
                    break;
                }
                visited[x] = cosetLeaders.size();
                coset.add(x);
            }
            if (notFound) {
                cosets[coset.e] = coset;
                syndrome[getSyndrome(coset.e, H)] = cosetLeaders.size();
                cosetLeaders.push_back(coset.e);
            }
            // if cosetLeaders.size() = 2^(n-m), then stop. 
            // Because there are only 2^(n-m) cosets.
            if (cosetLeaders.size() == (1 << (n - m))) break;
        }
    /*
    ............
    */
    printBit(Bm.d[NtoBm[cosetLeaders[visited[x]] ^ x]], m);

2.4 基于特征值的译码

基于 2.2 构造陪集表 求出的左陪集表，利用函数 $f_H:B^n \rightarrow B^r，f_H(x) = x\cdot H$ 求出每个陪集头对应的特征值 (Syndrome) ，并利用特征值，对$x_t \in B^n$ 进行译码.

2.4.1 输入输出格式

  • 输入：首行是译码个数N，整数格式；
  • 然后是N行比特串 $x_t \in B^n$ ，字符串格式；
  • 输出：N行译码所用的特征值、陪集头和译码结果，格式为 $f(x_t)$ = *** ，$\epsilon$ = ***，$d(xt) = b_1b_2…b_m$.

注：本程序采用的输出方式，对于每一个待解码二进制串输入，输出结果分为两行，第一行为采用极大似然译码计算的结果，结果前用序号 (a) 进行区分，第二行为基于特征值译码的结果，结果前用序号 (b) 区分，例如：

1
2
(a) d(10111) = 10
(b) f(10111) = 100, ε = 00100, d(10111)

2.4.2 设计思路及核心算法

  1. 输入 N ，然后输入待解码的 N 个比特串. 输入后需要进行字符串转整数操作；

  2. 通过函数 $f_H(x_t) = x_t \cdot H$ 计算 $x_t$ 特征值；

    1
    2
    3
    4
    5
    std::cout << " (b) f("; //特征值
    printBit(x, n);
    uint64_t s = getSyndrome(x, H);
    std::cout << ") = ";
    printBit(s, n - m);

    getSyndrome() 函数细节：

    1
    2
    3
    4
    5
    uint64_t getSyndrome(uint64_t x, BitMatrix &H) {
        BitMatrix X(1, H.n);
        X.d[0] = x;
        return (X * H).d[0];
    }

  3. 根据所得特征值在哈希表中对应到陪集头；

    1
    2
    std::cout << ", ε = "; //陪集头
    printBit(cosetLeaders[syndrome[s]], n);

  4. 陪集头与待译码比特串进行异或操作，得到的结果必然在群码对应的矩阵 $N$ 中找到，由此得到解码后结果.

    1
    2
    3
    4
    5
    std::cout << ", d("; //译码结果
    printBit(x, n);
    std::cout << ") = ";
    printBit(Bm.d[NtoBm[cosetLeaders[syndrome[s]] ^ x]], m);
    std::cout << std::endl;

3 代码优化

我们在四个地方优化了原本的代码，使得程序更加高效.

  1. 由于实验中满足 $1<m<n<7$，所以使用 uint64_t 替代数组，并使用位运算进行相应计算，同时降低了空间复杂度和时间复杂度；

    其中，矩阵乘法的核心代码如下：

    1
    2
    3
    4
    5
    // main.cpp, Line 30-33
    for (int i = 0; i < n; i++)
        for (int k = 0; k < m; k++)
            if (d[i] >> k & 1)
                res.d[i] ^= rhs.d[k];

    普通的矩阵乘法的时间复杂度是 $O(n^3)$，而使用位运算后，时间复杂度降低到了 $O(n^2)$.

  2. 当已经找到了所有的陪集时，就可以停止循环，不需要再继续遍历；

    1
    2
    // main.cpp, Line 132
    if (cosetLeaders.size() == (1 << (n - m))) break;

  3. 使用哈希表记录每一个码字所属的陪集，在解码的时候可以直接查找，而不需要再次遍历；

    1
    2
    // main.cpp, Line 123, 记录每一个码字所属的陪集
    visited[x] = cosetLeaders.size();
    1
    2
    // main.cpp, Line 155, 直接解码
    printBit(Bm.d[NtoBm[cosetLeaders[visited[x]] ^ x]], m);

    该哈希表还可以快速地判断当前陪集是否重复，减少了计算量；

    1
    2
    3
    4
    5
    // main.cpp, Line 119-122, 判断当前陪集是否重复
    if (visited.find(x) != visited.end()) {
        notFound = false;
        break;
    }

    考虑到 $1<m<n<7$，所以码字的范围最大不超过 $2^8 = 256$，所以实际上可以使用数组代替哈希表，降低部分常数.

  4. 特征值解码方式也使用了类似的方式；

    1
    2
    // main.cpp, Line 128, 记录每一个特征值对应的陪集(头)
    syndrome[getSyndrome(coset.e, H)] = cosetLeaders.size();
    1
    2
    // main.cpp, Line 155, 直接获取特征值对应的陪集头并进行异或运算
    printBit(Bm.d[NtoBm[cosetLeaders[syndrome[s]] ^ x]], m);

  5. 使用了 lowbit 快速计算最低位的 1 的位置，加快了 bitcount 的速度.

    1
    2
    // main.cpp, Line 17
    while (n) n -= n & -n, s++;

4 实验成果示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
$ ./main # 运行
2 5 # 输入m n
011
101 #这两行是一致性校验矩阵中的 H_{m*r} 部分
1. e(B^m):
  e(00) = 00000
  e(01) = 01101
  e(10) = 10011
  e(11) = 11110
2. Cosets:
  00000 | 00000 01101 10011 11110
  00001 | 00001 01100 10010 11111
  00010 | 00010 01111 10001 11100
  00100 | 00100 01001 10111 11010
  01000 | 01000 00101 11011 10110
  10000 | 10000 11101 00011 01110
  00110 | 00110 01011 10101 11000
  01010 | 01010 00111 11001 10100
Enter the number of words to decode: 3 # 在输出码字和陪集之后，在这里输入你要测试的次数
Enter the word #1: 10111
 (a) d(10111) = 10
 (b) f(10111) = 100, ε = 00100, d(10111) = 10 #输出译码的结果
······

附录：源代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
#include <iostream>
#include <string>
#include <vector>
#include <unordered_map>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cstdint>
#include <cassert>
void printBit(uint64_t x, int n) {
    for (int i = 0; i < n; i++)
        std::cout << (x >> i & 1);
}
template <typename T>
int bitcount(T n) {
    int s = 0;
    while (n) n -= n & -n, s++;
    return s;
}
class BitMatrix {
public:
    int n, m;
    uint64_t d[100];
    BitMatrix(int n, int m): n(n), m(m) {
        memset(d, 0, sizeof(uint64_t) * n);
    }
    BitMatrix operator*(const BitMatrix& rhs) const {
        assert (m == rhs.n);
        BitMatrix res(n, rhs.m);
        for (int i = 0; i < n; i++)
            for (int k = 0; k < m; k++)
                if (d[i] >> k & 1)
                    res.d[i] ^= rhs.d[k];
        return res;
    }
    void print() const {
        for (int i = 0; i < n; i++) {
            printBit(d[i], m);
            std::cout << std::endl;
        }
    }
};
class Coset {
public:
    std::vector<uint64_t> d; // coset
    uint64_t e; // coset leader
    Coset() : e(0) {}
    void add(uint64_t x) {
        if (d.empty())
            d.push_back(e = x);
        else {
            d.push_back(x);
            if (bitcount(x) < bitcount(e)) e = x;
        }
    }
};
uint64_t inputBit() {
    std::string line;
    std::cin >> line;
    uint64_t x = 0;
    int k = 0;
    for (size_t i = 0; i < line.size(); i++)
        if (line[i] == '1' || line[i] == '0')
            x |= (line[i] - '0') << k++;
    return x;
}
BitMatrix inputH(int m, int n) {
    BitMatrix H(n, n - m);
    for (int i = 0; i < m; i++) H.d[i] = inputBit();
    for (int i = m; i < n; i++) H.d[i] = 1 << (i - m);
    return H;
}
BitMatrix HtoH2(BitMatrix &H) {
    BitMatrix H2(H.n - H.m, H.n);
    for (int i = 0; i < H.n - H.m; i++) H2.d[i] = 1 << i;
    for (int i = 0; i < H.n - H.m; i++)
        for (int j = 0; j < H.m; j++)
            H2.d[i] |= (H.d[i] >> j & 1) << (j + H.n - H.m);
    return H2;
}
uint64_t getSyndrome(uint64_t x, BitMatrix &H) {
    BitMatrix X(1, H.n);
    X.d[0] = x;
    return (X * H).d[0];
}
int main() {
    int m, n;
    std::cin >> m >> n;
    assert (1 < m && m < n && n < 64);
    BitMatrix H = inputH(m, n);
    BitMatrix H2 = HtoH2(H), Bm(1 << m, m);
    for (int i = 0; i < (1 << m); i++)
        for (int j = 0; j < m; j++)
            Bm.d[i] |= (i >> j & 1) << (m - j - 1);
    BitMatrix N = Bm * H2;
    std::cout << "1. e(B^m):" << std::endl;
    for (int i = 0; i < (1 << m); i++) {
        std::cout << "  e(";
        printBit(Bm.d[i], m);
        std::cout << ") = ";
        printBit(N.d[i], n);
        std::cout << std::endl;
    }
    std::unordered_map<uint64_t, Coset> cosets;
    std::unordered_map<uint64_t, int> visited, NtoBm, syndrome;
    for (int i = 0; i < N.n; i++) NtoBm[N.d[i]] = i;
    std::vector<int> numbers;
    for (int i = 0; i < (1 << n); i++) numbers.push_back(i);
    std::sort(numbers.begin(), numbers.end(), [&](int a, int b) {
        int ca = bitcount(a), cb = bitcount(b);
        return ca == cb ? a > b : ca < cb;
    });
    std::vector<uint64_t > cosetLeaders;
    for (auto i : numbers) {
        Coset coset;
        bool notFound = true;
        for (int j = 0; j < N.n; j++) {
            uint64_t x = N.d[j] ^ uint64_t(i);
            if (visited.find(x) != visited.end()) {
                notFound = false;
                break;
            }
            visited[x] = cosetLeaders.size();
            coset.add(x);
        }
        if (notFound) {
            cosets[coset.e] = coset;
            syndrome[getSyndrome(coset.e, H)] = cosetLeaders.size();
            cosetLeaders.push_back(coset.e);
        }
        // if cosetLeaders.size() = 2^(n-m), then stop.
        // Because there are only 2^(n-m) cosets.
        if (cosetLeaders.size() == (1 << (n - m))) break;
    }
    std::cout << "2. Cosets:" << std::endl;
    for (auto i: cosetLeaders) {
        std::cout << "  ";
        printBit(i, n);
        Coset coset = cosets[i];
        std::cout << " |";
        for (size_t i = 0; i < coset.d.size(); i++) {
            std::cout << " ";
            printBit(coset.d[i], n);
        }
        std::cout << std::endl;
    }
    std::cerr << "Enter the number of words to decode: ";
    int k;
    std::cin >> k;
    for (int i = 0; i < k; i++) {
        std::cerr << "Enter the word #" << i + 1 << ": ";
        uint64_t x = inputBit();
        std::cout << " (a) d(";
        printBit(x, n);
        std::cout << ") = ";
        printBit(Bm.d[NtoBm[cosetLeaders[visited[x]] ^ x]], m);
        std::cout << std::endl;
        std::cout << " (b) f(";
        printBit(x, n);
        uint64_t s = getSyndrome(x, H);
        std::cout << ") = ";
        printBit(s, n - m);
        std::cout << ", ε = ";
        printBit(cosetLeaders[syndrome[s]], n);
        std::cout << ", d(";
        printBit(x, n);
        std::cout << ") = ";
        printBit(Bm.d[NtoBm[cosetLeaders[syndrome[s]] ^ x]], m);
        std::cout << std::endl;
    }
    return 0;
}

极大似然技术编码译码

http://m1yan.github.io/2023/12/20/极大似然技术编码译码/

作者

MiYan

发布于

2023-12-20

更新于

2023-12-20

许可协议

#

喜欢这篇文章？打赏一下作者吧

×