🥇Odbiór całego Bluetooth jednocześnie na SDR z CUDA? Łatwo

Ostatnio koledzy z branży niezależnie zaczęli mnie pytać, jak uzyskać z jednego odbiornika SDR jednocześnie wszystkie kanały Bluetooth? Pasmo pozwala, jest SDR z pasmem wyjściowym 80 MHz i więcej. Oczywiście można to zrobić na FPGA, ale czas realizacji takiego projektu będzie dość długi. Od dłuższego czasu wiadomo mi, że zrealizowanie tego na GPU jest dość proste, ale aby tak!

Standard Bluetooth definiuje poziom fizyczny w dwóch wersjach: Classic i Low Energy. Specyfikacja jest tutaj. Dokument jest strasznie obszerny, czytanie go w całości może być niebezpieczne dla mózgu. Na szczęście duże firmy produkujące sprzęt pomiarowy mają środki do tworzenia przystępnych dokumentów na ten temat. Tektronix i National Instruments, na przykład. Nie mam żadnych szans w konkurencji z nimi pod względem jakości prezentacji materiału. Osoby zainteresowane proszę o zapoznanie się z linkami.

Wszystko, co muszę wiedzieć o poziomie fizycznym do stworzenia filtrów wielokanałowych, to krok siatki częstotliwości i prędkość modulacji. Zostały one przedstawione w tabeli w jednym z wymienionych dokumentów:

W związku z tym musimy podzielić pasmo 80 MHz na 79 filtrów z krokiem ustawienia 1 MHz i jednocześnie na 40 filtrów z krokiem ustawienia 2 MHz. Częstotliwości próbkowania z wyjść filtrów powinny wynosić odpowiednio 1 MHz i 2 MHz.

W związku z tym potrzebujemy dwóch grzebieni filtrów.

Na początku wybierzemy parametry tych filtrów na podstawie pasm sygnałów Bluetooth Classic i Bluetooth Low Energy. Potrzebujemy ich charakterystyk impulsowych, aby obliczyć obciążenie obliczeniowe dla filtra. Tu od razu trzeba zaznaczyć, że długości charakterystyk impulsowych wybraliśmy na podstawie wymagań "szybkiego" algorytmu filtracji. Istota się od tego nie zmienia. A liczba współczynników charakterystyki impulsowej nie powinna być zbyt duża, aby filtr był realizowalny na rozsądnej aparaturze obliczeniowej.

Dla filtrów z krokiem 1 MHz wybierzemy pasmo przenoszenia LPF (połowa pasma przenoszenia filtru pasmowego) 500 kHz, długość charakterystyki impulsowej dostosujemy do 480 odprowadzeń. Dla filtrów z krokiem 2 MHz wybierzemy te parametry na 1 MHz i 240 odprowadzeń, odpowiednio. Typ okna wybieramy Kaisera. Obliczymy charakterystyki impulsowe w filterDesigner i wyeksportujemy je w formacie nagłówka C:

Zrzuty ekranu z filterDesigner

Można rozwiązać zadanie w sposób bezpośredni: zbudować odpowiednią tablicę DDC zgodnie z liczbą filtrów (Digital Down Converter). Takie podejście jest dobre dla FPGA, gdzie istnieje możliwość oszczędności przez zmniejszenie precyzji obliczeń pierwszych etapów. FPGA to także najbardziej energooszczędny sposób realizacji. Jednak nakład pracy w tym przypadku jest najwyższy.

Podczas realizacji filtrów grzebieniowych na popularnych obecnie GPU istnieje możliwość zaimplementowania bardziej zaawansowanego algorytmu: wielofazowej grzebieniowej banku filtrów opartych na FFT, który jest dostępny w CUDA z biblioteki. W literaturze zagranicznej algorytm nazywa się Polyphase lub WOLA (Weight, Overlap and Add) FFT Filterbank. Oparcie o rysunki nie pozwala mi na samodzielne wykonanie wizualnego wyjaśnienia. W Internecie jest wiele materiałów na ten temat, szczególnie wyraźny wykres można znaleźć tutaj na stronie 11 (ogromne dzięki szanownym autorom), oto on:

Postaram się wyjaśnić schemat przetwarzania swoimi słowami. Osoby wrażliwe proszę o nieczytanie.

Postaram się wyjaśnić schemat przetwarzania swoimi słowami w ramach moich metodologicznych możliwości. FFT jest splotem sygnału wejściowego ze wszystkimi spektrami ortogonalnych harmonicznych, które układają się na interwał impulsowej charakterystyki. Impulsowa charakterystyka filtra, na którą sygnał jest mnożony przed wejściem FFT, jest modulowana tymi spektrami harmonicznych. Innymi słowy, obwiednia impulsowych charakterystyk wynikowych filtrów grzebieniowych jest wyciągana poza nawias. Następnie, spektrum harmonicznych jest przerzedzane w pewnej liczbie razy, ze względu na rozszerzenie pasma przepustowego filtra w porównaniu do filtra w oknie prostokątnym. Na rysunku widzimy przerzedzenie o cztery. Innymi słowy, po rozszerzeniu pasma w oknie Keizera (z równoczesnym zwiększeniem tłumienia w paśmie zatrzymania), potrzebujemy już nie wszystkich filtrów, a jedynie ich czwartą część. Pozostałe są zbędne, ich charakterystyki częstotliwościowe nakładają się. Z czterech punktów FFT, które idą kolejno, wybieramy tylko zerowy, którego obliczenie polega na sumowaniu czterech
punktów wejściowych, wziętych w czasie równym jednej czwartej długości podstawowego FFT.

Sprzęt wybierzemy ten, który jest dostępny. Jest to karta wejściowa firmy "Instrumentalne Systemy" FMC126P. O niej już pisałem w jednym z poprzednich artykuł. W gnieździe FMC płyty włożony jest submoduł tej samej firmy z transceiverem AD9371 o przepustowości 100 MHz. Cały strumień z transceivera może być ciągle przesyłany do komputera w celu przetwarzania.

Wybierzemy kartę graficzną z GPU GTX 1050. (Nieprawda, to ona nas wybrała: to wszystko, co mieliśmy pod ręką, została wydarta z komputera obliczeniowego do obliczeń antenowych, ale tym bardziej zaskakujące było zobaczenie działającej grzebień). Przejdźmy do części programowej.

Niestety, z powodu licencji nie możemy opublikować pełnego kodu. Możemy pokazać tylko jądra GPU. A zresztą, reszta kodu nie jest zbyt interesująca.

Oto jądro, które wykonuje mnożenie sygnału przez okno i sumowanie, a także wrapper do jego wywołania:

__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<<>>(data, window, windowSize, result);
    return cudaGetLastError();
}

Kod przetwarzania sygnału, który wywołuje to jądro, ściśle powtarza schemat algorytmu przedstawionego na powyższym rysunku, dlatego nie widzę sensu przytaczać go tutaj.

Grzebienie były sprawdzane na wyjściowym widmie kanałów w czasie rzeczywistym. Do wejścia AD9371 podawany był sygnał generatora o częstotliwości 2450 MHz, selektywność filtrów odpowiadała obliczonej.