Еще немного про LDPC коды

от автора

Предисловие

Всем привет! Я тут недавно начал разбираться в одной очень интересной теме, связанной с обработкой цифровой информации. Объектом моего исследования стали помехоустойчивые коды. Когда я был студентом, я даже писал студенческую научную статью, в которой представил код на Си для кодирования информации кодом Хэмминга на 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 (взято из этой статьи):

 Базовая матрица из стандарта IEEE 802.11n

Базовая матрица из стандарта IEEE 802.11n

Для того чтобы получить желаемую матрицу H из базовой матрицы, сначала определимся с битрейтом. Я решил взять для эксперимента самый низкий битрейт 1/2 (324 информационных бита из 648). Ячейки в базовой таблице говорят нам, что мы должны взять единичную матрицу 27×27 (размерность для кода 324 x 648) и циклически сдвинуть ее на n разрядов вправо, где n — значение ячейки. То есть, если у нас число 0, то сдвигаем на 0 и получаем просто единичную матрицу. Сдвинутая единичная матрица заменяет нашу ячейку в колонке из базовой матрицы. Прочерк обозначает нулевую матрицу. В итоге я получил такую матрицу H:

 Полученная матрица H

Полученная матрица H

P.s. эту матрицу я набирал в блокноте вручную чтобы быть максимально уверенным 🙂

А вот теперь, берем в наши шаловливые ручки линейную алгебру и метод RU и за работу!

Сначала нам предлагается разделить матрицу H на сектора A,B,C,D,E и T по такому принципу (взято из этой статьи):

 Представление матрицы H

Представление матрицы H

Зная что M = 324, M = 648, а M-G = 297, я разделил таблицу следующим образом:

 Моя разделенная матрица H

Моя разделенная матрица H

Далее, для оптимизации дальнейших вычислений, воспользуемся методом Гаусса:

Тут нет ничего сверхестественного, мы, по сути, просто прибавили ко второй строчке первую, умноженную на -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. график взят из этой статьи.

В упомянутой статье приводится два способа “подкрутки” результатов:

  1. Minsum normalized. Ввести коэффициент a, который можно изменять в диапазоне (0, 1]. Изменяя этот коэффициент, можно измерить BER (Bit Error Rate) при фиксированном ОСШ (отношение сигнал/шум) и обнаружить при каком значении a показатель BER будет наилучшим:

  2. Minsum offset. Можно ввести коэффициент смещения b, который также может изменяться в диапазоне (0, 1]. Данный коэффициент подбирается аналогично как в случае с Minsum normalized:

  3. Комбинированный Minsum. Никто не мешает нам использовать коэффициенты a и b одновременно:

    На этом, пожалуй все, спасибо за внимание! По мере написания диссертации буду стараться писать больше новых статей. Всем добра!


ссылка на оригинал статьи https://habr.com/ru/articles/830150/