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