Панорама-FM или как увидеть все радиостанции сразу с помощью SDR
Привет, Хабр.
Наверное все, хоть немного интересующиеся радиосвязью, знают что с помощью SDR-приемника возможно принимать и обрабатывать широкую полосу спектра радиодиапазона. Собственно, отображением спектра в таких программах как HDSDR или SDR# никого не удивить. Я покажу как построить псевдо-3D спектр принимаемых станций с помощью RTL-SDR, GNU Radio и примерно 100 строк кода на языке Python.
Также мы возьмем приемник посерьезнее, и посмотрим на весь FM спектр 88–108.
Продолжение под катом.
Технически, задача довольно проста. SDR-приемник оцифровывает входящий радиосигнал с помощью довольно быстрого АЦП. На выходе мы получаем широкополосный IQ-сигнал в виде массива чисел, приходящих с АЦП с частотой, соответствующей частоте дискретизации. Эта же частота определяет максимальную ширину полосы, принимаемую приемником. Все тоже самое, что и в звуковой карте, только в секунду мы имеем не 22050, а 2000000 или даже 10000000 значений. Для того, чтобы вывести спектр на экран, мы должны выполнить над массивом данных преобразование Фурье. Дальше остается вывести все это на экран, и задача решена. Остается лишь вопрос, как обойтись минимумом кода, для чего нам поможет GNU Radio.
Для тестов нам также понадобится RTL-SDR приемник, цена которого составляет порядка 30$. Он позволяет принимать сигналы в диапазоне частот примерно от 70 до 1700МГц и шириной полосы до 2МГц:
Если кто захочет повторить тесты с RTL-SDR самостоятельно, приемник нужен именно такой как на фото. Есть более дешевые клоны, но они хуже качеством.
Ну, а мы приступим.
GNU Radio
Граф связей в GNU Radio показан на рисунке:
Как можно видеть, мы получаем данные с приемника, затем преобразуем непрерывный поток в набор «векторов» размером 1024, над этими блоками выполняем FFT, преобразуем значения из комплексных в вещественные, и отправляем данные по UDP. Разумеется, все это можно было бы сделать и непосредственно в Python с помощью SoapySDR и numpy, но объем кода был бы несколько больше.
Блок QT GUI Frequency Sink нужен для «отладки», с ним можно убедиться что сигналы радиостанций действительно принимаются и при необходимости настроить усиление приемника. Во время работы программы картинка должна быть примерно такой:
Если все работает, GUI Frequency Sink можно отключить, а в настройках проекта указать «No GUI», чтобы не тратить зря ресурсы. В принципе, эту программу дальше можно запускать как сервис, без какого-либо UI.
Отрисовка
Поскольку мы передаем данные по UDP, то принимать их мы можем любым клиентом, и даже на другом ПК. Будем использовать Python, для прототипа этого вполне достаточно.
Сначала мы получаем данные по UDP:
fft_size = 1024
udp_data = None
UDP_IP = "127.0.0.1"
UDP_PORT = 40868
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)
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. Рисование 3D-спектра делается просто, мы храним данные в массиве, и выводим линии сверху вниз, от более старых к более новым.
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", 101.8), ("102.1", 102.1), ("102.4", 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))
Собственно и все, объединяем все это вместе, запускаем обе программы одновременно, и получаем на экране real-time панораму, показывающую работающие в данный момент FM-станции:
Легко можно видеть, что разные станции вещают с разным уровнем, или по ширине полосы отличить моно-вещание от стерео.
Ну, а теперь обещанная панорама всего FM-диапазона. Для этого нам придется отложить RTL-SDR и взять приемник посерьезнее, например такой:
Я использую профессиональную модель от Ettus Research, но с точки зрения кода все тоже самое, достаточно поменять в GNU Radio один блок на другой. А вот так это выглядит на спектре при ширине полосы приема 24МГц:
Разумеется, принимать можно не только станции FM-диапазона, но и любые другие, в пределах рабочих частот SDR-приемника. Для примера, вот так выглядит авиадиапазон:
Можно видеть постоянно активные частоты (вероятно метеослужба ATIS) и периодически возникающий радиообмен. А вот так выглядит фрагмент спектра 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 data_2_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*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))
# 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
Sun Jun 7 13:17:58 2020
options
author
window_size
category
[GRC Hier Blocks]
comment
description
_enabled
True
_coordinate
(8, 8)
_rotation
0
generate_options
no_gui
hier_block_src_path
.:
id
top_block
max_nouts
0
qt_qss_theme
realtime_scheduling
run_command
{python} -u {filename}
run_options
prompt
run
True
sizing_mode
fixed
thread_safe_setters
title
placement
(0,0)
variable
comment
_enabled
True
_coordinate
(168, 12)
_rotation
0
id
fft_size
value
1024
variable
comment
_enabled
True
_coordinate
(256, 12)
_rotation
0
id
samp_rate
value
1800000
blocks_complex_to_mag_squared
alias
comment
affinity
_enabled
True
_coordinate
(648, 108)
_rotation
0
id
blocks_complex_to_mag_squared_0
maxoutbuf
0
minoutbuf
0
vlen
fft_size
blocks_nlog10_ff
alias
comment
affinity
_enabled
True
_coordinate
(832, 100)
_rotation
0
id
blocks_nlog10_ff_0
maxoutbuf
0
minoutbuf
0
vlen
fft_size
k
0
n
1
blocks_stream_to_vector
alias
comment
affinity
_enabled
True
_coordinate
(256, 108)
_rotation
0
id
blocks_stream_to_vector_0
type
complex
maxoutbuf
0
minoutbuf
0
num_items
fft_size
vlen
1
blocks_udp_sink
alias
comment
affinity
ipaddr
127.0.0.1
port
40868
_enabled
True
_coordinate
(1000, 76)
_rotation
0
id
blocks_udp_sink_0
type
float
psize
fft_size*4
eof
True
vlen
fft_size
fft_vxx
alias
comment
affinity
_enabled
True
fft_size
fft_size
forward
True
_coordinate
(424, 76)
_rotation
0
id
fft_vxx_0
type
complex
maxoutbuf
0
minoutbuf
0
nthreads
1
shift
True
window
window.blackmanharris(fft_size)
qtgui_freq_sink_x
autoscale
False
average
1.0
axislabels
True
bw
samp_rate
alias
fc
102.5e6
comment
ctrlpanel
False
affinity
_enabled
0
fftsize
1024
_coordinate
(256, 188)
gui_hint
_rotation
0
grid
False
id
qtgui_freq_sink_x_0
legend
True
alpha1
1.0
color1
"blue"
label1
width1
1
alpha10
1.0
color10
"dark blue"
label10
width10
1
alpha2
1.0
color2
"red"
label2
width2
1
alpha3
1.0
color3
"green"
label3
width3
1
alpha4
1.0
color4
"black"
label4
width4
1
alpha5
1.0
color5
"cyan"
label5
width5
1
alpha6
1.0
color6
"magenta"
label6
width6
1
alpha7
1.0
color7
"yellow"
label7
width7
1
alpha8
1.0
color8
"dark red"
label8
width8
1
alpha9
1.0
color9
"dark green"
label9
width9
1
maxoutbuf
0
minoutbuf
0
name
""
nconnections
1
showports
True
freqhalf
True
tr_chan
0
tr_level
0.0
tr_mode
qtgui.TRIG_MODE_FREE
tr_tag
""
type
complex
update_time
0.10
wintype
firdes.WIN_BLACKMAN_hARRIS
label
Relative Gain
ymax
10
ymin
-140
units
dB
rtlsdr_source
alias
ant0
bb_gain0
20
bw0
0
dc_offset_mode0
0
corr0
0
freq0
102.5e6
gain_mode0
False
if_gain0
20
iq_balance_mode0
0
gain0
20
ant10
bb_gain10
20
bw10
0
dc_offset_mode10
0
corr10
0
freq10
100e6
gain_mode10
False
if_gain10
20
iq_balance_mode10
0
gain10
10
ant11
bb_gain11
20
bw11
0
dc_offset_mode11
0
corr11
0
freq11
100e6
gain_mode11
False
if_gain11
20
iq_balance_mode11
0
gain11
10
ant12
bb_gain12
20
bw12
0
dc_offset_mode12
0
corr12
0
freq12
100e6
gain_mode12
False
if_gain12
20
iq_balance_mode12
0
gain12
10
ant13
bb_gain13
20
bw13
0
dc_offset_mode13
0
corr13
0
freq13
100e6
gain_mode13
False
if_gain13
20
iq_balance_mode13
0
gain13
10
ant14
bb_gain14
20
bw14
0
dc_offset_mode14
0
corr14
0
freq14
100e6
gain_mode14
False
if_gain14
20
iq_balance_mode14
0
gain14
10
ant15
bb_gain15
20
bw15
0
dc_offset_mode15
0
corr15
0
freq15
100e6
gain_mode15
False
if_gain15
20
iq_balance_mode15
0
gain15
10
ant16
bb_gain16
20
bw16
0
dc_offset_mode16
0
corr16
0
freq16
100e6
gain_mode16
False
if_gain16
20
iq_balance_mode16
0
gain16
10
ant17
bb_gain17
20
bw17
0
dc_offset_mode17
0
corr17
0
freq17
100e6
gain_mode17
False
if_gain17
20
iq_balance_mode17
0
gain17
10
ant18
bb_gain18
20
bw18
0
dc_offset_mode18
0
corr18
0
freq18
100e6
gain_mode18
False
if_gain18
20
iq_balance_mode18
0
gain18
10
ant19
bb_gain19
20
bw19
0
dc_offset_mode19
0
corr19
0
freq19
100e6
gain_mode19
False
if_gain19
20
iq_balance_mode19
0
gain19
10
ant1
bb_gain1
20
bw1
0
dc_offset_mode1
0
corr1
0
freq1
100e6
gain_mode1
False
if_gain1
20
iq_balance_mode1
0
gain1
10
ant20
bb_gain20
20
bw20
0
dc_offset_mode20
0
corr20
0
freq20
100e6
gain_mode20
False
if_gain20
20
iq_balance_mode20
0
gain20
10
ant21
bb_gain21
20
bw21
0
dc_offset_mode21
0
corr21
0
freq21
100e6
gain_mode21
False
if_gain21
20
iq_balance_mode21
0
gain21
10
ant22
bb_gain22
20
bw22
0
dc_offset_mode22
0
corr22
0
freq22
100e6
gain_mode22
False
if_gain22
20
iq_balance_mode22
0
gain22
10
ant23
bb_gain23
20
bw23
0
dc_offset_mode23
0
corr23
0
freq23
100e6
gain_mode23
False
if_gain23
20
iq_balance_mode23
0
gain23
10
ant24
bb_gain24
20
bw24
0
dc_offset_mode24
0
corr24
0
freq24
100e6
gain_mode24
False
if_gain24
20
iq_balance_mode24
0
gain24
10
ant25
bb_gain25
20
bw25
0
dc_offset_mode25
0
corr25
0
freq25
100e6
gain_mode25
False
if_gain25
20
iq_balance_mode25
0
gain25
10
ant26
bb_gain26
20
bw26
0
dc_offset_mode26
0
corr26
0
freq26
100e6
gain_mode26
False
if_gain26
20
iq_balance_mode26
0
gain26
10
ant27
bb_gain27
20
bw27
0
dc_offset_mode27
0
corr27
0
freq27
100e6
gain_mode27
False
if_gain27
20
iq_balance_mode27
0
gain27
10
ant28
bb_gain28
20
bw28
0
dc_offset_mode28
0
corr28
0
freq28
100e6
gain_mode28
False
if_gain28
20
iq_balance_mode28
0
gain28
10
ant29
bb_gain29
20
bw29
0
dc_offset_mode29
0
corr29
0
freq29
100e6
gain_mode29
False
if_gain29
20
iq_balance_mode29
0
gain29
10
ant2
bb_gain2
20
bw2
0
dc_offset_mode2
0
corr2
0
freq2
100e6
gain_mode2
False
if_gain2
20
iq_balance_mode2
0
gain2
10
ant30
bb_gain30
20
bw30
0
dc_offset_mode30
0
corr30
0
freq30
100e6
gain_mode30
False
if_gain30
20
iq_balance_mode30
0
gain30
10
ant31
bb_gain31
20
bw31
0
dc_offset_mode31
0
corr31
0
freq31
100e6
gain_mode31
False
if_gain31
20
iq_balance_mode31
0
gain31
10
ant3
bb_gain3
20
bw3
0
dc_offset_mode3
0
corr3
0
freq3
100e6
gain_mode3
False
if_gain3
20
iq_balance_mode3
0
gain3
10
ant4
bb_gain4
20
bw4
0
dc_offset_mode4
0
corr4
0
freq4
100e6
gain_mode4
False
if_gain4
20
iq_balance_mode4
0
gain4
10
ant5
bb_gain5
20
bw5
0
dc_offset_mode5
0
corr5
0
freq5
100e6
gain_mode5
False
if_gain5
20
iq_balance_mode5
0
gain5
10
ant6
bb_gain6
20
bw6
0
dc_offset_mode6
0
corr6
0
freq6
100e6
gain_mode6
False
if_gain6
20
iq_balance_mode6
0
gain6
10
ant7
bb_gain7
20
bw7
0
dc_offset_mode7
0
corr7
0
freq7
100e6
gain_mode7
False
if_gain7
20
iq_balance_mode7
0
gain7
10
ant8
bb_gain8
20
bw8
0
dc_offset_mode8
0
corr8
0
freq8
100e6
gain_mode8
False
if_gain8
20
iq_balance_mode8
0
gain8
10
ant9
bb_gain9
20
bw9
0
dc_offset_mode9
0
corr9
0
freq9
100e6
gain_mode9
False
if_gain9
20
iq_balance_mode9
0
gain9
10
comment
affinity
args
_enabled
1
_coordinate
(16, 96)
_rotation
0
id
rtlsdr_source_0
maxoutbuf
0
clock_source0
time_source0
clock_source1
time_source1
clock_source2
time_source2
clock_source3
time_source3
clock_source4
time_source4
clock_source5
time_source5
clock_source6
time_source6
clock_source7
time_source7
minoutbuf
0
nchan
1
num_mboards
1
type
fc32
sample_rate
samp_rate
sync
blocks_complex_to_mag_squared_0
blocks_nlog10_ff_0
0
0
blocks_nlog10_ff_0
blocks_udp_sink_0
0
0
blocks_stream_to_vector_0
fft_vxx_0
0
0
fft_vxx_0
blocks_complex_to_mag_squared_0
0
0
rtlsdr_source_0
blocks_stream_to_vector_0
0
0
rtlsdr_source_0
qtgui_freq_sink_x_0
0
0
Как обычно, всем удачных экспериментов.
