Онлайн
Alvaros
.
- Регистрация
- 14.05.16
- Сообщения
- 21.452
- Реакции
- 101
- Репутация
- 204
Привет, Хабр.
Наверное все, хоть немного интересующиеся радиосвязью, знают что с помощью 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 это не лучший выбор. Также алгоритм раскрашивания линий оставляет желать лучшего, в идеале, он может быть более гибким.
Для желающих поэкспериментировать самостоятельно, исходный код:
sdr_render.py
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
sdr_receiver.grc
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
Как обычно, всем удачных экспериментов.
Наверное все, хоть немного интересующиеся радиосвязью, знают что с помощью 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 это не лучший выбор. Также алгоритм раскрашивания линий оставляет желать лучшего, в идеале, он может быть более гибким.
Для желающих поэкспериментировать самостоятельно, исходный код:
sdr_render.py
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
sdr_receiver.grc
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
Как обычно, всем удачных экспериментов.