Software Defined Radio — как это работает? Часть 3

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

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

    Похожие публикации

    Комментарии 20

      0

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

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

          и.да, для метеор м2 существует скрипт под gnu-radi, а там уже lrpt
            0
            Спасибо. Штырек конечно плохо для NOAA, там и поляризация не совпадает, не говоря о КСВ.

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

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





          +1
          Я понимаю, что «некоторые», это «лучшее». Но впечатляет. Если бы вы поделились своим опытом, думаю, интересно было бы не только мне. В любом случае — респект!
            +2
            Спасибо. Думаю оформлю статью, как раз есть материал в виде фотографий, скриншотов и даже видео.
            0
            gnupredict управляет, а каким образом — под рукой линукса нет сейчас посмотреть.
            я, в свое время использовал, wxtoimg, qth антену, lna, и самодельный супергетеродин с полосой пропускания 30кгц. частота приема и время приема управлялось от wxtoimg, доплеровский сдвиг компенсировался програмно в wxtoimg.
              +1
              В gnupredict есть модуль «control radio», а в gqrx есть опция «remote control».
              Активируются оба модуля, программы соединяются через локальный сокет и дальше в отдельном окне контроллера радио можно задать параметры управления — настроить приемник на заданную частоту (можно выбрать из списка избранных спутников-транспондеров) и начать трекать частоту, в этом режиме как раз и происходит коррекция в зависимости расстояния от спутника.

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

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

              Но по сути ничего удивительного, может для каждого чипа себестоимость и 1$, но чтобы его сделать, нужно 5 лет работы лаборатории с бюджетом в миллион, соответственно, продавать этот чип будут за 50$ чтобы окупить затраты, да и на будущие разработки надо что-то отложить, иначе прогресса не будет.
              0
              Дмитрий, про GnuRadio напишите, желательно подробнее и для неподготовленных. Я ссылки на все статьи добавляю в группу
                0
                Интересная статья, благодарю. Пишите еще!
                  0
                  Добрый день. Отличная серия статей. Спасибо! Но Не до конца понятен процесс инсталляции PothosSDR. После его установки надо каким-то образом дать интерпретатору понять где он находится. Мой интерпретатор его не видит. Соответственно import SoapySDR выдаёт ошибку.
                    0
                    В инсталляторе есть флажок «Advanced/Python bindings».

                    В совсем крайнем случае можно собрать SoapySDR из исходников при помощи Visual Studio и CMake.
                    github.com/pothosware/SoapySDR
                      0

                      При установке флажок ставил. Но нет эффекта. ((

                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                Самое читаемое