-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdecoder.py
212 lines (171 loc) · 6.79 KB
/
decoder.py
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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
"""Decoding module."""
import numpy as np
import warnings
import utils
from numba import njit, int64, types, float64
def decode(H, y, snr, maxiter=1000):
"""Decode a Gaussian noise corrupted n bits message using BP algorithm.
Decoding is performed in parallel if multiple codewords are passed in y.
Parameters
----------
H: array (n_equations, n_code). Decoding matrix H.这个是正确的
y: array (n_code, n_messages) or (n_code,). Received message(s) in the
codeword space.为什么这里的矩阵维度是这样的呢？因为当初初始的信息就是个奇怪的形状
maxiter: int. Maximum number of iterations of the BP algorithm.
Returns
-------
x: array (n_code,) or (n_code, n_messages) the solutions in the
codeword space.这里面还包含着多余的信息啊。应该是得到(k_code, n_messages)才对，所以要有个切片
"""
m, n = H.shape
bits_hist, bits_values, nodes_hist, nodes_values = utils._bitsandnodes(H)
#bits_hist行重，nodes_hist列重
_n_bits = np.unique(H.sum(0))
_n_nodes = np.unique(H.sum(1))
if (_n_bits * _n_nodes == 1).all():#矩阵做比较
print("矩阵比较结果")
solver = _logbp_numba_regular #这里是啥意思啊？怎么操作译码器的呢？
bits_values = bits_values.reshape(n, -1)#这里的retgular应该是LDPC矩阵的regular和非regular
nodes_values = nodes_values.reshape(m, -1)
else:
solver = _logbp_numba#内置的选择函数运行的代码
var = 10 ** (-snr / 10)
if y.ndim == 1:
y = y[:, None]
# step 0: initialization
Lc = 2 * y / var#这里就是计算对数似然比了。
_, n_messages = y.shape
Lq = np.zeros(shape=(m, n, n_messages))
Lr = np.zeros(shape=(m, n, n_messages))
for n_iter in range(maxiter):
Lq, Lr, L_posteriori = solver(bits_hist, bits_values, nodes_hist,
nodes_values, Lc, Lq, Lr, n_iter)
x = np.array(L_posteriori <= 0).astype(int)
product = utils.incode(H, x)#如果hx==0就收敛了。提前结束
if product:
break
if n_iter == maxiter - 1:
warnings.warn("""Decoding stopped before convergence. You may want
to increase maxiter""")
return x.squeeze()
output_type_log2 = types.Tuple((float64[:, :, :], float64[:, :, :],
float64[:, :]))
@njit(output_type_log2(int64[:], int64[:], int64[:], int64[:], float64[:, :],
float64[:, :, :], float64[:, :, :], int64), cache=True)
def _logbp_numba(bits_hist, bits_values, nodes_hist, nodes_values, Lc, Lq, Lr,
n_iter):
"""Perform inner ext LogBP solver."""
m, n, n_messages = Lr.shape
# step 1 : Horizontal
bits_counter = 0
nodes_counter = 0
for i in range(m):
# ni = bits[i]
ff = bits_hist[i]
ni = bits_values[bits_counter: bits_counter + ff]
bits_counter += ff
for j in ni:
nij = ni[:]
X = np.ones(n_messages)
if n_iter == 0:
for kk in range(len(nij)):
if nij[kk] != j:
X *= np.tanh(0.5 * Lc[nij[kk]])
else:
for kk in range(len(nij)):
if nij[kk] != j:
X *= np.tanh(0.5 * Lq[i, nij[kk]])
num = 1 + X
denom = 1 - X
for ll in range(n_messages):
if num[ll] == 0:
Lr[i, j, ll] = -1
elif denom[ll] == 0:
Lr[i, j, ll] = 1
else:
Lr[i, j, ll] = np.log(num[ll] / denom[ll])
# step 2 : Vertical
for j in range(n):
# mj = nodes[j]
ff = nodes_hist[j]
mj = nodes_values[nodes_counter: nodes_counter + ff]
nodes_counter += ff
for i in mj:
mji = mj[:]
Lq[i, j] = Lc[j]
for kk in range(len(mji)):
if mji[kk] != i:
Lq[i, j] += Lr[mji[kk], j]
# LLR a posteriori:
L_posteriori = np.zeros((n, n_messages))
nodes_counter = 0
for j in range(n):
ff = nodes_hist[j]
mj = nodes_values[nodes_counter: nodes_counter + ff]
nodes_counter += ff
L_posteriori[j] = Lc[j] + Lr[mj, j].sum(axis=0)
return Lq, Lr, L_posteriori
@njit(output_type_log2(int64[:], int64[:, :], int64[:], int64[:, :],
float64[:, :], float64[:, :, :], float64[:, :, :],
int64), cache=True)
def _logbp_numba_regular(bits_hist, bits_values, nodes_hist, nodes_values, Lc,
Lq, Lr, n_iter):
"""Perform inner ext LogBP solver."""
m, n, n_messages = Lr.shape
# step 1 : Horizontal
for i in range(m):
ni = bits_values[i]
for j in ni:
nij = ni[:]
X = np.ones(n_messages)
if n_iter == 0:
for kk in range(len(nij)):
if nij[kk] != j:
X *= np.tanh(0.5 * Lc[nij[kk]])
else:
for kk in range(len(nij)):
if nij[kk] != j:
X *= np.tanh(0.5 * Lq[i, nij[kk]])
num = 1 + X
denom = 1 - X
for ll in range(n_messages):
if num[ll] == 0:
Lr[i, j, ll] = -1
elif denom[ll] == 0:
Lr[i, j, ll] = 1
else:
Lr[i, j, ll] = np.log(num[ll] / denom[ll])
# step 2 : Vertical
for j in range(n):
mj = nodes_values[j]
for i in mj:
mji = mj[:]
Lq[i, j] = Lc[j]
for kk in range(len(mji)):
if mji[kk] != i:
Lq[i, j] += Lr[mji[kk], j]
# LLR a posteriori:
L_posteriori = np.zeros((n, n_messages))
for j in range(n):
mj = nodes_values[j]
L_posteriori[j] = Lc[j] + Lr[mj, j].sum(axis=0)
return Lq, Lr, L_posteriori
def get_message(tG, x):#这个get_message是怎么回事，为什么译码的结果还要这个东西来得到，还有就是并行加速是怎么回事？
"""Compute the original `n_bits` message from a `n_code` codeword `x`.
Parameters
----------
tG: array (n_code, n_bits) coding matrix tG.
x: array (n_code,) decoded codeword of length `n_code`.
Returns
-------
message: array (n_bits,). Original binary message.
"""
n, k = tG.shape
rtG, rx = utils.gausselimination(tG, x)
message = np.zeros(k).astype(int)
message[k - 1] = rx[k - 1]
for i in reversed(range(k - 1)):
message[i] = rx[i]
message[i] -= utils.binaryproduct(rtG[i, list(range(i+1, k))],
message[list(range(i+1, k))])
return abs(message)