Наверное, каждый, кто хоть немного интересуется радио и связью, знает, что с помощью SDR-приемника можно одновременно принимать и обрабатывать широкий диапазон радиочастотного спектра. Отображение водопада в таких программах как HDSDR или SDR # неудивительно. Я покажу, как построить псевдо-3D спектр частотного диапазона FM, используя RTL-SDR, GNU Radio и около 100 строк кода Python.
Также возьмем более мощный приемник и рассмотрим весь FM диапазон 88–108 МГц.
Технически задача довольно простая. Приемник SDR оцифровывает входящий радиосигнал с помощью достаточно быстрого АЦП. На выходе мы получаем широкополосный IQ-сигнал в виде массива чисел, поступающий от АЦП с полосой пропускания, соответствующей частоте дискретизации АЦП. Частота АЦП определяет максимальную полосу пропускания, которую можно использовать. Это тот же процесс, что и в звуковой карте ПК, только не 22.050, а 2.000.000 или даже 10.000.000 отсчетов в секунду. Чтобы отобразить радиоспектр на экране, мы должны выполнить быстрое преобразование Фурье для массива данных, оно преобразует данные из так называемой «временной области» в «частотную область». Затем показываем данные на экране, и проблема решается. Я также постараюсь использовать минимум кода, чтобы GNU Radio могло помочь нам с обработкой данных.
Для тестов нам сначала понадобится приемник RTL-SDR, цена которого около 35 долларов. Это позволяет нам принимать радиосигналы в диапазоне частот 70..1700 МГц, используя полосу пропускания до 2 МГц:
Если кто-то хочет провести тесты с использованием RTL-SDR, рекомендуется получить приемник как на фото. Доступны более дешевые клоны, но они худшего качества.
Что ж, приступим.
Обработка радио GNU
Во-первых, нам нужно получить и обработать данные от приемника. График подключения GNU Radio показан на рисунке:
Как мы видим, мы получаем данные от приемника, конвертируем непрерывный поток данных в набор «векторов» размером 1024 значения, выполняем БПФ для этих векторов, преобразуем значения из комплексных в реальные и, наконец, отправляем данные через UDP. Конечно, все это можно сделать на чистом Python с использованием библиотек SoapySDR и numpy, но объем кода будет несколько больше.
Блок QT GUI Frequency Sink нужен только для «отладки», с помощью этого блока мы можем убедиться, что радиостанции видны, и при необходимости отрегулировать усиление приемника. Во время работы приложения картинка должна выглядеть примерно так:
Если все работает, блок Frequency Sink можно отключить, также в настройках проекта GNU Radio мы можем дополнительно указать режим «No GUI» и не тратить ресурсы пользовательского интерфейса. В принципе, эта программа может работать как служба без какого-либо пользовательского интерфейса.
Рендеринг
Поскольку мы передаем данные через UDP, мы можем получать их с любого клиента и даже на другом ПК. Буду использовать Python, для прототипа этого вполне достаточно.
Во-первых, давайте получим данные UDP:
fft_size = 1024 udp_data = None UDP_IP = "127.0.0.1" UDP_PORT = 40868sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) # UDP sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1) sock.bind((UDP_IP, UDP_PORT)) sock.settimeout(0.5)try: data, addr = sock.recvfrom(fft_size * 4) if len(data) == 4096: udp_data = np.frombuffer(data, dtype=np.float32) return True except socket.timeout: pass
Поскольку мы работаем с графикой, удобно использовать библиотеку pygame. Нарисовать трехмерный спектр просто: мы сохраняем данные в массиве и рисуем линии сверху вниз, от старых к новым.
fft_size = 1024 depth = 255 fft_data = np.zeros([depth, fft_size])def draw_image(screen, font): x_left, x_right, y_bottom = 0, img_size[0], img_size[1] - 5 # Draw spectrum in pseudo-3d for d in reversed(range(depth)): for x in range(fft_size - 1): d_x1, d_x2, d_y1, d_y2 = x + d, x + d + 1, y_bottom - int(y_ampl*fft_data[d][x]) - y_shift - d, y_bottom - int(y_ampl*fft_data[d][x+1]) - y_shift - d if d_y1 > y_bottom - 34: d_y1 = y_bottom - 34 if d_y2 > y_bottom - 34: d_y2 = y_bottom - 34 dim = 1 - 0.8*(d/depth) color = int(dim*data_2_color(fft_data[d][x])) pygame.draw.line(screen, (color//2,color,0) if d > 0 else (0, 250, 0), (d_x1, d_y1), (d_x2, d_y2), (2 if d == 0 else 1))
Мы также можем отображать частоты и названия станций на экране. Алгоритм преобразования Фурье дает на выходе 1024 точки, соответствующие полосе пропускания приемника. Мы знаем центральную частоту, поэтому вычислить положение пикселя можно по формуле для начальной школы.
stations = [("101.8 FM", 101.8), ("Rock FM", 102.4), ...]
for st_name, freq in stations:
x_pos = fft_size*(freq - center_freq)*1000000//sample_rate
textsurface = font.render(st_name, False, (255, 255, 0))
screen.blit(textsurface, (img_size[0]//2 + x_pos - textsurface.get_width()//2, y_bottom - 22))Собственно все, мы можем запускать обе программы одновременно, а на экране мы получаем панораму, показывающую в настоящее время работающие FM-станции в реальном времени:
Легко видеть, что разные станции имеют разные уровни сигнала, и мы даже можем различать моно и стереовещание.
Ну а теперь покажу обещанную панораму всего FM-диапазона. Для этого мы должны отложить RTL-SDR и использовать более качественное радио. Например, вот так:
Я использую SDR Ettus Research профессионального уровня, но с точки зрения кода все то же самое, просто нужно заменить один блок на другой в GNU Radio. И так выглядит спектр с полосой приема 24 МГц:
Интересно наблюдать за разнообразием силы сигнала от разных FM-станций.
Конечно, можно принимать не только FM-станции, но и любые другие в пределах рабочих частот SDR. Например, так выглядит воздушная повязка:
Мы можем видеть некоторые постоянно работающие частоты (вероятно, метеослужбу ATIS) и прерывистую радиосвязь между землей и пилотами. А вот так выглядит спектр GSM-диапазона (GSM-сигнал шире 24 МГц и не умещается полностью):
Вывод
Как видим, изучение радиоспектра довольно увлекательно, особенно в 3D. Конечно, здесь не было цели делать «еще один анализатор спектра», это просто прототип, сделанный для развлечения. Увы, рендеринг идет медленно, Python - не лучший выбор для отображения нескольких тысяч примитивов на экране. Алгоритм раскраски линий также может быть улучшен.
Как всегда, желаю всем читателям удачных экспериментов.
Полный исходный код рендеринга прилагается ниже.
import numpy as np
from matplotlib import pyplot as plt
from PIL import Image, ImageDraw
import sys
import pygame
from pygame.locals import *
from threading import Thread
import io
import cv2
import time
import socket
# FFT
receiver_name = "RTL-SDR"
center_freq = 102.5
sample_rate = 1800000
stations = [("101.8", 101.8), ("102.1", 102.1), ("102.4", 102.4), ("102.7", 102.7), ("103.0", 103.0), ("103.2", 103.2)]
# Load data from UDP
UDP_IP = "127.0.0.1"
UDP_PORT = 40868
udp_data = None
sock = None
# Panorama history
fft_size = 1024
depth = 255
fft_data = np.zeros([depth, fft_size])
# Canvas and draw
img_size = (fft_size, fft_size*9//16)
y_ampl = 90
color_ampl = 70
y_shift = 250
def udp_prepare():
global sock
sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) # UDP
sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
sock.bind((UDP_IP, UDP_PORT))
sock.settimeout(0.5)
def udp_getdata():
global sock, udp_data
try:
data, addr = sock.recvfrom(fft_size * 4)
if len(data) == 4096:
udp_data = np.frombuffer(data, dtype=np.float32)
return True
except socket.timeout:
pass
return False
def clear_data():
for y in range(depth):
fft_data[y, :] = np.full((fft_size,), -1024)
def add_new_line():
global udp_data, fft_data
# Shift old data up
for y in reversed(range(depth - 1)):
fft_data[y + 1, :] = fft_data[y, :]
# Put new data at the bottom line
if udp_data is not None:
fft_data[0, :] = udp_data
def datasock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) # UDP
sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
sock.bind((UDP_IP, UDP_PORT))
sock.settimeout(0.5)color(data):
c = -data + 2 # TODO: detect noise floor of the spectrum
color = 150 - int(color_ampl * c)
if color < 20:
color = 20
if color > 150:
color = 150
return color
def draw_image(screen, font):
x_left, x_right, y_bottom = 0, img_size[0], img_size[1] - 5
# Draw spectrum in pseudo-3d
for d in reversed(range(depth)):
for x in range(fft_size - 1):
d_x1, d_x2, d_y1, d_y2 = x + d, x + d + 1, y_bottom - int(y_ampl*fft_data[d][x]) - y_shift - d, y_bottom - int(y_ampl*fft_data[d][x+1]) - y_shift - d
if d_y1 > y_bottom - 34: d_y1 = y_bottom - 34
if d_y2 > y_bottom - 34: d_y2 = y_bottom - 34
dim = 1 - 0.8*(d/depth)
color = int(dim*datasock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) # UDP
sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
sock.bind((UDP_IP, UDP_PORT))
sock.settimeout(0.5)color(fft_data[d][x]))
pygame.draw.line(screen, (color//2,color,0) if d > 0 else (0, 250, 0), (d_x1, d_y1), (d_x2, d_y2), (2 if d == 0 else 1))
# Bottom line
pygame.draw.line(screen, (0,100,0), (x_left, y_bottom - 30), (x_right, y_bottom - 30), 2)
# Station names
for st_name, freq in stations:
x_pos = fft_size*(freq - center_freq)*1000000//sample_rate
textsurface = font.render(st_name, False, (255, 255, 0))
screen.blit(textsurface, (img_size[0]//2 + x_pos - textsurface.get_width()//2, y_bottom - 22))
text_mhz = font.render("MHz", False, (255, 255, 0))
screen.blit(text_mhz, (img_size[0] - 5 - text_mhz.get_width(), y_bottom - 22))
if __name__ == "__main__":
# UI init
screen = pygame.display.set_mode(img_size)
pygame.display.set_caption(receiver_name)
pygame.font.init()
font = pygame.font.SysFont('Arial Bold', 30)
# Subscribe to UDP
clear_data()
udp_prepare()
# Main loop
is_active = True
while is_active:
# Get new data
if udp_getdata():
add_new_line()
# Update screen
screen.fill((0, 0, 0))
draw_image(screen, font)
pygame.display.flip()
# Check sys events
for events in pygame.event.get():
if events.type == QUIT:
is_active = False
