Q65 — цифровой протокол, разработанный Джо Тейлором (K1JT) и его командой в 2021 году для проведения минимальных связей в условиях сложных трасс прохождения радиосигнала.
В предыдущей части части были рассмотрены общая структура протокола и алгоритмы формирования сигнала. В этой части статьи рассматриваются принципы цифровой обработки сигналов и синхронизации в протоколе Q65.
Статья может быть интересна радиолюбителям и людям, интересующимся темой цифровой обработки сигналов.
Общее представление
Ранее было рассмотрено, что процесс формирования сигнала состоит из этапов вычисления контрольных сумм, формирования QRA (Q-ary Repeat Accumulation), формирования FSK-сигнала с одним выделенным под синхронизацию тоном.
При получении сигнала, прежде чем решить обратную задачу по демодуляции и декодированию данных (рисунок 1), сигнал предварительно подвергается таким обработкам, как: компенсация эффекта доплера, подавление помех и борьба с быстрыми затуханиями (fast fading), возникающими при многолучевом распространении сигнала.
Анализ сигнала
Анализ сигнала производится в несколько этапов (рисунок 2):
выделение символов в спектре;
устранение эффекта Доплера;
нормализация спектра;
синхронизация;
извлечение информационных символов;
устранение локальных выбросов;
применение модели быстрых затуханий.
Выделение символов в спектре
Так как в протоколе Q65 используется FSK манипуляция, выделение символов осуществляется в частотной и временной областях сигнала; исходный сигнал разделяется на временные интервалы, к которым применяется быстрое преобразование Фурье (БПФ) для получения амплитудных значений в частотной области.
import typing
import numpy as np
import numpy.typing as npt
def symbol_spectra(self, num_bins: int, num_times: int) -> npt.NDArray[np.float64]:
fac = (1 / np.iinfo(np.int16).max) * 0.01
sym_spec = np.zeros((num_times, num_bins), dtype=np.float64)
for i in range(0, num_times, 2):
f_beg = i * self.sym_steps
f_end = f_beg + self.sym_samps
spec = np.fft.fft(self.signal[f_beg:f_end] * fac, n=self.fft_size)[:num_bins]
sym_spec[i, :] = np.abs(spec) ** 2
if self.smooth > 1:
for _ in range(self.smooth):
sym_spec[i] = smooth_121(sym_spec[i])
if i >= 2:
sym_spec[i - 1] = 0.5 * (sym_spec[i - 2] + sym_spec[i])
return sym_specФункция symbol_spectra, являющаяся частью класса Q65Monitor (более подробно будет описан далее), для дискретного сигнала self.signal осуществляет быстрое преобразование Фурье. Аргумент num_bins определяют на сколько частотных бинов будет раскладываться анализируемый сигнал; параметр num_times определяет из скольки временных отрезков сигнал будет разделен.
Квадрат модуля от результатов FFT позволяет получить амплитудные значения частот.
Примечание: анализ с шагом 2 и интерполяция промежуточных значений спектров используется для оптимизации.
В качестве результата функция symbol_spectra возвращает матрицу, содержащую спектры символов протокола в разрезе времени.
На рисунке 3 приведен спектр анализируемого сигнала, полученный как результат функции symbol_spectra. По оси ординат расположены частотные бины (нумерация относительная), по оси абсцисс — временные отсчеты.
Эффект Доплера
Для снижения шума и флуктуаций спектральной плотности, возникающим вследствие эффекта Доплера, в Q65 применяется сглаживание спектра.
Функция smooth_121 реализует сглаживание с использованием скользящего среднего, весами ¼, ½, ¼.
def smooth_121(a: np.ndarray):
v = np.array([0.25, 0.5, 0.25])
convolved = np.convolve(a, v, mode='same')
return np.concat([a[:1], convolved[1:-1], a[-1:]])Степень сглаживания (self.smooth) определяется шириной сигнала; для сигналов с широкой полосой применяется сглаживание спектра с большей степенью размытия.
В python размытие осуществляется функцией свертки convolve из библиотеки numpy, значения весов задаются массивом v.
Значение self.smooth задается как степень числа 2.
self.smooth = max(1, int(0.5 * self.q65_type ** 2))Так, для Q65B количество сглаживаний составляет 2, для Q65C — 4 и т.д..
Нормализация спектра
Для повышения чувствительности и снижения ложных срабатываний при детекции сигналов, в Q65 применяется нормализация спектра по 45 перцентилю. Применение нормализации с таким значением перцентиля позволяет обеспечить подавление шумов и выбросов, при этом обеспечивается устойчивость к импульсным помехам без сильного искажения спектра.
sym_spec = self.symbol_spectra(iz, jz)
for j in range(jz):
t_s = sym_spec[j, self.i0 - 64:self.i0 - 64 + LL]
if (base := np.percentile(t_s, 45)) == 0:
base = 0.000001
sym_spec[j, :] /= baseЗначение переменной jz (более детально о ней речь пойдет далее) определяет спектры символов в разрезе времени, значение t_s частотную область спектра символов FSK, в пределах которой будет рассчитываться перцентиль.
Значение перцентиля для символа определяется с помощью функции percentile из библиотеки numpy.
Выражение sym_spec[j, :] /= base осуществляет нормализацию.
Значения 45 в качестве перцентиля выведено эмпирическим путем, так как порядка 55% амплитуд приходится на шум, в то время как полезный сигнал использует оставшиеся 45% спектра.
Использование больших значений перцентиля повлияет на снижение чувствительности к небольшим локальным пикам, а меньшие значения могут привести к усилению шума.
Таким образом, использование значения перцентиля в 45 позволяет осуществлять масштабирование спектра минимизируя влияние локальных выбросов.
Нормализованный по 45 перцентилю сигнал проходит этап автоматической регулировки усиления (AGC).
for j in range(jz):
s1_max = 20.0
s_max = np.max(sym_spec[j, :])
if s_max > s1_max:
sym_spec[j, :] *= s1_max / s_maxКорректировка амплитуд осуществляется относительно значения s1_max. Значение 20.0 (s1_max) подобрано эмпирически.
Синхронизация
Как было рассмотрено в предыдущей статье, в протоколе Q65 используется 65-FSK, из которых 64 тона для передачи данных и 1 пилот-тон выделен для синхронизации, занимающий 22 (Q65_SYNC, Q65_SYNC_TONES_COUNT) FSK-символа из 85 (Q65_TONES_COUNT) передаваемых.
Q65_SYNC = np.array([1, 9, 12, 13, 15, 22, 23, 26, 27, 33, 35, 38, 46, 50, 55, 60, 62, 66, 69, 74, 76, 85])
Q65_SYNC_TONES_COUNT = 22
Q65_TONES_COUNT = 85
Q65_TONES_CENTER = 43Номера символов, содержащих пилот-тон и параметры передаваемого пакета данных Q65.
Задача этапа синхронизации — обнаружить в спектрах по двум координатом (времени и частоте) расположение пилот-тона, относительно которого будет осуществляется детекция информационных символов.
Для обнаружения пилот-тона в Q65 применяется кросс-корреляция спектра.
Функция sync_ccf реализует алгоритм кросс-корреляции спектра sym_spec и пилот-тонов Q65_SYNC с учетом частотного дрейфа. Аргумент f0 задает нижнюю границу частот, относительно которой производится анализ.
def sync_ccf(self, sym_spec: npt.NDArray[np.float64], f0: float)-> typing.Tuple[int, int, float, float]:
jz, iz = sym_spec.shape
dec_df = 50
snf_a = f0 - dec_df
snf_b = f0 + dec_df
...Значение dec_df расширяет границы частот (нижнюю snf_a и верхнюю snf_b) поиска в узком диапазоне (Search Narrow Frequency), относительно частоты f0, с учетом допплеровского дрейфа.
max_drift = min(100 / self.df, 60) if self.f_max_drift else 0
...Переменная max_drift определяет максимальное значение частотного дрейфа, возникающего вследствие эффекта Доплера, для его последующей компенсации при вычислении корреляции. Управление осуществляется флагом self.f_max_drift.
bin_start = int(max(self.nf_a, 100) / self.df)
bin_end = int(min(self.nf_b, 4900) / self.df)
...Значения bin_start и bin_end — номера бинов (частотных участков) в спектре, рассчитываемые на значениях частот nf_a и nf_b. Максимальная полоса частот составляет диапазон от 100 Гц до 4900 Гц.
Примечание: нижняя граница частот определена как 100 Гц для исключения низкочастотных помех (в том числе наводки из осветительной сети), возникающих в цепях приема и передачи сигнала от приемника к компьютеру.
Верхняя граница 4900 Гц — захватывает полосу частот с расположенной в них FSK-тонами, с учетом возможного частотного дрейфа и исключает диапазон частот Найквиста.
sym_spec_avg = np.sum(sym_spec[:, bin_start:bin_end], axis=0)
...Массив sym_spec_avg определяет общие энергии частот в спектрах sym_spec, вычитаемых из фонового шума при расчете корреляции.
ccf_best = 0.0
best = 0
lag_best = 0
for i in range(bin_start, bin_end):
ccf_max_s = 0
ccf_max_m = 0
lag_peak_s = 0
for lag in range(self.lag1, self.lag2 + 1):
for drift in range(-max_drift, max_drift + 1):
...Циклы i, lag и drift определяют позиционирование окна в разрезе частот, с учетом дрейфа, и временных отсчетов, в рамках которого производится расчет коэффициента кросс-корреляции для пилот-тона.
ccf_t = 0.0
for kk in range(Q65_SYNC_TONES_COUNT):
k = Q65_SYNC[kk] - 1
zz = drift * (k - Q65_TONES_CENTER)
ii = i + zz // Q65_TONES_COUNT
if ii < 0 or ii >= iz:
continue
n = NSTEP * k
j = n + lag + self.j0
if j > -1 and j < jz:
ccf_t += sym_spec[j, ii]
ccf_t -= (Q65_SYNC_TONES_COUNT / jz) * sym_spec_avg[i - bin_start]
...Переменная ccf_t определяет коэффициент кросс-корреляции пилот-тона; накопленное в цикле значение представляет собой сумму амплитуд 22 (Q65_SYNC_TONES_COUNT) тонов, расположенных в спектре на позициях, согласно номерам из Q65_SYNC.
Итоговый коэффициент кросс-корреляции формируется путем вычитания энергий фоновых шумов sym_spec_avg.
if ccf_t > ccf_max_s:
ccf_max_s = ccf_t
lag_peak_s = lag
if ccf_t > ccf_max_m and drift == 0:
ccf_max_m = ccf_t
f = i * self.df
if ccf_max_s > ccf_best and snf_a <= f <= snf_b:
ccf_best = ccf_max_s
best = i
lag_best = lag_peak_s
...Значения переменных ccf_best и lag_best определяют позицию окна, в рамках которого было достигнуто максимальное значение кросс-корреляции пилот-тона в спектре сигнала, то есть определяет позицию, в которой пилот-тон точнее всего располагается.
i_peak = best - self.i0
j_peak = lag_best
f0 += i_peak * self.df
dt = j_peak * self.dt_step
return i_peak, j_peak, f0, dtРезультатом функции sync_ccf является четверка значений с позицией сигнала в разрезе частоты (i_peak), времени (j_peak); также корректируется фактическая базовая частота (f0), относительно которой производился расчет кросс-корреляции пилот-тона; значение dt определяет опережение/отставание сигнала принятого во времени, относительно времени начала анализа.
Извлечение информационных символов
После этапа синхронизации, в Q65 из спектра исключается пилот-тон, оставляя тем самым только информационные символы.
Q65_DATA_TONES_COUNT = Q65_TONES_COUNT - Q65_SYNC_TONES_COUNT
i_peak, j_peak, ccf_freq, time_d = self.sync_ccf(sym_spec, f0)
data_sym = self.get_data_sym(sym_spec, i_peak, j_peak, LL)Извлечение спектров символов данных осуществляется функцией get_data_sym.
def get_data_sym(self, sym_spec: npt.NDArray[np.float64], i_peak: int, j_peak: int, LL: int) -> npt.NDArray[
np.float64]:
jz, iz = sym_spec.shape
i1 = self.i0 + i_peak + self.q65_type - 64
i2 = i1 + LL
i3 = i2 - i1
data_sym = np.zeros((Q65_DATA_TONES_COUNT, i3), dtype=np.float64)
if i1 > 0 and i2 < iz:
j = self.j0 + j_peak - 8
n = 0
for k in range(Q65_TONES_COUNT):
j += 8
if self.sync[k] > 0.0:
continue
if j > 0 and j < jz:
data_sym[n, :i3] = sym_spec[j, i1:i2]
n += 1
bzap(data_sym)
return data_symФункция get_data_sym принимает на вход матрицу-спектр исходного сигнала с метками синхронизации, в качестве результата возвращает уплотненную матрицу-спектр, в котором присутствуют только спектры символов данных.
Фактически, get_data_sym реализует фильтрацию массива. Переменные i1, i2, i3 определяют границы частотных бинов, копируемых в результирующий массив.
На рисунке 5 представлена визуализация матрицы data_sym. Отличительной особенностью является отсутствие пилот-тона синхронизации.
Устранение локальных выбросов в спектре
В Q65 для подавления помех, которые на спектре проявляют себя в виде локальных выбросов, применяется статистическая коррекция спектра.
def bzap(data_sym: npt.NDArray[np.float64], threshold: int = 15):
_, num_bins = data_sym.shape
hist = np.zeros(num_bins, dtype=np.int64)
peaks = np.argmax(data_sym, axis=1)
np.add.at(hist, peaks, 1)
zap_bins = np.where(hist > threshold)[0]
data_sym[:, zap_bins] = 1.0Функция bzap принимает на вход матрицу-спектр, в котором формируется гистограмма hist по позициям максимальной амплитуды в спектре. Относительно значений гистограммы hist происходит эквализация амплитуд частотных бинов до значений уровня 1, статистическая повторяемость которых превышает порог threshold (по умолчанию 15).
Функция модифицирует исходную матрицу data_sym, не возвращая результата.
Замирания
Затухания (fading) или замирания — изменение мощности радиосигнала, возникающее из-за многолучевого распространения, изменениям условий прохождения на трассе или перемещением приемника в пространстве.
В зависимости от скорости изменения и причин возникновения, затухания разделяются на две большие категории: быстрые и медленные.
Быстрые затухания (fast fading) проявляют себя в виде резких флуктуаций (изменений амплитуд и фаз сигнала) на коротких временных интервалах, сравнимых с длительностью одного символа. Они могут быть частотно-селективными, то есть на одних частотах проявляться сильнее, чем на других. Это приводит к искажениям формы импульсов, приводя к межсимвольной интерференции.
Основные причины возникновения — многолучевое распространение радиосигнала, относительно быстрое перемещение приемника, передатчика, объектов, способных отражать сигнал, в том числе доплеровский эффект.
Медленные затухания (slow fading) проявляются как медленное изменение или спад среднего уровня сигнала сравнительно больших интервалах времени. Такой тип затуханий, как правило, не является частотно-селективным.
Причинами возникновения медленных затуханий, в основном, являются изменения свойств прохождения трассы, например при движении облаков, изменении погодных условий, а также из-за эффектов "затенения" сигнала.
На рисунке 6 приведены демонстрационные графики амплитуд сигнала при быстрых и медленных затуханиях.
Модель быстрых затуханий в Q65
Так как протокол Q65 рассчитан на работу в условиях достаточно сложных трасс, на которых может проявляться Доплеровский эффект и присущи быстрые затухания, на принимающей стороне спектр сигнал становится "размытым", символы выглядят как область частот, в пределах которой может находится символ (рисунок 7), и, как следствие, снижается вероятность детектирования символа.
Для борьбы с перечисленными эффектами, в Q65 применяется модель быстрых затуханий, основанная на расчете апостериорных вероятностей символов, задача которой заключается в устранении эффекта размытия области частот символов, тем самым обеспечивая более высокую чувствительность детектора при соотношении сигнал/шум (SNR) ниже -20дБ.
Протокол Q65 реализует модели затуханий Гаусса и Лоренца.
Модель Гаусса достаточно узкополосная и учитывает энергии вблизи центральной частоты и применима в случаях, когда в сигнале присутствуют незначительные искажения; может быть применена для обработки сигналов при тропосферном и дождевом рассеянии.
Модель Лоренца учитывает энергии частот, расположенных намного дальше от центральной частоты и позволяет нивелировать сильные искажения сигнала и эффект Доплера; применим на более сложных трассах, таких как EME.
В основном используется модель быстрых затуханий Лоренца.
Модели затуханий заданы таблично.
lorentz01 = np.array([
0.0214, 0.9107
])
lorentz02 = np.array([
0.0244, 0.9030
])
...
lorentz64 = np.array([
0.0008, 0.0008, 0.0008, 0.0009, 0.0009, 0.0009, 0.0010, 0.0010, 0.0010, 0.0011, 0.0011, 0.0012, 0.0012, 0.0012,
0.0013, 0.0013, 0.0014, 0.0014, 0.0015, 0.0016, 0.0016, 0.0017, 0.0018, 0.0019, 0.0020, 0.0021, 0.0022, 0.0023,
0.0024, 0.0025, 0.0027, 0.0028, 0.0030, 0.0031, 0.0033, 0.0035, 0.0038, 0.0040, 0.0043, 0.0046, 0.0049, 0.0052,
0.0056, 0.0061, 0.0066, 0.0071, 0.0077, 0.0084, 0.0091, 0.0100, 0.0109, 0.0120, 0.0132, 0.0145, 0.0159, 0.0175,
0.0192, 0.0209, 0.0228, 0.0246, 0.0264, 0.0279, 0.0291, 0.0299, 0.0301
])
fm_tab_lorentz = [
lorentz01, lorentz02, ..., lorentz64
]
...
gauss01 = np.array([
0.0296, 0.9101
])
gauss02 = np.array([
0.0350, 0.8954
])
...
gauss64 = np.array([
0.0028, 0.0029, 0.0030, 0.0031, 0.0032, 0.0034, 0.0035, 0.0036, 0.0037, 0.0039, 0.0040, 0.0041, 0.0043, 0.0044,
0.0046, 0.0047, 0.0048, 0.0050, 0.0051, 0.0053, 0.0054, 0.0056, 0.0057, 0.0059, 0.0060, 0.0062, 0.0063, 0.0065,
0.0066, 0.0068, 0.0069, 0.0071, 0.0072, 0.0074, 0.0075, 0.0077, 0.0078, 0.0079, 0.0081, 0.0082, 0.0083, 0.0084,
0.0086, 0.0087, 0.0088, 0.0089, 0.0090, 0.0091, 0.0092, 0.0093, 0.0094, 0.0094, 0.0095, 0.0096, 0.0097, 0.0097,
0.0098, 0.0098, 0.0099, 0.0099, 0.0099, 0.0099, 0.0100, 0.0100, 0.0100
])
fm_tab_gauss = [
gauss01, gauss02, ..., gauss64
]На рисунке 8 приведены графики весов моделей Гаусса и Лоренца.
Подробнее с матрицами весов моделей Гаусса и Лоренца можно ознакомиться здесь и здесь.
Примечание: формулы, по которым были рассчитаны коэффициенты, не приведены и, возможно, их значения были подобраны эмпирически.
Функция q65_intrinsics_fastfading (синоним q65_intrinsics_ff далее) реализует вычисление апостериорных вероятностей символов на основе модели быстрых затуханий.
class FadingModel(Enum):
Gaussian = 0
Lorentzian = 1
TS_QRA64 = 0.576
def q65_intrinsics_fastfading(
codec: Q65Codec,
input_energies: npt.NDArray[np.float64],
sub_mode: int,
B90Ts: float,
fading_model: FadingModel
) -> npt.NDArray[np.float64]:
...Функция принимает на вход input_energies — матрица спектров символов. Параметр sub_mode определяет режим протокола Q65 (0 — A, 1 — B, …, 4 — E); B90Ts — полоса пропускания сигнала 90% энергии (большие значения соответствуют большему растеканию энергии сигнала).
B90 = B90Ts / TS_QRA64
h_idx = int(np.log(B90) / np.log(1.09) - 0.499)
h_idx = min(63, max(0, h_idx))
if fading_model == FadingModel.Gaussian:
weights = fm_tab_gauss[h_idx]
elif fading_model == FadingModel.Lorentzian:
weights = fm_tab_lorentz[h_idx]
weights_count = len(weights)
...На основе значений параметров B90Ts и fading_model определяются весовые коэффициенты моделей.
EsNo_deg = 8.0 * np.log(B90) / np.log(240.0)
EsNo_metric = codec.decoderEsNoMetric * np.pow(10.0, EsNo_deg / 10.0)
...Корректировка оптимальной Es/No (отношение энергии на символ к спектральной плотности шума) характеристики декодера (EsNo_metric) при частотном дрейфе, исходя из предположения, что порог декодирования 50% равен Es/No(dB) = Es/No(AWGN)(dB) + 8*log(B90)/log(240)(dB).
При максимальном частотном дрейфе 240 Гц, вызываемым эффектом Доплера, наблюдается ухудшение Es/No приблизительно на 8 дБ (EsNo_deg) при измерениях на аддитивном Гауссовском шуме (AWGN).
M = codec.qra_code.alphabet_size
N = codec.qra_code.codeword_length
tone_span = 1 << sub_mode
sym_span = M * (2 + tone_span)
noise_var = np.mean(input_energies)
noise_var = noise_var / (1.0 + EsNo_metric / sym_span)
...Определение шумовой картины (noise variance).
codec.NoiseVar = noise_var
codec.EsNoMetric = EsNo_metric
codec.BinsPerTone = tone_span
codec.BinsPerSymbol = sym_span
codec.WeightsCount = weights_count
EsNo_gain = EsNo_metric * weights
weight = EsNo_gain / (EsNo_gain + 1) / noise_var
codec.FastFadingWeights = weight
sym_probs = np.zeros((N, M), dtype=np.float64)
tap_half = weights_count - 1
tap_center = 2 * tap_half
for n in range(N):
sym_start = M - weights_count + 1
for k in range(M):
log_prob = 0
for j in range(tap_half):
log_prob += weight[j] * (
input_energies[n, sym_start + j] + input_energies[n, sym_start + tap_center - j]
)
log_prob += weight[tap_half] * input_energies[n, sym_start + tap_half]
sym_probs[n, k] = log_prob
sym_start += tone_span
...Для каждого частотного бина в символе вычисляется логарифм вероятности (log_prob) с учетом весов модели затуханий.
log_prob_max = np.max(sym_probs[n])
sum_prob = 0
for k in range(M):
rel_log = sym_probs[n, k] - log_prob_max
rel_log = min(85.0, max(-85.0, rel_log))
rel_prob = np.exp(rel_log)
sym_probs[n, k] = rel_prob
sum_prob += rel_prob
sym_probs[n, :] *= 1.0 / sum_prob
return sym_probsЛогарифмические значения вероятностей нормализуются. В качестве результата функция q65_intrinsics_fastfading возвращает матрицу апостериорных вероятностей для каждого символа Q65.
На рисунке 9 приведен спектр сигнала, после применения к нему модели быстрых затуханий; частоты с символами более контрастны на фоне общей шумовой картины.
Структура Q65Codec представляет собой датакласс с параметрами QRA-декодера (в данной статье не рассматривается. В q65_intrinsics_fastfading происходит заполнение полей для последующего декодирования результатов вычисления апостериорных вероятностей.
Монитор
Часть вышеописанных функций объединена в класс Q65Monitor, задача которого предоставить интерфейс загрузки данных и декодирования сообщений.
Q65_TONES = 64
class Q65Monitor(AbstractMonitor):
__slots__ = ["sample_rate", "signal"]
CCF_OFFSET_R = 70
CCF_OFFSET_C = 7000
def __init__(self, eme_delay: bool, q65_type: typing.Literal[1, 2, 3, 4],
period: typing.Optional[typing.Literal[15, 30, 60, 120, 300]] = None):
self.signal = np.zeros(0, dtype=np.float64)
self.q65_type = q65_type
if period == 30:
self.sym_samps = 3600
elif period == 60:
self.sym_samps = 7200
elif period == 120:
self.sym_samps = 16000
else:
self.sym_samps = 3600
...Определение количества семплов на символ, в зависимости от периода протокола Q65.
self.fft_size = self.sym_samps
self.df = 12000.0 / self.fft_size
self.smooth = max(1, int(0.5 * self.q65_type ** 2))
...Параметры расчета БПФ; self.smooth — фактор сглаживания спектра, для снижения эффекта Доплера, определяется шириной полосы сигнала.
self.sym_steps = self.sym_samps // NSTEP
self.dt_step = self.sym_samps / (NSTEP * 12000.0)
self.lag1 = int(-1.0 / self.dt_step)
self.lag2 = int((5.5 if self.sym_samps >= 3600 and eme_delay else 1) / self.dt_step + 0.9999)
...Определение количество отсчетов на один символ и определение временных лагов для сигналов подверженных сильному Доплеровскому влиянию.
self.i0 = 0
self.j0 = int((1 if self.sym_samps >= 7200 else 0.5) / self.dt_step)
self.NN = 63
...Количество информационных символов в сигнале, определяет размер матриц сигнала.
self.max_iters = 100
self.nf_a = 214.333328
self.nf_b = 2000.000000
...Значения полей self.nf_a, self.nf_b определяют полосу частот, в которой будет производиться декодирование сигнала.
self.bw_a = 1
self.bw_b = 11
...Границы рассчета B90Ts — ширины спектра с 90% энергии, для расчета в модели быстрых замираний.
self.f_max_drift = False
self.q65_codec = q65_init()
self.sync = np.full(Q65_TONES_COUNT, -22.0 / 63.0, dtype=np.float64)
self.sync[Q65_SYNC - 1] = 1
...В self.sync записывается пилот-тон используемый в расчете кросс-корреляции при синхронизации.
def symbol_spectra(self, num_bins: int, num_times: int) -> npt.NDArray[np.float64]:
...
def sync_ccf(self, sym_spec: npt.NDArray[np.float64], f0: float) -> typing.Tuple[int, int, float, float]:
...
def get_data_sym(self, sym_spec: npt.NDArray[np.float64], i_peak: int, j_peak: int, LL: int) -> npt.NDArray[
np.float64]:
...
def decode_2(self, s3: npt.NDArray[np.float64], sub_mode: int, b90ts: float) -> typing.Optional[
typing.Tuple[float, npt.NDArray]
]:
s3prob = q65_intrinsics_ff(self.q65_codec, s3, sub_mode, b90ts, fading_model=FadingModel.Lorentzian)
...Метод decode_2 формирует матрицу вероятностей символов на основе модели быстрых затуханий, для ее последующей передачи в QRA-декодер (в рамках текущей статье не рассматривается).
def decode_q012(self, s3: npt.NDArray[np.float64]) -> typing.Optional[typing.Tuple[float, npt.NDArray]]:
if self.q65_type == 2:
sub_mode = 1
elif self.q65_type == 4:
sub_mode = 2
elif self.q65_type == 8:
sub_mode = 3
else:
sub_mode = 0
baud = 12000.0 / self.sym_samps
for bw in range(self.bw_a, self.bw_b + 1):
b90 = 1.72 ** bw
b90ts = b90 / baud
if (decoded := self.decode_2(s3, sub_mode, b90ts)) is not None:
EsNo_dB, payload = decoded
snr = EsNo_dB - dB(2500.0 / baud) + 3.0
return snr, payload
return None
...Метод decode_q012 реализует декодирование Q65 для режимов q1, q2 и q3. В цикле bw подбираются значения B90Ts для модели быстрых затуханий до тех пор, пока декодирование сигнала не даст положительный результат.
def decode_0(self, f0: float) -> typing.Tuple[float, float, npt.NDArray[np.int64], float]:
LL = Q65_TONES * (2 + self.q65_type)
txt = Q65_TONES * self.sym_samps / 12000.0 + (2 if self.sym_samps >= 6912 else 1)
iz = int(5000.0 / self.df)
jz = int(txt * 12000.0 / self.sym_steps)
...LL — число частотных бинов в спектре; txt — размер временного окна анализа сигнала; iz — общее количество частотных бинов в спектре с верхней границей в 5 КГц; jz — количество временных отсчетов в окне длительностью txt.
self.i0 = min(max(int(f0 / self.df), Q65_TONES), iz + Q65_TONES - LL)
...Значение self.i0 определяет центральную частоту, относительно которой будет производиться нормализация и выделение пилот-тона.
sym_spec = self.symbol_spectra(iz, jz)
for j in range(jz):
t_s = sym_spec[j, self.i0 - Q65_TONES:self.i0 - Q65_TONES + LL]
if (base := np.percentile(t_s, 45)) == 0:
base = 0.000001
sym_spec[j, :] /= base
for j in range(jz):
s1_max = 20.0
s_max = np.max(sym_spec[j, :])
if s_max > s1_max:
sym_spec[j, :] *= s1_max / s_max
i_peak, j_peak, ccf_freq, time_d = self.sync_ccf(sym_spec, f0)
data_sym = self.get_data_sym(sym_spec, i_peak, j_peak, LL)
snr, data = self.decode_q012(data_sym)
return time_d, ccf_freq, data, snr
...Метод decode_0 реализует верхнеуровневый интерфейс декодирования Q65, ранее рассмотренные построение спектра символов, рассчет 45 перцентиля, выделение символов из спектра и передачу его в декодер q1-q2-q3.
def decode(self, f0: int, **kwargs) -> typing.Generator[LogItem, None, None]:
dT, ccf_freq, payload, snr = self.decode_0(f0=f0)
yield LogItem(snr, dT, ccf_freq - f0, payload.tobytes(), 0)
def monitor_process(self, frame: npt.NDArray[np.int16]) -> None:
self.signal = np.concat([self.signal, frame])Заключение
В продолжении цикла статей, посвященных протоколу Q65, в этой статье были рассмотрены механизмы цифровой обработки сигналов, устранение влияния эффекта Доплера, механизм синхронизации пилот-тона, нормализация и фильтрация выбросов в спектре, определение вероятностей символов в условиях быстрых затуханий.
В статью не вошел алгоритм декодирования QRA, по причине достаточно большой темы Q-ary Repeat Accumulation кодирования, объем которой сопоставим с отдельной статьей.
Декодер Q65 может представлять интерес с точки зрения цифровой обработки сигналов в условиях Доплеровских искажений и быстрых затуханий.