MSK144 — цифровой протокол, разработанный Джо Тейлором (K1JT) и его командой в 2016 году для проведения связей через метеорное рассеивание.
В предыдущей части были рассмотрены общие характеристики протокола, а также алгоритмы, формирующие сигнал. В этой части статьи рассматриваются механизмы детектирования и декодирования сигналов в протоколе MSK144.
Статья может быть интересна радиолюбителям и людям, интересующимся темой цифровой обработки сигналов.
Общее представление
Как было рассмотрено в предыдущей статье, формирование сигнала, отправляемого в радиоэфир, состоит из этапов упаковки бит данных, расчета контрольной суммы, вычислению LDPC-кодов, после чего, к данным добавляются маркеры (или якоря) синхронизации, состоящие из последовательности бит 01110010, прежде чем всё это преобразовать в дискретный сигнал звуковой частоты.
Обратная задача, при приеме сигнала, заключается в поиске маркера синхронизации с последующей фазовой и частотной коррекцией сигнала, после чего, выполняя IQ преобразования, выделить аналоговое представление бит данных, которые далее проходят процесс декодирования, сверки контрольных сумм и так далее, пока на выходе не останутся только данные сообщения.
Анализ сигнала
Анализ сигнала состоит из нескольких этапов (рисунок 3):
фильтрация сигнала;
детекция сигнала;
декодирования фрейма.
Фильтрация сигнала
Перед тем, как начать анализировать сигнал, осуществляется фильтрация сигнала полосовым фильтром, задача которого исключить незадействованные области частот анализируемого сигнала и оставить только ту полосу, в которой присутствуют интересующие сигналы с данными.
В MSK144 применяется полосовой фильтр на основе преобразований Фурье. Отличительными особенностями фильтров такого типа является то, что фильтрация осуществляется в частотной области, что позволяет добиться высокой степени аттенюации, конфигурировать АЧХ фильтра с любой степенью наклона и задание множества полос частот; а также, помимо аттенюации, осуществлять усиление сигнала (эквализация).
На рисунке 4 показана схема фильтра.
Принцип работы фильтра заключается в том, что исходный сигнал подвергается быстрому преобразованию Фурье, результатом которого является представление сигнала в комплексной форме, в виде амплитуд и фаз в разрезе частот; после умножения амплитудных значений на требуемую характеристику, осуществляется обратное преобразование Фурье, позволяющее обратно получить аналоговый сигнал.
В Python реализовать фильтр на основе преобразования Фурье можно используя функционал библиотеки numpy.
import numpy as np
import numpy.typing as npt
def fourier_bpf(signal: npt.NDArray[np.float64], n_fft: int,
response: npt.NDArray[np.float64]) -> npt.NDArray[np.complex128]:
freq_d = np.fft.fft(signal, n=n_fft)
response_len = len(response)
freq_d[:response_len] *= response
freq_d[0] = 0.5 * freq_d[0]
freq_d[response_len: n_fft] = complex(0, 0)
time_d = np.fft.ifft(freq_d, n=n_fft)
return time_dНа вход функции fourier_bpf подается исходный сигнал и АЧХ фильтра, задаваемая аргументом response. Аргумент n_fft определяет размер блока данных для преобразования Фурье; значение задается степенью числа 2.
Частоты сигналов MSK144 в основном располагаются в диапазоне от 1.0 КГц до 2.0 КГц, фильтр должен подавлять все остальные частоты, находящиеся за пределами диапазона (рисунок 6).
На рисунке 5 представлена АЧХ фильтра для фильтрации сигналов MSK144. Так как сигналы подвержены частотному смещению, полоса пропускания фильтра имеет несколько большую ширину захвата.
За формирование кривой АЧХ фильтра отвечает функция filter_response.
MSK144_CENTER_FREQ = 1500
def filter_response(n_fft: int, sample_rate: int) -> npt.NDArray[np.float64]:
t = 1 / 2000
beta = 0.1
df = sample_rate / n_fft
nh = (n_fft // 2) + 1
response = np.full(nh, 1, dtype=np.float64)
for i in range(nh):
f = abs(df * i - MSK144_CENTER_FREQ)
if (1 + beta) / (2 * t) >= f > (1 - beta) / (2 * t):
response[i] = 0.5 * (1 + np.cos((np.pi * t / beta) * (f - (1 - beta) / (2 * t))))
elif f > (1 + beta) / (2 * t):
response[i] = 0
return responseФункция принимает на вход размер выборки преобразования Фурье n_fft, и частоту дискретизации сигнала sample_rate. Результатом является массив значений от 0 до 1.
Кривая АЧХ описывается уравнением:
Точки, соответствующие аттенюации в 3дБ (что соответствует 0.5 от исходного значения амплитуды) соответствуют частотам 0.5 КГц и 2.5 КГц соответственно. Гладкость характеристики в зонах аттенюации определяется функцией косинуса а ширина значением константы beta, чем меньше значение, тем сильнее кривая будет стремиться к прямоугольной форме.
На рисунке 5 изображены спектры сигналов до и после применения фильтра с характеристикой filter_response.
Недостатками фильтров такого типа являются требования в вычислительным ресурсам, а также эффект Гиббса — то есть появление нежелательных осцилляций на границах сигнала, что требует применения к сигналу сглаживающих оконных функций.
Частотное смещение
Сигналы протокола MSK144 формируются в слышимом человеческим ухом диапазоне частот, с тем расчетом, что их прием и передача в радиоэфире осуществляется с помощью радиолюбительского трансивера, рассчитанного на передачу голоса; так как опорные генераторы в любительских трансиверах имеют некоторую погрешность в точности частоты, то отправляемые и принимаемые сигналы имеют частотный дрейф, в связи с этим, принятый сигнал необходимо подстраивать под частоту, на которой будет происходить декодирование.
Смещение частоты сигнала реализуется функцией shift_freq.
def shift_freq(complex_signal: npt.NDArray[np.complex128],
freq: float, sample_rate: int) -> npt.NDArray[np.complex128]:
phi = 2 * np.pi * freq / sample_rate
step = complex(np.cos(phi), np.sin(phi))
phasor = complex(1, 1)
signal = np.zeros(len(complex_signal), dtype=np.complex128)
for i, cs_val in enumerate(complex_signal):
phasor *= step
signal[i] = phasor * cs_val
return signalФункция shift_freq принимает на вход отфильтрованный преобразованиями Фурье сигнал в комплексном представлении, аргумент freq определяет на сколько герц необходимо сместить сигнал; sample_rate — частота дискретизации сигнала.
Функция реализует фазовое вращение сигнала в каждой точке на угол phi, соответствующему смещению на заданную частоту. Фазор реализует комплексную экспоненту.
В общем виде смещение может быть записано как: 
.
Детекция сигнала
Задача детекции сигнала заключается в коррекции фазовых и частотных ошибок сигнала и синхронизация по маркерам.
Функция детектора состоит из этапов:
поиск сигналов-кандидатов по амплитуде;
коррекция фазы и частотной ошибки;
синхронизация;
выделение фрейма.
Этап фазовой коррекции применяется для снижения межсимвольных интерференций сигнала, связанных с тем, что в большинстве трансиверов используется узкополосные фильтры, способные вносить фазовые изменения сигнала.
Маркер синхронизации
Для выделения маркера с помощью автокорреляционной функции, формируется эталонный, искомый сигнал.
MSK144_SYNC = [0, 1, 1, 1, 0, 0, 1, 0] # 0x72 in binary
SAMPLES_PER_WORD = 12
SAMPLES_PER_HALF_WORD = SAMPLES_PER_WORD // 2
WORD_SAMPLES = np.sin([i * np.pi / SAMPLES_PER_WORD for i in range(SAMPLES_PER_WORD)])
SYNC_WORDS = np.array([2 * b - 1 for b in MSK144_SYNC])
SYNC_I = np.array([
sample * SYNC_WORDS[j * 2 + 1]
for j in range(4)
for sample in WORD_SAMPLES
])
SYNC_Q = np.array([
pp * SYNC_WORDS[j * 2]
for j in range(4)
for pp in (WORD_SAMPLES[SAMPLES_PER_HALF_WORD:] if j == 0 else WORD_SAMPLES)
])
SYNC_WAVEFORM = np.array([complex(SYNC_I[i], SYNC_Q[i]) for i in range(42)])
SYNC_WAVEFORM_LEN = len(SYNC_WAVEFORM)
SYNC_WAVEFORM_CONJ = np.conj(SYNC_WAVEFORM)В MSK144 для маркера синхронизации используется последовательность бит 01110010 (0x72), задаваемая в MSK144_SYNC.
WORD_SAMPLES определяет дискретное представление участка сигнала для формирования символа MSK.
В SYNC_WORDS формируется последовательность сдвигов частоты, в соответствии с последовательностью бит в MSK144_SYNC.
Формирование сигнала с MSK манипуляцией осуществляется в комплексной форме, через IQ преобразование, четные биты записываются в SYNC_Q, нечетные в SYNC_I соответственно (рисунок 7).
Результирующий сигнал SYNC_WAVEFORM записывается в комплексном виде, где синфазная составляющая часть сигнала записывается в вещественную часть значений, а квадратурная в мнимую.
SYNC_WAVEFORM_CONJ хранит комплексно-сопряженное представление SYNC_WAVEFORM, используемое в вычислении автокорреляции.
Поиск кандидатов
При поиске кандидатов происходит анализ спектра сигнала в полосе захвата, оценка силы сигнала, его соотношения сигнал/шум и определение погрешности частоты.
Детектор сигналов реализуется функцией detect_signals, являющаяся методом класса MSK144Monitor.
MSK144_NSPM = 864
MSK144_NPTS = 3 * MSK144_NSPM
MSK144_NFFT = 864
MSK144_FREQ_CENTER = 1500
MSK144_FREQ_HI = MSK144_FREQ_CENTER + 500
MSK144_FREQ_LO = MSK144_FREQ_CENTER - 500
SMOOTH_WINDOW_LEN = SAMPLES_PER_WORD
SMOOTH_WINDOW = (1 - np.cos([i * np.pi / SMOOTH_WINDOW_LEN for i in range(SMOOTH_WINDOW_LEN)])) / 2
def detect_signals(self, signal: npt.NDArray[np.complex128], signal_len: int, start: float,
max_cand: int = 16, tolerance: int = 150) -> typing.Generator[LogItem, None, None]:
...Функция detect_signals принимает на вход отфильтрованный сигнал signal в комплексном виде, параметр signal_len определяет общую длину исходного сигнала (так как после фильтра на основе преобразований Фурье результирующие массивы по длине отличаются от исходного сигнала, параметр передается отдельно).
Параметр start определяет временную позицию в сигнале, относительно которой необходимо осуществлять анализ (например для ручного анализа участка сигнала).
Параметр tolerance позволяет задать ширину полосы частот анализируемого сигнала.
dt = 1 / self.sample_rate
df = self.sample_rate / MSK144_NFFT
f_hi = MSK144_FREQ_HI * 2
f_lo = MSK144_FREQ_LO * 2
f_hi_id_top = int((f_hi + 2 * tolerance) / df)
f_hi_id_bot = int((f_hi - 2 * tolerance) / df)
f_lo_id_top = int((f_lo + 2 * tolerance) / df)
f_lo_id_bot = int((f_lo - 2 * tolerance) / df)
f_2Khz = int(f_lo / df)
f_4Khz = int(f_hi / df)
...Переменные dt и df определяют минимальный временной шаг, относительно частоты дискретизации, и минимальный шаг изменения частоты в спектре после преобразования Фурье.
Значения f_hi и f_lo определяют минимальную и максимальную частоты MSK модуляции.
Группы переменных "f_hi_id_top, f_hi_id_bot" и "f_lo_id_top, f_lo_id_bot" определяют полосу частот, в рамках которых будет производиться анализ.
step_size = MSK144_NSPM // 4
steps_count = (signal_len - MSK144_NPTS) // step_size
det_amp = np.zeros(steps_count)
det_snr = np.zeros(steps_count)
det_f_err = np.full(steps_count, -999.99)
...Блок массивов "det_amp, det_snr, det_f_err" аккумулирует параметры амплитуд, SNR и отклонения частот на каждой итерации анализа участка сигнала.
steps_real = 0
for step in range(steps_count):
part_start = step_size * step
part_end = part_start + MSK144_NSPM
if part_end > signal_len:
break
part = signal[part_start: part_end]
part = part ** 2
part[:SMOOTH_WINDOW_LEN] *= SMOOTH_WINDOW
part[MSK144_NSPM - SMOOTH_WINDOW_LEN:MSK144_NSPM] *= SMOOTH_WINDOW[::-1]
spec = np.fft.fft(part, MSK144_NFFT)
amps = np.abs(spec) ** 2
...В цикле до steps_count из общего сигнала выделяется участок part размером MSK144_NSPM, который далее возводится в степень 2 (part ** 2), для выделения несущих тонов в спектре и повышения их амплитуд относительно шума (по этой причине ранее объявленные переменные со значениями частот умножались на 2).
Затем к границам сигнала применяется сглаживающая функция SMOOTH_WINDOW, чья характеристика соответствует оконной функции Ханна, изображенная на рисунке 8.
Быстрое преобразование Фурье позволяет получить спектр (spec) и АЧХ (amps) анализируемого участка сигнала.
# Hi freq
peak_hi = f_hi_id_bot + np.argmax(amps[f_hi_id_bot:f_hi_id_top])
delta_hi = -((spec[peak_hi - 1] - spec[peak_hi + 1]) / (
2 * spec[peak_hi] - spec[peak_hi - 1] - spec[peak_hi + 1])).real
amp_hi = amps[peak_hi]
amp_avg_hi = (np.sum(amps[f_hi_id_bot:f_hi_id_top]) - amp_hi) / (f_hi_id_top - f_hi_id_bot)
snr_hi = amp_hi / (amp_avg_hi + 0.01)
...На участке частот спектра в окрестностях частоты 4.0 КГц ищется максимальное значение peak_hi. Значение delta_hi определяет смещение частоты пика. Значение переменной amp_avg_hi определяет среднее значение амплитуд шума для расчета значения SNR в переменной snr_hi.
Аналогичным образом осуществляется поиск и оценка шумовой картины для частоты 2.0 КГц.
# Low freq
peak_lo = f_lo_id_bot + np.argmax(amps[f_lo_id_bot:f_lo_id_top])
delta_lo = -((spec[peak_lo - 1] - spec[peak_lo + 1]) / (
2 * spec[peak_lo] - spec[peak_lo - 1] - spec[peak_lo + 1])).real
amp_lo = amps[peak_lo]
amp_avg_lo = (np.sum(amps[f_lo_id_bot:f_lo_id_top]) - amp_lo) / (f_lo_id_top - f_lo_id_bot)
snr_lo = amp_lo / (amp_avg_lo + 0.01)
# Errors
f_err_hi = (peak_hi + delta_hi - f_4Khz) * df / 2
f_err_lo = (peak_lo + delta_lo - f_2Khz) * df / 2
f_err = f_err_hi if amp_hi >= amp_lo else f_err_lo
...Значения переменных f_err_hi и f_err_lo определяют разность между искомой и найденной частотами.
Тройка массивов "det_amp, det_snr, det_f_err" для каждого участка проанализированного сигнала содержит амплитуды сигналов, значение SNR и значение частотной погрешности.
Следующий этап заключается в отборе сигналов-кандидатов, на основе ранее результатов детектора частот.
indices = np.argsort(det_amp[:steps_real])
amp_median = det_amp[indices[steps_real // 4]]
if amp_median == 0:
amp_median = 1
det_amp[:steps_real] /= amp_median
...Амплитудные значения во временных окнах нормализуются относительно медианного значения.
time_arr = np.zeros(max_cand, dtype=np.float64)
freq_arr = np.zeros(max_cand, dtype=np.float64)
snr_arr = np.zeros(max_cand, dtype=np.float64)
count_cand = 0
for ip in range(max_cand):
il = np.argmax(det_amp[:steps_real])
if det_amp[il] < 3.5:
break
if abs(det_f_err[il]) <= tolerance:
time_arr[count_cand] = ((il - 0) * step_size + MSK144_NSPM / 2) * dt
freq_arr[count_cand] = det_f_err[il]
snr_arr[count_cand] = 12 * np.log10(det_amp[il]) / 2 - 9
count_cand += 1
det_amp[il] = 0
...В цикле из общего скопа задетектированных частот отбираются только те, которые не превышают заданную ошибку несущей частоты и имеют достаточное амплитудное значение, превышающее порог в 3.5.
Пороговое значение в 3.5 выбрано как минимально приемлемое значение для решения задач обработки акустических сигналов, так как меньшие значения считаются близкими к уровню шума и плохо пригодными для достоверного извлечения полезного сигнала.
if count_cand < 3:
for ip in range(max_cand - count_cand):
if ip >= max_cand - count_cand:
break
il = np.argmax(det_snr[:steps_real])
if det_snr[il] < 12.0:
break
if abs(det_f_err[il]) <= tolerance:
time_arr[count_cand] = ((il - 0) * step_size + MSK144_NSPM / 2) * dt
freq_arr[count_cand] = det_f_err[il]
snr_arr[count_cand] = 12 * np.log10(det_snr[il]) / 2 - 9
count_cand += 1
det_snr[il] = 0
...Если отобранных кандидатов по критерию амплитуды сигнала нашлось менее 3, дополнительно осуществляется поиск кандидатов по SNR. Такой механизм позволяет обнаруживать слабые сигналы при тропосферном прохождении.
if count_cand > 0:
indices = np.argsort(time_arr[:count_cand])
hashes = set()
for cand in range(count_cand):
with suppress(NextCand):
cand_idx = indices[cand]
frame_mid = int(time_arr[cand_idx] * self.sample_rate)
if frame_mid < MSK144_NPTS / 2:
frame_mid = MSK144_NPTS // 2
if frame_mid > signal_len - MSK144_NPTS / 2:
frame_mid = signal_len - MSK144_NPTS // 2
t0 = time_arr[cand_idx] + dt * (start)
part = signal[frame_mid - MSK144_NPTS // 2: frame_mid + MSK144_NPTS // 2]
f_err = freq_arr[cand_idx]
snr = 2 * int(snr_arr[cand_idx] / 2)
snr = max(-4.0, min(24.0, snr))
part = self.shift_freq(part, -(MSK144_FREQ_CENTER + f_err), self.sample_rate)
...Для каждого кандидата выбирается участок сигнала из найденных временных границ, производится коррекция ошибки несущей частоты функцией shift_freq.
Синхронизация
Следующий этап — синхронизация с найденными сигналами.
sync1 = np.zeros(MSK144_NPTS, dtype=np.complex128)
sync2 = np.zeros(MSK144_NPTS, dtype=np.complex128)
for i in range(MSK144_NPTS - (56 * 6 + 42)):
sync1[i] = sum(part[i: i + SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
sync2[i] = sum(part[i + 56 * 6: i + 56 * 6 + SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
sync_corr = abs(sync1) * abs(sync2)
...В sync1 и sync2 считаются значения АКФ сигнала и ранее рассмотренного маркера синхронизации SYNC_WAVEFORM_CONJ.
Значение sync_corr определяет общее значение АКФ для обоих маркеров синхронизации в сигнале.
peaks = []
for i in range(6):
peak = np.argmax(sync_corr)
sync_corr[max(0, peak - 7): min(MSK144_NPTS - 56 * 6 - 42, peak + 7)] = 0.0
peaks.append(peak)
...Формируется список шести пиков АКФ маркеров синхронизации для последующей синхронизации.
for peak in peaks:
peak_corr = np.zeros(6, dtype=np.complex128)
for i in range(6):
j = peak + i
if peak + 11 + MSK144_NSPM < MSK144_NPTS:
s1 = slice(j + 6, j + 6 + MSK144_NSPM, 6)
s2 = slice(j, j + MSK144_NSPM, 6)
else:
s1 = slice(j + 6, MSK144_NPTS, 6)
s2 = slice(j, MSK144_NPTS - 6, 6)
peak_corr[i] = np.sum((part[s1] * np.conj(part[s2])) ** 2)
peak_corr_idx = np.argmax(np.abs(peak_corr))
if peak_corr_idx <= 2:
peak_corr_idx -= 1
if peak_corr_idx > 2:
peak_corr_idx -= 7
...Для найденных шести кипов вычисляется корреляция для тонкой подстройки джиттера, чтобы попасть в центр глазковой диаграммы.
for id in range(3):
if id == 1:
sign = -1
elif id == 2:
sign = 1
else:
sign = 0
ic = peak + peak_corr_idx + sign
if ic < 0:
ic = ic + MSK144_NSPM
cca = np.sum(part[ic:ic + SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
if ic + 56 * 6 + 42 < MSK144_NPTS:
ccb = np.sum(part[ic + 56 * 6:ic + 56 * 6 + SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
cfac = ccb * np.conj(cca)
f_err_2 = np.atan2(cfac.imag, cfac.real) / (2 * np.pi * 56 * 6 * dt)
else:
ccb = np.sum(part[ic - 88 * 6:ic - 88 * 6 + SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
cfac = ccb * np.conj(cca)
f_err_2 = np.atan2(cfac.imag, cfac.real) / (2 * np.pi * 88 * 6 * dt)
freq_est = int(MSK144_FREQ_CENTER + f_err + f_err_2)
...В freq_est вычисляется итоговое значение ошибки отклонения несущей частоты.
Выделение фрейма
Финальный этап детекции — выделение фрейма с сигналом для передачи его в декодер.
for idf in range(5):
if idf == 0:
delta_f = 0
elif idf % 2 == 0:
delta_f = idf
else:
delta_f = -(idf + 1)
subpart = self.shift_freq(part, -(f_err_2 + delta_f), self.sample_rate)
subpart = np.roll(subpart, -(ic - MSK144_NSPM))
for avg_pattern in range(8):
if avg_pattern == 0:
frame = subpart[MSK144_NSPM: MSK144_NSPM + MSK144_NSPM]
elif avg_pattern == 1:
frame = subpart[MSK144_NSPM - 432: MSK144_NSPM - 432 + MSK144_NSPM]
frame = np.roll(frame, 432)
elif avg_pattern == 2:
frame = subpart[2 * MSK144_NSPM - 432: 2 * MSK144_NSPM - 432 + MSK144_NSPM]
frame = np.roll(frame, 432)
elif avg_pattern == 3:
frame = subpart[:MSK144_NSPM]
elif avg_pattern == 4:
frame = subpart[2 * MSK144_NSPM: 2 * MSK144_NSPM + MSK144_NSPM]
elif avg_pattern == 5:
frame = subpart[:MSK144_NSPM] + subpart[MSK144_NSPM:MSK144_NSPM + MSK144_NSPM]
elif avg_pattern == 6:
frame = subpart[MSK144_NSPM: MSK144_NSPM + MSK144_NSPM] + subpart[2 * MSK144_NSPM: 2 * MSK144_NSPM + MSK144_NSPM]
elif avg_pattern == 7:
frame = subpart[:MSK144_NSPM] + subpart[MSK144_NSPM:MSK144_NSPM + MSK144_NSPM] + subpart[2 * MSK144_NSPM: 2 * MSK144_NSPM + MSK144_NSPM]
...Перебор значений delta_f реализует частотный джиттеринг, используемый для тонкой коррекции частотной ошибки, которая могла быть недостаточно точно скомпенсирована на предыдущих этапах.
В subpart содержится скорректированная по частоте и выровненная по времени часть анализируемого сигнала.
Значение avg_pattern определяет схему выборки фрейма для декодирования.
Скорректированный по частоте и подогнанный под временные границы фрейм отправляется в функцию декодирования. Цикл джиттеринга и перебора схем прерывается в случае успешного декодирования фрейма и извлечения из него данных.
if x := self.decode_fame(frame, self.LDPC_ITERATIONS):
status, payload, eye_opening, bit_errors = x
if status.crc_extracted in hashes:
raise NextCand
hashes.add(status.crc_extracted)
df_hv = freq_est - MSK144_FREQ_CENTER
log_item = LogItem(
snr=snr,
dT=t0,
dF=df_hv,
p_avg=avg_pattern + 1,
freq=freq_est,
eye_open=eye_opening,
bit_err=bit_errors,
payload=payload,
crc=status.crc_extracted
)
yield log_item
raise NextCand
if len(hashes) >= 3:
breakРезультат декодирования заносится в LogItem, хранящий в себе информацию о принятом сигнале, его качестве, а также декодированное сообщение.
При декодировании трех и более сообщений в сигнале, цикл декодирования прерывается.
Описание датакласса LogItem.
from dataclasses import dataclass
@dataclass
class LogItem:
snr: float
dT: float
dF: float
p_avg: int
freq: int
eye_open: float
bit_err: int
payload: typing.ByteString
crc: intОписание класса-исключения для управления циклом обработки кандидатов.
class NextCand(Exception):
passДекодирование фрейма
Разлагая участки сигнала, в которых, предположительно, должны содержаться биты данных, на синфазную и квадратурную составляющие и сравнивая значения с эталонным сигналом, получая тем самым коэффициент правдоподобия.
Декодирование осуществляется функцией decode_fame.
Функция принимает на вход данные декодируемого фрейма, извлеченные в детекторе. Аргумент max_iterations определяет сколько раз декодер LDPC будет осуществлять попыток восстановить данные, прежде, чем они будут отбракованы.
def decode_fame(self, frame: npt.NDArray[np.complex128],
max_iterations: int) -> typing.Optional[typing.Tuple[DecodeStatus, typing.ByteString, float, int]]:
sync_corr_1 = sum(frame[:SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
sync_corr_2 = sum(frame[56 * 6: 56 * 6 + SYNC_WAVEFORM_LEN] * SYNC_WAVEFORM_CONJ)
sync_corr = sync_corr_1 + sync_corr_2
phase_0 = np.atan2(sync_corr.imag, sync_corr.real)
corr_factor = complex(np.cos(phase_0), np.sin(phase_0))
frame *= corr_factor.conjugate()
...Для фазовой коррекции сигнала, происходит вычисление корреляции по маркерным сигналам (sync_corr, sync_corr_1, sync_corr_2). Фазовая коррекция осуществляется умножением сигнала на комплексное сопряженное значение фазы corr_factor.
soft_bits = np.zeros(MSK144_BITS_COUNT, dtype=np.float64)
hard_bits = np.zeros(MSK144_BITS_COUNT, dtype=np.int64)
for i in range(6):
soft_bits[0] += frame[i].imag * WORD_SAMPLES[i + 6]
soft_bits[0] += frame[i + (MSK144_NSPM - 5 - 1)].imag * WORD_SAMPLES[i]
for i in range(12):
soft_bits[1] += frame[i].real * WORD_SAMPLES[i]
for i in range(1, 72):
soft_sum = 0
for j in range(12):
soft_sum += frame[i * 12 - 6 + j].imag * WORD_SAMPLES[j]
soft_bits[2 * i] = soft_sum
soft_sum = 0
for j in range(12):
soft_sum += frame[i * 12 + j].real * WORD_SAMPLES[j]
soft_bits[2 * i + 1] = soft_sum
if soft_bits[2 * i] >= 0:
hard_bits[2 * i] = 1
if soft_bits[2 * i + 1] >= 0:
hard_bits[2 * i + 1] = 1
if soft_bits[0] >= 0:
hard_bits[0] = 1
if soft_bits[1] >= 0:
hard_bits[1] = 1
...В soft_bits аккумулируются коэффициенты правдоподобия бит данных для синфазной и квадратурной составляющих сигнала (real — I, imag — Q).
bad_sync_1 = 0
bad_sync_2 = 0
for i, sw in enumerate(SYNC_WORDS):
bad_sync_1 += (2 * hard_bits[i] - 1) * sw
bad_sync_2 += ((2 * hard_bits[i + 57 - 1] - 1) * sw)
bad_sync_1 = (MSK144_SYNC_LEN - bad_sync_1) // 2
bad_sync_2 = (MSK144_SYNC_LEN - bad_sync_2) // 2
bad_sync = bad_sync_1 + bad_sync_2
if bad_sync > 4:
return
...Переменная hard_bits используется для оценки точности синхронизации bad_sync.
soft_bits_std = np.std(soft_bits)
soft_bits /= soft_bits_std
...Коэффициенты правдоподобия бит нормализуются через стандартное отклонение.
sigma = 0.60
soft_bits_128 = np.concat([soft_bits[8: 9 + 47], soft_bits[64: 65 + 80 - 1]])
log128 = 2 * soft_bits_128 / sigma ** 2
...В soft_bits_128 конкатенируются информационные биты, тем самым исключаются биты маркеров синхронизации.
Нормализованные логарифмические значения коэффициентов правдоподобия передаются в декодер LDPC. Декодирование LDPC алгоритмом Belief Propagation, возвращает последовательность двоичных бит.
ldpc_errors, plain128 = bp_decode(log128, max_iterations)
if ldpc_errors > self.MAX_LDPC_ERRORS:
return None
...В plain128 содержатся декодированные методом Belief Propagation двоичные значения бит.
Значение ldpc_errors определяет количество ошибок бит остались неразрешенными в процессе декодирования. Если значение превышает порог MAX_LDPC_ERRORS, декодируемый фрейм исключается.
Финальный этап — проверка контрольных сумм и формирование отчета о качестве сигнала.
if not mskx_check_crc(plain128):
return None
payload = self.pack_bits(plain128, MSKX_LDPC_K)
crc_extracted = mskx_extract_crc(payload)
tones = list(msk144_encode(payload))
msg_bits = np.zeros(MSK144_BITS_COUNT)
for i in range(MSK144_BITS_COUNT - 1):
j = i + 1
if tones[i] == 0:
msg_bits[j] = msg_bits[i] if j % 2 == 1 else (msg_bits[i] + 1) % 2
else:
msg_bits[j] = (msg_bits[i] + 1) % 2 if j % 2 == 1 else msg_bits[i]
eye_top = 1
eye_bot = -1
for i in range(MSK144_BITS_COUNT):
if msg_bits[i] == 1:
eye_top = min(soft_bits[i], eye_top)
else:
eye_bot = max(soft_bits[i], eye_bot)
eye_opening = eye_top - eye_bot
bit_errors = np.count_nonzero(hard_bits != msg_bits)
return DecodeStatus(ldpc_errors, crc_extracted), payload, eye_opening, bit_errorsВ payload содержатся данные сообщения в байтовом представлении.
crc_extracted — контрольная сумма сообщения.
Формирование tones и msg_bits используется для оценки BER (bit error rate).
На основе значений soft_bits производится оценка сигнала по глазковой диаграмме (рисунки 9, 10).
Проверка CRC-13
При извлечении двоичных бит данных сообщения, осуществляется проверка корректности контрольной суммы принятого пакета.
MSKX_CRC_WIDTH = 13
MSKX_LDPC_K = 90
def mskx_crc(msg_bits: npt.NDArray[np.uint8], msg_len: int) -> npt.NDArray[np.uint8]:
div = [1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1]
msg = np.zeros(msg_len + MSKX_CRC_WIDTH, dtype=np.uint8)
for i in range(msg_len + MSKX_CRC_WIDTH):
if i < 77:
msg[i] = msg_bits[i]
for i in range(msg_len):
if msg[i] != 0:
for j, d in enumerate(div):
msg[i + j] = msg[i + j] ^ d
return msg[msg_len:msg_len + MSKX_CRC_WIDTH]
def mskx_check_crc(msg_bits: npt.NDArray[np.uint8]) -> bool:
crc = mskx_crc(msg_bits, 83)
return np.array_equal(crc, msg_bits[MSKX_LDPC_K - MSKX_CRC_WIDTH:MSKX_LDPC_K])Предикат mskx_check_crc на вход принимает биты данных сообщения, вычисляет CRC-13 (mskx_crc), если биты вычисленной контрольной суммы совпадает с битами контрольной суммы пакета данных, тогда пакет считается корректным и может быть обработан дальше.
Вспомогательная функция извлечения контрольной суммы из пакета данных сообщения.
def mskx_extract_crc(msg: typing.ByteString) -> int:
return ((msg[9] & 0x07) << 10) | (msg[10] << 2) | (msg[11] >> 6)Монитор
Описанные выше функции объединены в один класс MSK144Monitor, задача которого предоставить интерфейс загрузки данных и инкапсулировать пайплайн и логику декодера.
class MSK144Monitor(AbstractMonitor):
MAX_LDPC_ERRORS = 18
LDPC_ITERATIONS = 10
@staticmethod
@cache
@jit(nopython=True)
def filter_response(n_fft: int, sample_rate: int) -> npt.NDArray[np.float64]:
...
@staticmethod
def fourier_bpf(signal: npt.NDArray[np.float64], n_fft: int,
response: npt.NDArray[np.float64]) -> npt.NDArray[np.complex128]:
...
@staticmethod
@jit(nopython=True)
def shift_freq(complex_signal: npt.NDArray[np.complex128],
freq: float, sample_rate: int) -> npt.NDArray[np.complex128]:
...
def __init__(self, sample_rate: int):
self.sample_rate = sample_rate
self.signal = np.zeros(0, dtype=np.float64)
def decode_fame(self, frame: npt.NDArray[np.complex128],
max_iterations: int) -> typing.Optional[typing.Tuple[DecodeStatus, typing.ByteString, float, int]]:
...
def detect_signals(self, signal: npt.NDArray[np.complex128], signal_len: int, start: float,
max_cand: int = 16, tolerance: int = 150) -> typing.Generator[LogItem, None, None]:
...
def decode(self, tm_slot_start: float) -> typing.Generator[LogItem, None, None]:
signal_len = min(len(self.signal), 30 * self.sample_rate)
signal_part = self.signal[:signal_len]
rms = np.sqrt(np.mean(signal_part ** 2))
if rms == 0 or np.isnan(rms):
rms = 1
part_norm = signal_part / (rms / 2)
n = int(np.log(signal_len) / np.log(2) + 1)
nfft = int(min(2 ** n, 1024 ** 2))
filter_response = self.filter_response(nfft, self.sample_rate)
part_filtered = self.fourier_bpf(part_norm, nfft, filter_response)
yield from self.detect_signals(part_filtered, signal_len, tm_slot_start)
def monitor_process(self, frame: npt.NDArray):
wave = frame * 0.000390625
avg = np.mean(wave)
wave -= avg
self.signal = waveMAX_LDPC_ERRORS определяет пороговое значение количества ошибок, свыше которого декодированный сигнал игнорируется.
LDPC_ITERATIONS задает максимальное количество итераций декодирования Belief Propagation алгоритма.
В __init__ задается частота семплирования сигнала sample_rate.
Метод monitor_process сохраняет дискретный сигнал для декодирования. Значение 0.000390625 используется для перехода от int16 к float64 значениям.
Генераторная функция decode осуществляет фильтрацию сигнала и передачу его в функцию-детектор. Генератор выдает данные, инкапсулированные в LogItem.
Полный пример декодера с использованием класса MSK144Monitor:
from scipy.io.wavfile import read
from decode_msk144 import MSK144Monitor
from message import message_decode
def main():
sample_rate, data = read("examples/signal.wav")
mon = MSK144Monitor(sample_rate=sample_rate)
mon.monitor_process(data)
for it in mon.decode(0):
msg = message_decode(it.payload)
msg_text = " ".join(it for it in msg if isinstance(it, str))
print(
f"dB: {it.snr:.3f}\t"
f"T: {it.dT:.3f}\t"
f"DF: {it.dF:.3f}\t"
f"Averaging pattern: {it.p_avg}\t"
f"Freq: {it.freq}\t"
f"Eye opening: {it.eye_open:.3f}\t"
f"Bit errors: {it.bit_err}\t"
f"CRC: {it.crc}\t"
f"Message text: {msg_text}"
)
if __name__ == '__main__':
main()Результат декодирования тестового сигнала:
dB: 2.000 T: 0.234 DF: -1.000 Averaging pattern: 1 Freq: 1499 Eye opening: 0.768 Bit errors: 1 CRC: 1509 Message text: CQ R9FEU LO87Заключение
В продолжении цикла статей, посвященных протоколу MSK144, в этой статье были рассмотрены принципы демодуляции и декодирования сигналов с минимальным частотным сдвигом. Были рассмотрены механизмы цифровой обработки сигналов, фазовая и частотная коррекции.
В статью не вошли описание алгоритмов декодирования LDPC, кодирования сообщений в текст, формирование эталонного сигнала и прочее, так как они были описаны в предыдущих статьях, посвященных FT8 и кодированию MSK144.
Декодер MSK144 может представлять интерес с точки зрения цифровой обработки сигналов.