Предисловие
Всем привет! Я тут недавно начал разбираться в одной очень интересной теме, связанной с обработкой цифровой информации. Объектом моего исследования стали помехоустойчивые коды. Когда я был студентом, я даже писал студенческую научную статью, в которой представил код на Си для кодирования информации кодом Хэмминга на Arduino. Только вот коды Хемминга вряд ли можно применить в каких-нибудь сложных каналах связи по типу WiFi или LTE, поэтому я начал изучать другие коды. Немного погуглив, я понял что мейнстримом среди помехоустойчивых кодов являются LDPC коды.
На Хабре была статья на тему помехоустойчивого кодирования и LDPC кодов. В ней автор очень круто описал основные принципы обработки информации, закодированной LDPC кодом, и даже привел пример декодирования методом SPA и некоторые мысли о том, как это дело можно оптимизировать. Я решил привнести свою лепту и подготовил свою небольшую статью в которой расскажу про кодирование информации на примере метода Ричардсона-Урбанке (Richardson — Urbanke method), а также рассмотрю вариант декодирования информации методом minsum и различные способы оптимизации этого метода.
Зачем?
Дело в том, что я пишу диссертацию, связанную с разработкой цифровых систем связи для роботехнических систем. На моей кафедре мне посоветовали углубиться в тему SDR приемников. Мне сказали что SDR — это сейчас последний писк моды в мире роботехнических систем). На самом деле, технология и правда очень интересная и перспективная да и сама возможность создавать собственные приемники путем написания кода, а не пайки радиоэлементов у меня вызывала восторг) Теперь пайка транзисторов и других радиоэлементов кажется архаизмом (особенно учитывая то что катушку с добротностью, необходимой для системы связи на частотах 900 МГц я никогда не намотаю).
Я решил написать код на 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). Ячейки в базовой таблице говорят нам, что мы должны взять единичную матрицу 27×27 (размерность для кода 324 x 648) и циклически сдвинуть ее на n разрядов вправо, где n — значение ячейки. То есть, если у нас число 0, то сдвигаем на 0 и получаем просто единичную матрицу. Сдвинутая единичная матрица заменяет нашу ячейку в колонке из базовой матрицы. Прочерк обозначает нулевую матрицу. В итоге я получил такую матрицу H:
P.s. эту матрицу я набирал в блокноте вручную чтобы быть максимально уверенным 🙂
А вот теперь, берем в наши шаловливые ручки линейную алгебру и метод RU и за работу!
Сначала нам предлагается разделить матрицу H на сектора A,B,C,D,E и T по такому принципу (взято из этой статьи):
Зная что M = 324, M = 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 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 одновременно:
На этом, пожалуй все, спасибо за внимание! По мере написания диссертации буду стараться писать больше новых статей. Всем добра!
ссылка на оригинал статьи https://habr.com/ru/articles/830150/
Добавить комментарий