24.10.2020 11:00
DmitrySpb79 17 4100 Источник
Привет, Хабр.

Читатели старшего поколения, заставшие дома радиоприемники средних, длинных и коротких волн, наверное помнят, что разные длины волн по-разному распространяются в различное время суток. Но как действительно это работает?



Я покажу как с помощью SDR-приемника и 50 строк кода на Python получить визуализацию сигналов радиостанций с точностью до долей герца, и увидеть довольно-таки любопытные атмосферные эффекты.

Продолжение под катом.

Общий принцип


Большинство АМ-радиостанций работает круглосуточно, что позволяет весьма наглядно изучить передаваемые ими радиосигналы. Для этого мы запишем сигнал радиостанции в формате WAV и построим его спектр при помощи FFT (Быстрого Преобразования Фурье). FFT позволяет из сигнала во «временной области» (time domain) получить изображение в «частотной области» (frequency domain), проще говоря, спектр сигнала. Чем больше размер окна преобразования, тем большее разрешение по частоте мы можем получить.

Как известно, сигналы АМ-станций выглядит в эфире следующим образом:



Собственно содержимое передачи нас интересовать не будет (кто-то еще слушает вообще радио?), для нас важна несущая — центр сигнала. Она является хорошим маркером, по которой удобно контролировать сигнал станции на спектре.

Запись


Для записи нам потребуется практически любой радиоприемник, умеющий принимать сигнал в формате боковой полосы (USB, Upper Side Band). SDR в этом плане наиболее удобен, но теоретически, даже обычный китайский Tecsun/Degen может подойти, если подключить его к линейному входу ПК.

Важный момент: станция передает в АМ, но но нам нужен сигнал до демодуляции, поэтому в настройках нужно выставить режим USB а не AM. Выберем ширину записи так, чтобы несущая радиостанции попадала в середину спектра. Это важно, т.к. при обработке мы будем вырезать из спектра именно середину.

Чем длиннее запись, тем интереснее результаты, основное ограничение тут в размере получаемого файла. Я использовал прием в режиме USB с полосой 4 КГц, формат записи WAV 8000 семплов/с, частоту SDR выбрал так, чтобы несущая частота радиостанции была в середине полосы фильтра. При таких настройках запись длительностью 24ч занимает в WAV около 1.3 ГБайт (на всякий случай напомню, что MP3 или другое сжатие с потерями для анализа сигналов использовать нельзя). Мои настройки HDSDR при записи выглядят так (важные моменты обозначены цифрами 1, 2, 3):



Обработка


Исходный код на языке Python приведен под спойлером. Мы последовательно читаем данные из WAV-файла, применяем к каждому блоку данных FFT + оконную функцию, и сохраняем результат в виде изображения. Яркость спектра можно варьировать в коде с помощью изменения параметра k_brightness, размер блока FFT передается в виде параметра командной строки. При использовании больших размеров FFT, например 4194304, мы не можем создать изображение такого же размера, поэтому из спектра сохраняется только центр (именно поэтому важно, чтобы несущая была по центру, хотя при желании смещение можно скорректировать вручную в коде).

spectrum.py 
import numpy as np
import matplotlib.pyplot as plt
import wave
from PIL import Image
import sys
import time


def wav_process(filename: str, fft_size: int):
    # FFT size can be a number like 1024 or power of 2, like 2**20
    if fft_size < 512:
        fft_size = 2**fft_size

    w = wave.open(filename, 'r')
    num_columns = w.getnframes() // fft_size
    if num_columns > 4096:
        num_columns = 4096

    # Spectum parameters
    width, height = num_columns, 2048
    k_brightness = 16
    palette_r, palette_g, palette_b = 4, 1, 4

    # Output image
    data_out = np.zeros([height, width, 3])
    if fft_size < 2*height:
        fft_size = 2*height

    print("WAV Sample width:", w.getsampwidth())
    print("WAV Frame rate:", w.getframerate())
    print("Image size: {}x{}".format(width, height))
    print("Columns max. count:", w.getnframes() // fft_size)
    print("FFT Size:", fft_size)
    print("Spectrum width, Hz:", height*w.getframerate()/fft_size)
    print()

    pos_top = 0
    if fft_size > 2 * height:
        # Auto align vertical
        spectrum_center_y = height//2
        pos_top = spectrum_center_y*fft_size//(2*height) - height//2

    # Draw
    for x in range(num_columns):
        if (x % 50) == 0:
            print("Column {} of {}".format(x, num_columns))

        data = w.readframes(fft_size)
        raw_data = np.frombuffer(data, dtype=np.int16)

        if raw_data.shape[0] != fft_size:
            break
        data_h = np.hanning(fft_size)*raw_data
        raw_fft = np.fft.fft(data_h, n=fft_size, norm="ortho")[1:fft_size]
        raw_abs = np.absolute(raw_fft)
        raw_int = (raw_abs/k_brightness)  # (raw_abs/k).astype(int)

        for p1 in range(height):
            lum_val = raw_int[pos_top + p1]
            v = lum_val if lum_val < 255 else 255
            data_out[p1][x] = [v/palette_r, v/palette_g, v/palette_b]

    # Numpy array to PIL image
    img = Image.fromarray(np.uint8(data_out), 'RGB')
    # img.save('spectrum_{}x{}.png'.format(width, height))
    
    # Display
    plt.figure(tight_layout=True)
    time_total = num_columns*fft_size//w.getframerate()
    hz_total = height*w.getframerate()//fft_size
    extent = [0, time_total/3600, -hz_total//2, hz_total//2]
    plt.xlabel('Time, Hr')
    plt.ylabel('Freq, Hz')
    plt.imshow(img, extent=extent, aspect='auto')
    plt.show()


if __name__ == '__main__':
    print("WAV Spectrum Builder v0.1. (c) DmitrySpb79 for habr.com")
    if len(sys.argv) < 2:
        print("Usage: python spectrum.py file.wav [fft_size]")
        exit()

    wav_process(sys.argv[1], fft_size=int(sys.argv[2]) if len(sys.argv) == 3 else 0)

    print("App done")

Запуск программы производится из командной строки, в качестве параметров указываются имя файла и размер блока FFT, который должен быть равен степени двух. Например:

python spectrum.py C:\HDSDR\HDSDR_20201023_083057Z_6068kHz_AF.wav 524288

Чтобы читателям не считать степени двойки вручную, приведу значения здесь: 4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576, 2097152, 4194304, 8388608, 16777216,… Чем больше размер блока, тем выше разрешение по частоте, но ниже разрешение по времени, соответственно больше отсчетов требуется для отображения результата. Так, при размере блока FFT в 4194304, мы получаем разрешение по вертикали 0.002 Гц на пиксел, но всего лишь 70 пикселов спектра из 8-часовой записи. В коде нет никаких оптимизаций, возможно картинку можно улучшить перекрытиями спектра или варьированием вида оконной функции, но в разы лучше вряд ли будет, по сути все ограничивается длиной записи.

Результаты


Несколько примеров работы программы.

Запись на частоте 894 КГц. Небольшой размер блока FFT (4096 отсчетов), файл с большой шириной полосы записи, мы видим сигналы АМ-станций практически в том виде, в каком они передаются в эфире. Собственно, на картинке можно наблюдать сразу две станции, стандартный шаг в АМ между станциями 9 КГц.



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

Рассмотрим запись с большей детализацией. Та же частота, по центру картинки расположена несущая, размер FFT 262144 семпла. На экране ~5ч записи.



Запись производилась ночью, видно что к утру сигналы станции стали слабее и затем совсем пропали. Также видны другие, даже более интересные эффекты. Рядом с несущей, около нулевой частоты, видны другие линии, расположенные совсем близко по частоте. Что это, я не знаю. Могу предположить, что это сигналы радиостанций других регионов, которые слишком далеко для полноценного приема, но их несущая слегка «пробивается».

Одна такая линия начинается на отметке «1ч» и пропадает на отметке «4.4ч», при увеличении вполне четко виден симметричный спектр АМ-сигнала:



Скорее всего, здесь мы видим сигнал какого-то далекого передатчика, который принимался отражаясь от ионосферы, затем прохождение завершилось. Он расположен на той же частоте 894 КГц, но настройка не является идеальной, и реальная частота отличается от первого передатчика примерно на 20 Гц, что и позволило увидеть их раздельно. На слух уловить такое, разумеется, невозможно. Как и в астрономии, такие эффекты видны лишь при накоплении сигнала очень большой длительностьи.

И последний пример. Короткие волны, станция на частоте 6070 КГц. Запись продолжительностью сутки:



Как можно видеть, станция проходила наиболее сильно с 10 утра до 10 вечера, ночью сигнал сильно ослаб, а к утру на пару часов совсем исчез. Также интересен медленный дрифт частоты несущей в пределах 10Гц, но относится ли он к передатчику или к приемнику, мне неизвестно. Интересны также периодические моменты «раздвоения» несущей, возможно это какая-то турбулентность слоев ионосферы, отражение которых вызывает допплеровский сдвиг частоты. Возможно, это связано с текущей активностью Солнца и солнечного ветра. Теоретически, большое количество работающих АМ-станций представляет собой некое подобие «пассивного радара», что позволяет, зная частоты и местоположение станций, фактически «бесплатно», без затрат на излучение, анализировать состояние ионосферы. Проводились ли такие работы, мне впрочем, неизвестно.

Кстати, как бонус отображения спектра — по нему можно видеть моменты, когда станция проходила наиболее сильно. В связи с этим, у любителей DX может возникнуть вопрос, как же прослушать саму передачу. Для этого лучше записывать сигнал с большей частотой дискретизации, например 16 КГц, а центр несущей расположить на частоте 4 КГц. Далее, сдвинуть спектр в ноль несложно в GNU Radio с помощью такого графа соединений:



При запуске блок-схемы в GNU Radio будет создан сконвертированный WAV-файл, прослушать который можно любым плеером.

Заключение


Как можно видеть, обработка сигналов дальних станций с помощью спектрального анализа может быть довольно интересной. Любопытно и то, что метод не требует дорогостоящего оборудования. Основная сложность лишь в длительности записи, но наверное, чтобы не держать ПК постоянно включенным, хватит и ресурсов Raspberry Pi 4.

Общие закономерности распространения радиоволн разумеется, давно известны, но увидеть это «вживую» гораздо интереснее, чем просто прочитать в учебнике. Такой способ анализа практически не используется среди радиолюбителей (мне известна только одна ветка на сайте radioscanner, где похожим методом искали моменты прохождения сигнала дальних станций, но геофизические аспекты участников того форума не интересовали, да и исходники там не выложены), так что наверно здесь еще могут найтись любопытные закономерности.

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