Русский
Русский
English
Статистика
Реклама

Fft

Прием всего Bluetooth разом на SDR с CUDA? Легко

31.07.2020 16:11:33 | Автор: admin

В последнее время коллеги по "цеху" независимо друг от друга стали спрашивать меня: как получить c одного SDR-приемника одновременно все каналы Bluetooth? Полоса ведь позволяет, есть SDR с выходной полосой 80 МГц и более. Можно, конечно, сделать это на ПЛИС, но время такой разработки будет довольно большим. Мне давно было известно, что сделать такое на GPU довольно просто, но чтобы так!


Стандарт Bluetooth определяет физический уровень в двух версиях: Classic и Low Energy. Спецификация есть здесь. Документ ужасно большой, читать его целиком опасно для мозга. К счастью, большие компании, производящие измерительную технику, имеют средства для создания наглядных документов по теме. Tektronix и National Instruments, например. У меня совершенно нет шансов в конкуренции с ними по качеству представления материала. Интересующихся прошу изучать по ссылкам.


Все, что мне нужно знать о физическом уровне для создания многоканального фильтра, это шаг сетки частот и скорость модуляции. Они сведены в таблицу в одном из указанных документов:


image


Таким образом, нам нужно нарезать полосу 80 МГц на 79 фильтров с шагом настройки 1 МГц и, одновременно с этим, на 40 фильтров с шагом настройки 2 МГц. Частоты следования отсчетов с выходов фильтров должны быть 1 МГц и 2 МГц, соответственно.


Таким образом, нам необходимо две гребенки фильтров.


Для начала выберем параметры этих фильтров исходя из полос сигналов Bluetooth Classic и Bluetooth Low Energy. Нам нужны их импульсные характеристики, чтобы рассчитать нагрузку на вычислительное устройство фильтра. Здесь сразу стоит оговориться, что длины импульсных характеристик мы выбрали исходя из требований "быстрого" алгоритма фильтрации. Суть от этого не меняется. И число коэффициентов импульсной характеристики не должно быть слишком большим, чтобы фильтр был реализуем на вменяемой вычислительной аппаратуре.


Для фильтров с шагом 1 МГц выберем полосу пропускания ФНЧ (половина полосы пропускания полосового фильтра) 500 кГц, длину импульсной характеристики подгоним к 480 отводам. Для фильтров с шагом 2 МГц эти параметры выберем 1 МГц и 240 отводов, соответственно. Тип окна выбираем Кайзера. Рассчитаем импульсные характеристики в filterDesigner и выгрузим их в формате С-header:


Скриншоты из filterDesigner

image
image
image
image


Можно решать задачу лобовым способом: строить соответствующий по числу фильтров массив DDC (Digital Down Converter). Такой подход хорош для ПЛИС, где возможна экономия за счет уменьшения разрядности вычислителей первых ступеней. Также ПЛИС это наиболее энергетически-эффективный способ реализации. Но трудозатраты при этом способе наиболее высоки.


При выполнении гребенки фильтров на популярных нынче GPU появляется возможность реализации более хитрого алгоритма: полифазной гребенки фильтров на основе БПФ, который на CUDA доступен из библиотеки. В заграничной литературе алгоритм называется Polyphase or WOLA (Weight, Overlap and Add) FFT Filterbank. Лень к рисованию не дает мне возможности выполнить наглядное объяснение самостоятельно. В Сети есть много материалов по теме, особенно наглядный график выполнен здесь на странице 11 (огромное спасибо уважаемым авторам), вот он:


image


Я попробую пояснить схему обработки своими словами. Слабонервных прошу не читать.

Я попробую пояснить схему обработки своими словами в рамках моих методических возможностей. БПФ есть свертка входного сигнала со всем спектром комплексных ортогональных гармоник, укладывающихся на интервале импульсной характеристики. Импульсная характеристика фильтра, на которую умножается сигнал перед входом БПФ, модулируется этим спектров гармоник. Другими словами, огибающая импульсных характеристик результирующих фильтров гребенки выносится за скобки. Далее, спектр гармоник прореживается в некоторое число раз, ввиду расширения полосы пропускания фильтра относительно фильтра в прямоугольным окном. На рисунке мы видим прореживание на четыре. Другими словами, после расширения полосы окном Кайзера (с одновременным увеличением затухания в полосе задержания), нам уже нужны не все фильтры а только четвертая их часть. Остальные избыточны, их частотные характеристики перекрываются. Из четырех точек БПФ, идущих подряд, выбираем только нулевую, вычисление которой есть суммирование четырех
входных точек, взятых через время, равное четверти длительности исходного БПФ.


Железо выберем то, которое есть под рукой. Это плата ввода компании "Инструментальные Системы" FMC126P. О ней я уже писал в одной предыдущей статье. В разъем FMC платы вставлен субмодуль той же компании с трансивером AD9371 с полосой 100 МГц. Весь поток с трансивера может быть непрерывно передан в компьютер для обработки.


Выберем видео-карту с GPU GTX 1050. (Соврал, это она нас выбрала: это все что было под рукой, была выдрана из вычислителя для расчета антенн, но тем удивительней было увидеть рабочую гребенку). Перейдем к программной части.


К сожалению, из-за лицензий мы не можем опубликовать полный код. Можем показать только ядра GPU. Впрочем, остальной код и не особо интересен.


Вот ядро, которое выполняет перемножение сигнала на окно и сложение, и обертка для его вызова:


__global__ void cuComplexMultiplyWindowKernel(const cuComplex *data, const float *window, size_t windowSize, cuComplex *result) {    __shared__ cuComplex multiplicationResult[480];    multiplicationResult[threadIdx.x] = cuComplexMultiplyFloat(data[threadIdx.x + windowSize / 4 * blockIdx.x], window[threadIdx.x]);    __syncthreads();    cuComplex sum;    sum.x = sum.y = 0;    if (threadIdx.x < windowSize / 4) {        for(int i = 0; i < 4; i++) {            sum = cuComplexAdd(sum, multiplicationResult[threadIdx.x + i * windowSize / 4]);        }        result[threadIdx.x + windowSize / 4 * blockIdx.x] = sum;    }}cudaError_t cuComplexMultiplyWindow(const cuComplex *data, const float *window, size_t windowSize, cuComplex *result, size_t dataSize, cudaStream_t stream) {    size_t windowStep = windowSize / 4;    cuComplexMultiplyWindowKernel<<<dataSize / windowStep - 3, windowSize, 1024, stream>>>(data, window, windowSize, result);    return cudaGetLastError();}

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


Гребенки были проверены на выходном спектре каналов в реальном времени. На вход AD9371 подавался сигнал генератора сигналов 2450 МГц, селективность фильтров соответствовала расчетной.


image


В планах: адаптация софта на плату XRTX и реализация поиска пакетов, если это кому-нибудь окажется нужно или будет свободное время.


Всю работу по софту выполнил gaudima, слава ему!

Подробнее..

Изучаем распространение радиосигналов в ионосфере с помощью SDR

24.10.2020 14:21:16 | Автор: admin
Привет Хабр.

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



Я покажу как с помощью 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 npimport matplotlib.pyplot as pltimport wavefrom PIL import Imageimport sysimport timedef 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, где похожим методом искали моменты прохождения сигнала дальних станций, но геофизические аспекты участников того форума не интересовали, да и исходники там не выложены), так что наверно здесь еще могут найтись любопытные закономерности.

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

Категории

Последние комментарии

  • Имя: Макс
    24.08.2022 | 11:28
    Я разраб в IT компании, работаю на арбитражную команду. Мы работаем с приламы и сайтами, при работе замечаются постоянные баны и лаги. Пацаны посоветовали сервис по анализу исходного кода,https://app Подробнее..
  • Имя: 9055410337
    20.08.2022 | 17:41
    поможем пишите в телеграм Подробнее..
  • Имя: sabbat
    17.08.2022 | 20:42
    Охренеть.. это просто шикарная статья, феноменально круто. Большое спасибо за разбор! Надеюсь как-нибудь с тобой связаться для обсуждений чего-либо) Подробнее..
  • Имя: Мария
    09.08.2022 | 14:44
    Добрый день. Если обладаете такой информацией, то подскажите, пожалуйста, где можно найти много-много материала по Yggdrasil и его уязвимостях для написания диплома? Благодарю. Подробнее..
© 2006-2024, personeltest.ru