Декодируем сигнал метеоспутника NOAA в 50 строчек кода

Привет, Хабр.

Продолжим тему декодирования различных видов радиосигналов. Одним из интересных форматов, разработанных еще в прошлом веке, является APT (Automatic Picture Transmission). Он используется для передачи изображений Земли из космоса, и что для нас гораздо более интересно и актуально, прием сигналов APT прост и доступен радиолюбителям.

Спутник NOAA © en.wikipedia.org/wiki/NOAA-15

Для тех, кому интересно, продолжение под катом.

Общая информация

Интересующие нас метеоспутники NOAA относятся к серии TIROS (Television InfraRed Observation Satellite), первый из которых был запущен еще в 1960 г. В настоящее время функционируют 3 спутника (NOAA-15, NOAA-18 и NOAA-19, самый старый из которых, NOAA-15 работает с 1998 г). Спутники вращаются вокруг Земли на высоте 850 км и делают один оборот примерно за 1.5 часа. На борту имеются различные сенсоры, нас же будет интересовать прием метеоизображений. И тут есть два варианта. Самый простой и общедоступный — аналоговый сигнал уже упомянутого формата APT на частоте 137 МГц. Кроме него, спутники передают изображения также и в формате высокого разрешения HRPT (High-resolution Picture Transmission) на частоте 1.7 ГГц. Декодеры HRPT доступны, но приема уже желательна поворотная антенна с неплохим усилением, что гораздо более сложно и дорого.

Согласно международным конвенциям, все метеорологические данные являются открытыми, и принимать сигналы NOAA может любой желающий. Чем мы сейчас и займемся.

По приему сигналов NOAA уже было несколько статей на Хабре (например 1, 2), но все они обычно сводятся к «подключили приемник, запустили программу, получили картинку», без особого разъяснения того, как же эта картинка передается. Попробуем опуститься на уровень ниже, и посмотреть, как же сигнал устроен внутри.

Прием

Самое первое и очевидное. Спутник не является геостационарным, он движется по небу, поэтому время приема нужно «ловить». Самое простое — воспользоваться онлайн-сервисом n2yo.com, который позволяет рассчитать время пролета для любого спутника. Вот пример для NOAA 19:

Хороший прием, когда спутник высоко над горизонтом, бывает примерно раз в день, сам пролет спутника по небу длится около 10 минут, так что может быть полезно настроить запись по расписанию в нужное время.

Разумеется, для приема радиосигнала нужен приемник. Подойдет самый дешевый RTL-SDR V3 за 35$, я использовал более качественный SDRPlay со следующими настройками:

Как можно видеть, я установил значение decimation на максимум, что позволяет получить максимальный динамический диапазон. Уровни усиления LNA и Gain стоит подобрать в зависимости от антенны. Спутники NOAA 15, 18 и 19 передают сигналы на частотах 137.620, 137.9125 и 137.100 МГц соответственно. Сам сигнал имеет ширину полосы около 50 КГц, и если все было сделано правильно, в заданное время сигнал должен появиться на спектре:

Интересно отметить наклон линий из-за эффекта Допплера — спутник пролетает мимо нас, а с учетом времени пролета, это еще одно хорошее доказательство шарообразности Земли 😉

Кстати о приеме, спутники NOAA хороши тем, что их прием весьма доступен для начинающих. Сигнал будет в принципе слышен на любую антенну, даже на «усы» от телевизора, но для хорошей картинки качество антенны весьма критично. С плохой антенны (как у меня:) картинка будет менее контрастной, но для демонстрации работы декодера этого достаточно. Специализированную антенну на 137 МГц можно купить или сделать из сантехнических и медных трубок, пример можно посмотреть здесь. Однако, тема статьи все же обработка сигналов, а не столярные работы, желающие могут изучить тему постройки антенн самостоятельно.

Итак, смотрим время пролета спутника, выбираем в программе режим записи FM с шириной полосы 50 КГц и записываем сигнал. Результатом должен быть WAV-файл продолжительностью около 10 минут, на слух должны быть четко различимы периодически повторяющиеся импульсы. Теперь можно приступать к декодированию.

Декодирование

Шаг 1. Загрузим файл с помощью numpy и выведем его на экран. Я вывожу только фрагмент с 20 по 21ю секунду, иначе отрисовка будет слишком долгой.

import scipy.io.wavfile as wav
import scipy.signal as signal
import numpy as np
import matplotlib.pyplot as plt


fs, data = wav.read('HDSDR_20201227_070306Z_137100kHz_AF.wav')
data_crop = data[20*fs:21*fs]

plt.figure(figsize=(12,4))

plt.plot(data_crop)
plt.xlabel("Samples")
plt.ylabel("Amplitude")
plt.title("Signal")
plt.show()

Должна быть видна четкая периодичность сигнала:

Шаг 2. Для ускорения декодирования уменьшим частоту дискретизации в 4 раза, отбросив лишние значения:

resample = 4
data = data[::resample]
fs = fs//resample

Шаг 3. Изображение передается в амплитудной модуляции, для конвертации в АМ достаточно применить преобразование Гильберта:

def hilbert(data):
    analytical_signal = signal.hilbert(data)
    amplitude_envelope = np.abs(analytical_signal)
    return amplitude_envelope

data_am = hilbert(data)

Можно вывести картинку на экран и убедиться, что мы получили огибающую сигнала:

Шаг 4. Вывод. Собственно, декодирование уже закончено. Сами данные передаются в аналоговом формате, так что цвет пиксела зависит от уровня сигнала. Мы можем «развернуть» картинку в 2D-изображение, из описания формата известно, что одна линия передается за 0.5с:

from PIL import Image

frame_width = int(0.5*fs)
w, h = frame_width, data_am.shape[0]//frame_width
image = Image.new('RGB', (w, h))

px, py = 0, 0
for p in range(data_am.shape[0]):
    lum = int(data_am[p]//32 - 32)
    if lum < 0: lum = 0
    if lum > 255: lum = 255
    image.putpixel((px, py), (0, lum, 0))
    px += 1
    if px >= w:
        if (py % 50) == 0:
            print(f"Line saved {py} of {h}")
        px = 0
        py += 1
        if py >= h:
            break

Функция putpixel не самый быстрый способ работы с изображениями, и код можно ускорить раз в 10, используя numpy.reshape и Image.fromarray, но построчный способ более наглядный. Для конвертации амплитуды сигнала в диапазон яркости 0..255, значения делятся на 32, для другой антенны возможно, значение придется изменить.

Для удобства просмотра растягиваем картинку по вертикали и выводим на экран:

image = image.resize((w, 4*h))
plt.imshow(image)
plt.show()

Если все было сделано правильно, мы должны получить изображение примерно такого вида:

Теплый и ламповый зеленый цвет был выбран просто для красоты, желающие могут изменить цветовую гамму по своему усмотрению, изменяя параметры функции putpixel. Что мы видим на экране? В формате APT передаются два канала. Дальний ИК передается на одной половине кадра, ближний/средний ИК передается на другой, режим выбирается в зависимости от того, ночную или дневную картинку передает спутник. В изображении также передаются маркеры синхронизации и телеметрия, желающие могут посмотреть описание формата APT более подробно. Более продвинутые программы используют эти маркеры для выравнивания картинки, но могут не работать на слабых сигналах. Вышеприведенный код никогда не теряет синхронизацию, т.к. её здесь просто нет, даже самый слабый и зашумленный сигнал будет виден, хоть и менее контрастно.

Заключение

Как уже говорилось ранее, в мире осталось не так уж много действующих систем связи, сигнал которых можно декодировать с помощью 20 строчек кода. Спутникам NOAA примерно 20 лет, и когда они выйдут из строя, новый формат уже скорее всего будет цифровым и более сложным. Так что желающим попробовать что-то простое и несложное для декодирования, можно посоветовать поторопиться.

Разумеется, необязательно декодировать сигналы NOAA таким способом, существуют и готовые декодеры, которые обладают большими возможностями, вроде нанесения контуров стран и городов, формирования псевдоцветов из монохромных каналов и пр. Но сделать это самому с нуля имхо гораздо интереснее.

Как обычно, всем удачных экспериментов.