Предисловие
Всем привет! Я тут недавно начал разбираться в одной очень интересной теме, связанной с обработкой цифровой информации. Объектом моего исследования стали помехоустойчивые коды. Когда я был студентом, я даже писал студенческую научную статью, в которой представил код на Си для кодирования информации кодом Хэмминга на Arduino. Только вот коды Хемминга вряд ли можно применить в каких‑нибудь сложных каналах связи по типу WiFi или LTE, поэтому я начал изучать другие коды. Немного погуглив, я понял что мейнстримом среди помехоустойчивых кодов являются LDPC коды.
На Хабре была статья на тему помехоустойчивого кодирования и LDPC кодов. В ней автор очень круто описал основные принципы обработки информации, закодированной LDPC кодом, и даже привел пример декодирования методом SPA и некоторые мысли о том, как это дело можно оптимизировать. Я решил привнести свою лепту и подготовил свою небольшую статью в которой расскажу про кодирование информации на примере метода Ричардсона‑Урбанке (Richardson — Urbanke method), а также рассмотрю вариант декодирования информации методом minsum и различные способы оптимизации этого метода.
Зачем?
Дело в том, что я пишу диссертацию, связанную с разработкой цифровых систем связи для роботехнических систем. На моей кафедре мне посоветовали углубиться в тему SDR приемников. Мне сказали что SDR — это сейчас последний писк моды в мире роботехнических систем). На самом деле, технология и правда очень интересная и перспективная да и сама возможность создавать собственные приемники путем написания кода, а не пайки радиоэлементов у меня вызывала восторг) Теперь пайка схемок и транзисторов и других радиоэлементов кажется архаизмом.
Я решил написать код на Python»е для кодирования и раскодирования информации для того чтобы использовать его как блок в среде GnuRadio. Для тех кто не знает: GnuRadio — это мощный программный пакет для обработки сигналов, получаемых с SDR приемников. Он позволяет преобразовывать сигналы (усиливать, переносить спектры, демодулировать) при помощи готовых блоков, а также писать свои на Python и C++. Программный комплекс GnuRadio может быть запущен на одноплатном ПК RuspberryPi, что позволяет внедрять разработанное решения в различные системы, требующие беспроводной связи.
Кодируем информацию
Давайте сначала рассмотрим классический метод кодирования. Предположим что у нас есть некоторая матрица H:
Для того чтобы закодировать информацию, нам необходимо получить из матрицы H порождающую матрицу G:
Кодирование информации сводится к матричному перемножению порождающей матрицы G и вектора с информационными битами x:
Чтобы проверить правильность принятых бит нам нужно перемножить матрицу H с вектором закодированного сообщения։
Если при перемножении получаем нулевой вектор - информация получена верно!
Однако, существуют различные алгоритмы, ускоряющие процесс кодирования информации. Ниже будет рассмотрен метод Ричардсона-Урбанке (В этой научной статье вы можете ознакомиться с ним подробнее). Этот метод предполагает замену перемножения матриц G и вектора s на некоторые нехитрые операции, которые мы рассмотрим далее.
Еще одно достоинство данного метода заключается в том, что его удобно использовать со стандартами, такими как 802.11 (наш любимый WiFi) где в самом стандарте нам предоставляется способ вычисления матрицы H, которую удобно применять в методе RU. Например, так нам предлагает вычислить матрицу H стандарт 802.11 (взято из этой статьи):
Для того чтобы получить желаемую матрицу H из базовой матрицы, сначала определимся с битрейтом. Я решил взять для эксперимента самый низкий битрейт 1/2 (324 информационных бита из 648). Ячейки в базовой таблице говорят нам, что мы должны взять единичную матрицу 27x27 (размерность для кода 324 x 648) и циклически сдвинуть ее на n разрядов вправо, где n - значение ячейки. То есть, если у нас число 0, то сдвигаем на 0 и получаем просто единичную матрицу. Сдвинутая единичная матрица заменяет нашу ячейку в колонке из базовой матрицы. Прочерк обозначает нулевую матрицу. В итоге я получил такую матрицу H:
P.s. эту матрицу я набирал в блокноте вручную чтобы быть максимально уверенным :)
А вот теперь, берем в наши шаловливые ручки линейную алгебру и метод RU и за работу!
Сначала нам предлагается разделить матрицу H на сектора A,B,C,D,E и T по такому принципу (взято из этой статьи):
Зная что M = 324, N = 648, а M-G = 297, я разделил таблицу следующим образом:
Далее, для оптимизации дальнейших вычислений, воспользуемся методом Гаусса:
Тут нет ничего сверхестественного, мы, по сути, просто прибавили ко второй строчке первую, умноженную на -ET-1.
Далее, для упрощения, выполним небольшую подстановку:
Мы уже близки к финалу! Предположим что наше сообщение является вектором x, где x = [s, p1, p2], где s — наша исходная (не закодированная информация), а p1 и p2 parity check биты, которые нам нужно вычислить. Нам известно что H@x = 0. Тогда для нахождения p1 и p2 нам нужно просто решить небольшое уравнение:
Я решил опустить некоторые математические выкладки и сразу вывел решение. Кстати, существует еще одно требование: D’ не должно быть сингулярным или, говоря по-русски, D’ должно быть таким, чтобы определитель матрицы D’ не был равен нулю. Если же D’ сингулярна, то нужно выполнить еще одно преобразование строк матрицы H’ так чтобы выполнялось det(D’) ≠ 0.
Теперь попробуем написать код на Python:
Код на Python
import numpy as np
from numpy.linalg import inv
H = np.loadtxt("./H.txt")
T = H[0:297, 351:648]
A = H[0:297, 0:324]
B = H[0:297, 324:351]
C = H[297:324, 0:324]
D = H[297:324, 324:351]
E = H[297:324, 351:648]
ET = ((-1*E)@inv(T))%2
ETA = (ET@A)%2
ETB = (ET@B)%2
ETT = (ET@T)%2
C_1 = (ETA + C)%2
D_1 = (ETB + D) % 2
E_1 = (ETT + E) % 2
#Тестируем
s = np.random.randint(0, 2, size = 324)#генерируем рандомное сообщение
s_T = np.transpose(s)
Cs_T = (C_1@s_T)%2
As_T = (A@s_T)%2
p1 = (inv(D_1)@Cs_T)%2
Bp_T = (B@np.transpose(p1))%2
p2 = ((-inv(T)%2)@((As_T + Bp_T)%2))%2
x = np.concatenate([s, p1, p2], axis=0)
#Проверка
print("Проверяем произведение H на x:\n", (H@np.transpose(x)) % 2)
print("Полученное сообщение:\n", x)
np.savetxt('text.txt', H);
np.savetxt('msg.txt', x);В ответ получаем следующее:
Проверяем произведение H на x:
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Полученное сообщение:
[1. 1. 0. 1. 0. 0. 0. 1. 0. 0. 1. 0. 0. 0. 1. 1. 0. 1. 1. 1. 1. 0. 1. 1.
0. 0. 0. 0. 1. 1. 0. 0. 0. 1. 1. 0. 1. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 0.
0. 1. 1. 0. 1. 0. 0. 1. 0. 1. 1. 1. 1. 1. 0. 1. 0. 1. 0. 1. 1. 1. 1. 0.
1. 1. 1. 1. 0. 1. 0. 1. 1. 0. 1. 1. 0. 1. 1. 1. 0. 1. 0. 0. 1. 0. 1. 1.
0. 0. 1. 0. 0. 1. 1. 1. 0. 1. 0. 0. 0. 0. 0. 1. 1. 0. 1. 0. 0. 1. 0. 0.
0. 1. 1. 1. 0. 0. 1. 0. 0. 0. 1. 1. 1. 0. 0. 1. 0. 1. 0. 1. 0. 1. 1. 0.
1. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 0. 0. 1. 1. 1. 0. 1. 0. 1. 0. 0. 0.
1. 1. 0. 1. 0. 1. 0. 1. 1. 1. 0. 1. 1. 1. 0. 1. 1. 1. 0. 0. 1. 1. 1. 0.
0. 0. 1. 1. 1. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0. 1. 0. 1. 1. 0. 1. 1. 0. 0.
1. 1. 0. 1. 0. 1. 1. 1. 0. 1. 1. 0. 1. 0. 0. 1. 1. 1. 0. 1. 1. 1. 1. 1.
1. 0. 1. 1. 1. 0. 0. 0. 0. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 1. 1. 0.
1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 0. 1. 0. 0. 1. 1. 1. 0. 1. 1. 1. 0. 0. 0.
1. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1. 0. 1. 1. 1. 0. 1. 0.
1. 1. 0. 1. 0. 0. 1. 1. 0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0. 1. 0. 0. 1. 0.
1. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 1. 0. 1. 0. 0. 1. 1. 0. 0. 1. 0. 0. 0.
1. 1. 0. 0. 0. 1. 0. 1. 1. 1. 1. 1. 0. 1. 0. 0. 0. 1. 0. 1. 1. 1. 1. 1.
1. 1. 0. 0. 1. 1. 1. 1. 1. 0. 1. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 1.
1. 1. 0. 1. 0. 0. 0. 1. 1. 1. 0. 1. 0. 1. 0. 1. 1. 1. 0. 0. 1. 0. 1. 1.
1. 0. 0. 1. 1. 1. 0. 0. 1. 0. 0. 1. 1. 0. 1. 0. 1. 0. 0. 1. 0. 1. 0. 1.
1. 0. 0. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 1. 0. 1. 0. 0. 1. 0. 1. 1. 1. 1.
1. 0. 1. 1. 1. 0. 1. 1. 1. 0. 1. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 1. 0. 0.
1. 1. 0. 1. 1. 1. 1. 0. 1. 0. 0. 0. 0. 1. 1. 0. 1. 1. 1. 0. 1. 1. 1. 1.
0. 0. 1. 0. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 1. 1. 1. 0. 0. 1. 1. 0. 0. 1.
0. 1. 0. 1. 1. 0. 1. 0. 1. 1. 0. 1. 0. 0. 0. 0. 1. 1. 0. 1. 1. 1. 0. 0.
0. 1. 0. 0. 0. 0. 0. 1. 1. 1. 1. 0. 1. 1. 0. 0. 1. 0. 1. 0. 1. 1. 1. 0.
0. 1. 1. 0. 1. 0. 1. 1. 0. 0. 0. 1. 1. 0. 1. 0. 0. 0. 0. 1. 1. 0. 0. 0.
1. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]Вуаля! Все работает так как надо!
Декодирование информации
Так как автор статьи рассмотрел метод SPA, я решил попробовать на основе результатов автора воплотить метод minsum. Если вы не знакомы с методом SPA, прочитайте пожалуйста эту статью чтобы понять что к чему. Данный метод по сути призван ускорить метод SPA. Он отличается тем, что операция нахождения тангенса при вычислении матрицы E:
заменяется на более простую операцию с перемножением элементов и поиском минимума:
У меня получился следующий алгоритм:
Алгоритм minsum на Python
import numpy as np
class CustomMinSum:
def __init__(self, H):
self.H = H
def nrz(self, arr):
for i in range(np.shape(arr)[0]):
if arr[i] <= -1:
arr[i] = 1
elif arr[i] > -1:
arr[i] = 0
elif arr[i] >= 1:
arr[i] = 0
else:
arr[i] = 1
return arr
def encode(self, r):
M = np.zeros(np.shape(self.H))
E = np.zeros(np.shape(M))
self.r = r
while True:
for j in range(np.shape(H)[0]):
M[j, :] = r*H[j, :]
E = self.evaluate_E(M, E)
l = self.evaluate_l(E)
l = self.nrz(l)
s = (H@np.transpose(l)) % 2
if np.prod(s == np.zeros(np.size(s))) == 1:
return l
else:
M = self.evaluate_M(E, M)
def evaluate_M(self, E, M):
for j in range(np.shape(self.H)[0]):
for i in range(np.shape(self.H)[1]):
if self.H[j,i] != 0:
M[j,i] = np.sum(E[:, i]) - E[j,i] + self.r[i]
M = M * self.H
return M
def evaluate_E(self, M, E):
for i in range(np.shape(M)[0]):
non_zero_mask = H[i,:] != 0
for j in range(np.shape(M)[1]):
if H[i, j] != 0:
E[i, j] = (np.prod(np.sign(M[i,:])[non_zero_mask]) / np.sign(M[i,j]))
exclude = np.arange(len(M[i,:])) != j
min = np.min(np.abs(M[i,:])[exclude & non_zero_mask])
E[i, j] = E[i, j]*min
return E
def evaluate_l(self, E):
return self.r + np.sum(E, axis=0)
H = np.loadtxt('H.txt')
r = np.loadtxt('msg.txt')
r = r*(-1)
res = CustomMinSum(H).encode(r)
print("Результат:", res)
#Проверяем что декодированные и изначально закодированные символы совпадают:
true_message = np.loadtxt('true_message.txt')
true_message[true_message == -1.34] = 0
true_message[true_message == 1.34] = 1
print("Совпадение: ", np.prod(res == true_message) == 1)Тестирование декодирования
В файле msg.txt у меня хранится результат кодирования алгоритмом RU из прошлого раздела. Я заменил в нем все единицы на 1.34, а нули на -1.34 чтобы эмитировать значения в виде логарифмированных коэффициентов правдоподобия LLR (предполагаем что наш демодулятор возвращает мягкие значения модуляции вместо единиц и нулей). Для эксперимента я случайным образом подменил 10 бит чтобы эмулировать ошибки, а в файл true_message сохранил сообщение без ошибок чтобы сравнить их и удостовериться что все OK. Получился такой результат:
Результат: [1. 0. 0. 0. 1. 0. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 0. 1. 0.
0. 1. 0. 0. 0. 1. 1. 0. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1.
0. 0. 0. 1. 1. 1. 0. 1. 1. 1. 1. 1. 0. 0. 1. 0. 0. 1. 0. 0. 0. 1. 0. 1.
0. 1. 1. 0. 1. 1. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 0. 1.
0. 0. 1. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 1. 0. 1. 0. 0. 1. 0. 0. 0. 0.
1. 0. 0. 0. 0. 0. 1. 0. 1. 0. 0. 1. 1. 1. 0. 0. 0. 1. 1. 0. 1. 1. 1. 0.
0. 1. 0. 1. 1. 1. 0. 0. 1. 1. 1. 0. 1. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1.
0. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0. 1.
1. 1. 1. 0. 1. 1. 1. 1. 0. 0. 0. 0. 1. 0. 1. 1. 0. 0. 1. 1. 1. 0. 0. 1.
0. 0. 1. 1. 0. 1. 0. 1. 1. 1. 1. 0. 1. 1. 0. 0. 1. 0. 1. 0. 1. 0. 1. 0.
0. 0. 1. 1. 0. 1. 0. 1. 1. 1. 1. 0. 1. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1.
1. 0. 1. 1. 0. 1. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 0. 1. 0. 1. 1. 1.
1. 1. 1. 1. 1. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 0. 0. 0.
1. 1. 1. 0. 1. 1. 0. 1. 0. 1. 0. 0. 1. 0. 0. 0. 1. 0. 1. 0. 0. 0. 1. 1.
0. 0. 1. 0. 1. 1. 1. 1. 1. 0. 0. 1. 0. 0. 1. 1. 0. 0. 0. 1. 0. 0. 0. 0.
1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 0. 0. 1. 0. 1. 1. 0. 1. 0. 1. 1. 1. 0. 0.
0. 0. 1. 1. 1. 1. 1. 0. 1. 0. 1. 1. 0. 1. 0. 1. 0. 0. 1. 0. 1. 0. 0. 1.
0. 1. 0. 0. 0. 1. 1. 1. 1. 1. 0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1.
1. 0. 0. 0. 0. 0. 0. 1. 0. 1. 1. 0. 1. 0. 0. 1. 0. 1. 0. 0. 0. 1. 0. 0.
1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 1. 0. 1. 0. 0. 1. 1. 1. 0. 0. 1. 1. 1. 0.
1. 0. 0. 1. 0. 0. 0. 1. 1. 1. 0. 1. 0. 0. 1. 1. 0. 0. 1. 0. 1. 0. 0. 1.
1. 0. 0. 1. 0. 1. 0. 0. 0. 0. 1. 0. 1. 0. 1. 0. 0. 1. 1. 0. 1. 1. 0. 1.
1. 1. 1. 1. 1. 1. 1. 0. 0. 1. 0. 1. 1. 1. 0. 1. 1. 0. 1. 0. 1. 1. 1. 0.
1. 0. 0. 1. 0. 1. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1.
1. 1. 0. 0. 1. 0. 1. 0. 0. 0. 0. 1. 0. 0. 1. 1. 1. 0. 1. 1. 0. 0. 1. 1.
0. 0. 1. 0. 0. 0. 1. 0. 1. 0. 1. 1. 0. 1. 0. 0. 0. 0. 1. 0. 1. 0. 0. 1.
1. 1. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 1. 0.]
Совпадение: TrueКак видно, все биты совпали!
Если верить этой статье, то алгоритм minsum проигрывает spa по энергетическому выигрышу кодирования на 0.2 — 0.5 дБ. Но это поправимо, ведь мы можем «подкрутить» наши вычисления так, чтобы они были ближе к тому, как если бы мы вычисляли E через tanh.
Если взглянуть на график BER, то можно увидеть, что проигрыш алгоритма MinSum перед SPA действительно составляет примерно 0.2 — 0.3 дБ.
График BER от SNR
P.s. график взят из этой статьи.
В упомянутой абзацем ранее статье приводится два способа «подкрутки» результатов:
Minsum normalized. Ввести коэффициент a, который можно изменять в диапазоне (0, 1]. Изменяя этот коэффициент, можно измерить BER (Bit Error Rate) при фиксированном ОСШ (отношение сигнал/шум) и обнаружить при каком значении a показатель BER будет наилучшим:
Minsum offset. Можно ввести коэффициент смещения b, который также может изменяться в диапазоне (0, 1]. Данный коэффициент подбирается аналогично как в случае с Minsum normalized:
Комбинированный Minsum. Никто не мешает нам использовать коэффициенты a и b одновременно:
На этом, пожалуй все, спасибо за внимание! По мере написания диссертации буду стараться писать больше новых статей. Всем добра!
