Предисловие

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

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:

заменяется на более простую операцию с перемножением элементов и поиском минимума:

У меня получился следующий алгоритм:

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)

Результат: [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 одновременно:

   

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