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

Во второй части были рассмотрены практические аспекты использования SDR. В этой части мы разберемся, как принять данные метеоспутника NOAA с помощью Python и недорогого (30$) приемника RTL-SDR. Рассмотренный код будет работать везде — на Windows, OSX, Linux и даже на Raspberry Pi.



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

SoapySDR


Производителей различных SDR-устройств довольно много, и поддерживать каждое по отдельности было бы весьма неудобно, да и дорого в плане покупки «железа» для тестирования. В принципе, для унифицированного доступа существует две библиотеки, которые по сути стали стандартом. Первая — это уже довольно старый интерфейс ExtIO DLL, которому наверно не менее 10 лет, вторая — более современная библиотека SoapySDR, которую мы и рассмотрим.

SoapySDR — это набор кроссплатформенных библиотек, написанных на C++, предоставляющих унифицированный доступ к SDR-устройствам, как к приемникам, так и трансиверам. Если производитель делает такой интерфейс, то его устройство «автоматом» будет работать с достаточно большим числом популярных программ (GQRX, GNU Radio, CubicSDR и пр). Практически все адекватные производители, кроме некоторых, (пользуясь случаем, передаю привет компании ЕЕ) имеют поддержку SoapySDR, список поддерживаемых устройств можно посмотреть на странице проекта. Как можно видеть, он довольно большой, и включает HackRF, USRP, SDRPlay, LimeSDR, RTL-SDR, Red Pitaya и многие другие.

Библиотека SoapySDR является кроссплатформенной, т.е. написанный под неё код будет работать под Windows, OSX, Linux, и даже на Raspberry Pi. Для Windows нужные библиотеки входят в состав пакета PothosSDR, для остальных платформ скомпилировать SoapySDR придется самостоятельно. Нужно скомпилировать две части — собственно библиотеку, и «плагин» для нужного приемника, в нашем случае это будет SoapyRTLSDR (под Windows библиотеку тоже можно собрать из исходников, для этого понадобится Visual Studio, Cmake, и SWIG). Теперь все готово, и можно писать код.

Импортируем библиотеку и получаем список приемников:

from __future__ import print_function
import SoapySDR

# Enumerate devices
print("SDR devices:")
for d in SoapySDR.Device.enumerate(''):
    print(d)
print()

Подключаем приемник, запускаем код и видим список устройств, среди которых есть наш rtlsdr.



Остальные устройства это звуковые карты, как мы помним, исторически первые SDR работали именно через линейный вход ПК, и библиотека их тоже поддерживает. Получаем информацию об устройстве — число доступных каналов, частотный диапазон и пр:

soapy_device = "rtlsdr"
device = SoapySDR.Device(dict(driver = soapy_device))

channels = list(range(device.getNumChannels(SoapySDR.SOAPY_SDR_RX)))
print("Channels:", channels)

ch = channels[0]

sample_rates = device.listSampleRates(SoapySDR.SOAPY_SDR_RX, ch)
print("Sample rates:\n", sample_rates)

bandwidths = list(map(lambda r: int(r.maximum()), device.getBandwidthRange(SoapySDR.SOAPY_SDR_RX, ch)))
print("Bandwidths:\n", bandwidths)

print("Gain controls:")
for gain in device.listGains(SoapySDR.SOAPY_SDR_RX, ch):
    print("  %s: %s" % (gain, device.getGainRange(SoapySDR.SOAPY_SDR_RX, ch, gain)))

frequencies = device.listFrequencies(SoapySDR.SOAPY_SDR_RX, ch)
print("Frequencies names:", frequencies)

frequency_name = frequencies[0]
print("Frequency channel name:", frequency_name)

print("Frequency range:", device.getFrequencyRange(SoapySDR.SOAPY_SDR_RX, ch, frequency_name)[0])

Запускаем программу и видим информацию о приемнике:



Мы видим что приемник имеет один входной канал с названием «RF», возможные частоты дискретизации [250000.0, 1024000.0, 1536000.0, 1792000.0, 1920000.0, 2048000.0, 2160000.0, 2560000.0, 2880000.0, 3200000.0] и частотный диапазон 24МГц-1.7ГГц.

Лайфхак — те же данные можно получить и из командной строки, набрав команду SoapySDRUtil --probe=«driver=rtlsdr».

Зная это, мы можем записать поток данных в WAV. Как говорилось в предыдущей части, данные с SDR представлены потоком сигналов, называемых I и Q, представляющих собой отсчеты с АЦП, грубо их можно представить как RAW-данные с фотокамеры. Кому интересно подробнее, могут почитать например здесь. Для нас достаточно знать, что мы эти данные можем записать, а другие SDR-программы потом могут с ними работать.

Сама по себе запись довольно проста — функция readStream заполняет буфер если есть данные, если данных еще нет, то вернется -1. Ниже показан код записи 10 отсчетов (несущественные части кода опущены).

device.setFrequency(SoapySDR.SOAPY_SDR_RX, channel, "RF", frequency)
device.setGain(SoapySDR.SOAPY_SDR_RX, channel, "TUNER", gain)
device.setGainMode(SoapySDR.SOAPY_SDR_RX, channel, False)
device.setSampleRate(SoapySDR.SOAPY_SDR_RX, channel, sample_rate)

# Number of blocks to save
block, max_blocks = 0, 10

block_size = device.getStreamMTU(stream)
print("Block size:", block_size)

buffer_format = np.int8
buffer_size = 2*block_size  # I+Q
buffer = np.empty(buffer_size, buffer_format)

while True:
    d_info = device.readStream(stream, [buffer], buffer_size)
    if d_info.ret > 0:
        wav.write(buffer[0:2*d_info.ret])
        print("Bytes saved:", 2*d_info.ret)
        block += 1
        if block > max_blocks:
            break

Результат на скриншоте:



Как можно видеть, мы получаем от устройства блоки данных, размер одного блока составляет 131072 байта, что при частоте дискретизации 250000 дает нам длительность примерно полсекунды. В общем, кто работал ранее со звуковыми картами в Windows, найдет много общего.

Для теста записываем файл и проверяем, что все нормально — его можно воспроизвести в SDR#. Тут есть еще одна хитрость — чтобы SDR# корректно показал частоты станций, имя файла нужно записать в формате, совместимом с HDSDR, вида «HDSDR_20190518_115500Z_101000kHz_RF.wav» (как нетрудно догадаться, в начале идет дата и время в GMT, затем частота в килогерцах). Это нетрудно записать на Python:

frequency = 101000000
file_name = "HDSDR_%s_%dkHz_RF.wav" % (datetime.datetime.utcnow().strftime("%Y%m%d_%H%M%SZ"), frequency/1000)

Для начала проверяем на FM-диапазоне. Все нормально, станции видно, музыка играется, RDS работает.



Можно приступать к записи NOAA.

Прием NOAA


Итак, у нас есть приемник и есть программа записи. Нам будут интересны метеоспутники NOAA 15, NOAA 18 и NOAA 19, передающие изображения поверхности Земли на частотах 137.620, 137.9125 и 137.100МГц. Основная сложность здесь в том, что нужно «поймать» момент, когда спутник пролетает над нами. Узнать время пролета можно онлайн по ссылкам https://www.n2yo.com/passes/?s=25338, https://www.n2yo.com/passes/?s=28654 и https://www.n2yo.com/passes/?s=33591 соответственно.



Чтобы не сидеть у компьютера, добавим в программу ожидание нужного времени. Это позволит также запускать программу на Raspberry Pi, без дисплея и клавиатуры.

import datetime

def wait_for_start(dt):
    # Wait for the start
    while True:
        now = datetime.datetime.now()
        diff = int((dt - now).total_seconds())
        print("{:02d}:{:02d}:{:02d}: Recording will be started after {}m {:02d}s...".format(now.hour, now.minute, now.second, int(diff / 60), diff % 60))
        time.sleep(5)
        if diff <= 1:
            break

wait_for_start(datetime.datetime(2019, 5, 18, 21, 49, 0))


Кстати, чтобы запустить скрипт на Raspberry Pi и оставить его работать после закрытия консоли, нужно ввести команду «nohup python recorder.py &».

Все готово, запускаем скрипт и можем заниматься другими делами, запись длится примерно 20 минут. Параллельно может возникнуть вопрос — можно ли увидеть пролет спутника невооруженным глазом? Согласно таблице, его максимальная яркость составляет порядка 5.5м звездной величины, предел человеческого глаза в идеальных условиях 6м. Т.е. при реально темном небе, далеко за городом, пролет спутника NOAA теоретически заметить можно, в городе конечно, шансов нет (как писали на Хабре, уже выросло поколение людей, никогда в жизни не видевших Млечный Путь).

Результатом работы скрипта является записанный wav-файл, его спектр показан на скриншоте.



Мы видим вполне различимый сигнал, хотя конечно со специальной антенной для приема NOAA качество было бы гораздо лучше. Формат сигнала называется APT (Automatic Picture Transmission), из него можно получить изображение земной поверхности, если кому интересно, можно отдельно рассмотреть его декодирование. Но есть разумеется, и готовые программы, декодировать такие сигналы можно с помощью WxToImg или MultiPSK.

Интересно видеть допплеровский сдвиг на спектре, который происходит оттого, что спутник пролетает мимо нас. Наверное, несложно посчитать и его скорость, зная время и сдвиг частоты в герцах.

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

Исходный код
from __future__ import print_function
import SoapySDR
import numpy as np
import struct
import sys
import time
import datetime


def wait_for_start(dt):
    # Wait for the start
    while True:
        now = datetime.datetime.now()
        diff = int((dt - now).total_seconds())
        print("{:02d}:{:02d}:{:02d}: Recording will be started after {}m {:02d}s...".format(now.hour, now.minute, now.second, int(diff / 60), diff % 60))
        time.sleep(5)
        if diff <= 1:
            break


def sdr_enumerate():
    # Enumerate SDR devices
    print("SDR devices:")
    for d in SoapySDR.Device.enumerate(''):
        print(d)
    print()


def sdr_init():
    soapy_device = "rtlsdr"
    device = SoapySDR.Device({"driver": soapy_device})

    channels = list(range(device.getNumChannels(SoapySDR.SOAPY_SDR_RX)))
    print("Channels:", channels)

    ch = channels[0]

    sample_rates = device.listSampleRates(SoapySDR.SOAPY_SDR_RX, ch)
    print("Sample rates:\n", sample_rates)

    print("Gain controls:")
    for gain in device.listGains(SoapySDR.SOAPY_SDR_RX, ch):
        print("  %s: %s" % (gain, device.getGainRange(SoapySDR.SOAPY_SDR_RX, ch, gain)))

    frequencies = device.listFrequencies(SoapySDR.SOAPY_SDR_RX, ch)
    print("Frequencies names:", frequencies)

    frequency_name = frequencies[0]
    print("Frequency channel name:", frequency_name)

    print("Frequency range:", device.getFrequencyRange(SoapySDR.SOAPY_SDR_RX, ch, frequency_name)[0])

    print()
    return device


def sdr_record(device, frequency, sample_rate, gain, blocks_count):
    print("Frequency:", frequency)
    print("Sample rate:", sample_rate)
    print("Gain:", gain)

    channel = 0  # Always for RTL-SDR
    device.setFrequency(SoapySDR.SOAPY_SDR_RX, channel, "RF", frequency)
    device.setGain(SoapySDR.SOAPY_SDR_RX, channel, "TUNER", gain)
    device.setGainMode(SoapySDR.SOAPY_SDR_RX, channel, False)
    device.setSampleRate(SoapySDR.SOAPY_SDR_RX, channel, sample_rate)

    data_format = SoapySDR.SOAPY_SDR_CS8 # if 'rtlsdr' in soapy_device or 'hackrf' in soapy_device else SoapySDR.SOAPY_SDR_CS16
    stream = device.setupStream(SoapySDR.SOAPY_SDR_RX, data_format, [channel], {})
    device.activateStream(stream)

    block_size = device.getStreamMTU(stream)
    print("Block size:", block_size)
    print("Data format:", data_format)
    print()

    # IQ: 2 digits ver variable
    buffer_format = np.int8
    buffer_size = 2 * block_size # I + Q samples
    buffer = np.empty(buffer_size, buffer_format)

    # Number of blocks to save
    block, max_blocks = 0, blocks_count

    # Save to file
    file_name = "HDSDR_%s_%dkHz_RF.wav" % (datetime.datetime.utcnow().strftime("%Y%m%d_%H%M%SZ"), frequency/1000)
    print("Saving file:", file_name)
    with open(file_name, "wb") as wav:
        # Wav data info
        bits_per_sample = 16
        channels_num, samples_num = 2, int(max_blocks * block_size)
        subchunk_size = 16  # always 16 for PCM
        subchunk2_size = int(samples_num * channels_num * bits_per_sample / 8)
        block_alignment = int(channels_num * bits_per_sample / 8)

        # Write RIFF header
        wav.write('RIFF'.encode('utf-8'))
        wav.write(struct.pack('<i', 4 + (8 + subchunk_size) + (8 + subchunk2_size)))  # Size of the overall file
        wav.write('WAVE'.encode('utf-8'))
        # Write fmt subchunk
        wav.write('fmt '.encode('utf-8'))  # chunk type
        wav.write(struct.pack('<i', subchunk_size))  # subchunk data size (16 for PCM)
        wav.write(struct.pack('<h', 1))  # compression type 1 - PCM
        wav.write(struct.pack('<h', channels_num))  # channels
        wav.write(struct.pack('<i', int(sample_rate)))  # sample rate
        wav.write(struct.pack('<i', int(sample_rate * bits_per_sample * channels_num/ 8)))  # byte rate
        wav.write(struct.pack('<h', block_alignment))  # block alignment
        wav.write(struct.pack('<h', bits_per_sample))  # sample depth
        # Write data subchunk
        wav.write('data'.encode('utf-8'))
        wav.write(struct.pack('<i', subchunk2_size))
        while True:
            d_info = device.readStream(stream, [buffer], buffer_size)
            if d_info.ret > 0:
                data = buffer[0:2*d_info.ret]
                fileData = data
                if data_format == SoapySDR.SOAPY_SDR_CS8:
                   fileData = data.astype('int16')
                wav.write(fileData)
                print("Block %d saved: %d bytes" % (block, 2*d_info.ret))
                block += 1
                if block > max_blocks:
                    break

    device.deactivateStream(stream)
    device.closeStream(stream)

if __name__ == "__main__":
    print("App started")

    # Forecast for active NOAA satellites
    # NOAA 15: 137.620, https://www.n2yo.com/passes/?s=25338
    # NOAA 18: 137.9125, https://www.n2yo.com/passes/?s=28654
    # NOAA 19: 137.100, https://www.n2yo.com/passes/?s=33591
    # Wait for the start: 18-May 21:49 21:49:
    wait_for_start(datetime.datetime(2019, 5, 18, 21, 49, 0))

    device = sdr_init()

    t_start = time.time()

    sdr_record(device, frequency=137912500, sample_rate=250000, gain=35, blocks_count=2100)

    print("Recording complete, time = %ds" % int(time.time() - t_start))
    print()


Плюс SoapySDR в том, что эта же программа с минимальными изменениями будет работать и с другими приемниками, например с SDRPlay или HackRF. Ну и про кроссплатформенность тоже уже упоминалось.

Если у читателей еще остался интерес к теме радиоприема, можно рассмотреть пример использования SDR с GNU Radio на примере создания нескольких виртуальных приемников на базе одного «железного».

Комментарии (14)


  1. ptica_filin
    19.05.2019 00:34

    Какой антенной принимали?


    1. DmitrySpb79 Автор
      19.05.2019 00:36

      Антенна никакая по сути, просто штырек на окне.


      1. I-denis
        19.05.2019 10:16

        на штырек у окна сигнал прилетает с периодическими замираниями, так что декодировать более менее получается отдельные строки — в итоге будет зашумленная зебра, окно лучше всего восточное или западное, и выбирать тот пролет, что бы спутник был в максимальном возвышении 60-70 градусов, иначе дом экранирует.

        и.да, для метеор м2 существует скрипт под gnu-radi, а там уже lrpt


        1. DmitrySpb79 Автор
          19.05.2019 13:40

          Спасибо. Штырек конечно плохо для NOAA, там и поляризация не совпадает, не говоря о КСВ.

          «Метеор» пока ни разу не пробовал принять, все как-то некогда было, да и я не настолько фанат «космоса», чтобы всерьез этим заморачиваться. Мне больше интересна математика чем столярные работы с крепежами и трубками ;)

          Если есть хорошие записи с NOAA в wav, можно поэскпериментировать с Python-декодером на них.


  1. elenbert
    19.05.2019 01:10

    Прошлым летом занимался активным приемом сигналов спутников NOAA-18 и NOAA-19.
    Антенна — квадрифиляр из трубок от кондиционера. Приемник — такой же RTL.
    Софт — Gnupredict на Linux и GQRX. Первая программа удобно отображает пролеты всех спутников и умеет управлять приемником GQRX, корректируя частоту, в зависимости от допплеровского сдвига.
    Постобработка записанного сигнала в Audacity (стерое в моно, компрессия и ресемплинг) и получение финального изображения в WxToImg.

    Некоторые результаты:






    1. wtpltd
      19.05.2019 01:46

      Я понимаю, что «некоторые», это «лучшее». Но впечатляет. Если бы вы поделились своим опытом, думаю, интересно было бы не только мне. В любом случае — респект!


      1. elenbert
        19.05.2019 10:33

        Спасибо. Думаю оформлю статью, как раз есть материал в виде фотографий, скриншотов и даже видео.


    1. I-denis
      19.05.2019 10:26

      gnupredict управляет, а каким образом — под рукой линукса нет сейчас посмотреть.
      я, в свое время использовал, wxtoimg, qth антену, lna, и самодельный супергетеродин с полосой пропускания 30кгц. частота приема и время приема управлялось от wxtoimg, доплеровский сдвиг компенсировался програмно в wxtoimg.


      1. elenbert
        19.05.2019 10:37

        В gnupredict есть модуль «control radio», а в gqrx есть опция «remote control».
        Активируются оба модуля, программы соединяются через локальный сокет и дальше в отдельном окне контроллера радио можно задать параметры управления — настроить приемник на заданную частоту (можно выбрать из списка избранных спутников-транспондеров) и начать трекать частоту, в этом режиме как раз и происходит коррекция в зависимости расстояния от спутника.

        Вот так это выглядит


        1. I-denis
          19.05.2019 10:51

          спасибо!


  1. Balling
    19.05.2019 13:26

    Слушайте, я хотел посоветовать вам одну статью. В свое время она меня так поразила. Там про себестоимость трансивера SDR. Это ж просто ужас какой-то. habr.com/ru/company/zeptobars/blog/359144


    1. DmitrySpb79 Автор
      19.05.2019 13:35

      Да, интересно.

      Но по сути ничего удивительного, может для каждого чипа себестоимость и 1$, но чтобы его сделать, нужно 5 лет работы лаборатории с бюджетом в миллион, соответственно, продавать этот чип будут за 50$ чтобы окупить затраты, да и на будущие разработки надо что-то отложить, иначе прогресса не будет.


  1. R820T2
    19.05.2019 23:03

    Дмитрий, про GnuRadio напишите, желательно подробнее и для неподготовленных. Я ссылки на все статьи добавляю в группу


  1. afiskon
    20.05.2019 13:03

    Интересная статья, благодарю. Пишите еще!