Procesamiento digital de señales multidimensional con aceleración por GPU

Una solución provisional para cumplir con los requisitos de tiempo real en aplicaciones DSP multidimensionales consiste en utilizar una frecuencia de muestreo más baja, lo que permite reducir eficazmente el número de muestras a procesar simultáneamente y, por lo tanto, disminuir el tiempo total de procesamiento. Sin embargo, esto puede provocar el problema del aliasing debido al teorema de muestreo y resultados de baja calidad. En algunas aplicaciones, como radares militares e imágenes médicas, se requieren resultados altamente precisos y exactos. En estos casos, utilizar una frecuencia de muestreo más baja para reducir la cantidad de cálculos en el dominio DSP multidimensional no siempre es viable.

Los procesadores de señales digitales están diseñados específicamente para procesar operaciones vectoriales. Se han utilizado ampliamente en cálculos DSP durante décadas. Sin embargo, la mayoría de los procesadores de señales digitales solo son capaces de manipular unas pocas operaciones en paralelo. Este tipo de diseño es suficiente para acelerar el procesamiento de audio (señales unidimensionales) y el procesamiento de imágenes (señales bidimensionales). No obstante, con un gran número de muestras de datos de señales multidimensionales, aún no son lo suficientemente potentes como para obtener resultados de cálculo en tiempo real.

Para acelerar los cálculos DSP multidimensionales, en ciertas circunstancias, como en la predicción meteorológica y los radares militares, se requiere el uso de supercomputadoras o clústeres de computadoras . Sin embargo, el uso de supercomputadoras destinadas únicamente a realizar operaciones DSP implica un costo económico y un consumo energético considerables. Además, no resulta práctico ni adecuado para todas las aplicaciones DSP multidimensionales.

Las GPU se diseñaron originalmente para acelerar el procesamiento de imágenes y la renderización de secuencias de vídeo. Además, dado que las GPU modernas tienen una buena capacidad para realizar cálculos numéricos en paralelo con un coste relativamente bajo y una mayor eficiencia energética, se están convirtiendo en una alternativa popular para reemplazar a las supercomputadoras que realizan DSP multidimensional. ^{[ 1 ]}

CUDA es la interfaz estándar para programar las GPU de Nvidia . Nvidia también proporciona muchas bibliotecas CUDA para admitir la aceleración DSP en los dispositivos GPU de Nvidia. ^{[ 4 ]}

OpenCL es un estándar industrial propuesto originalmente por Apple Inc. y mantenido y desarrollado actualmente por el Grupo Khronos . ^{[ 5 ]} OpenCL proporciona API similares a C++ para programar diferentes dispositivos de forma universal, incluidas las GPGPU.

La siguiente figura ilustra el flujo de ejecución de un programa OpenCL en una GPU. La CPU primero detecta los dispositivos OpenCL (la GPU en este caso) y luego invoca un compilador justo a tiempo para traducir el código fuente de OpenCL al binario de destino. A continuación, la CPU envía datos a la GPU para realizar los cálculos. Mientras la GPU procesa los datos, la CPU puede dedicarse a sus propias tareas.

C++ AMP es un modelo de programación propuesto por Microsoft . C++ AMP es una biblioteca basada en C++ diseñada para programar procesadores SIMD ^{[ 6 ].}

OpenACC es un estándar de programación para computación paralela desarrollado por Cray , CAPS, NVIDIA y PGI. ^{[ 7 ]} OpenACC está orientado a la programación para sistemas heterogéneos de CPU y GPU con extensiones de C , C++ y Fortran .

Supongamos que A y B son dos matrices m × m y que queremos calcular C = A × B.

${\displaystyle \mathbf {A} ={\begin{pmatrix}A_{11}&A_{12}&\cdots &A_{1m}\\A_{21}&A_{22}&\cdots &A_{2m}\\\vdots &\vdots &\ddots &\vdots \\A_{m1}&A_{m2}&\cdots &A_{mm}\\\end{pmatrix}},\quad \mathbf {B} ={\begin{pmatrix}B_{11}&B_{12}&\cdots &B_{1m}\\B_{21}&B_{22}&\cdots &B_{2m}\\\vdots &\vdots &\ddots &\vdots \\B_{m1}&B_{m2}&\cdots &B_{mm}\\\end{pmatrix}}}$

${\displaystyle \mathbf {C} =\mathbf {A} \times \mathbf {B} ={\begin{pmatrix}C_{11}&C_{12}&\cdots &C_{1m}\\C_{21}&C_{22}&\cdots &C_{2m}\\\vdots &\vdots &\ddots &\vdots \\C_{m1}&C_{m2}&\cdots &C_{mm}\\\end{pmatrix}},\quad C_{ij}=\sum _{k=1}^{m}A_{ik}B_{kj}}$

Para calcular cada elemento en C se requieren m multiplicaciones y ( m – 1 ) sumas. Por lo tanto, con una implementación en CPU, la complejidad temporal para lograr este cálculo es Θ(n³ ⁾ en el siguiente ejemplo de C. Sin embargo, sabemos que los elementos en C son independientes entre sí. Por lo tanto, el cálculo puede paralelizarse completamente mediante procesadores SIMD, como los dispositivos GPGPU. Con una implementación en GPGPU, la complejidad temporal se reduce significativamente a Θ(n) al desenrollar el bucle for como se muestra en el siguiente ejemplo de OpenCL .

// Multiplicación de matrices MxM en Cmatriz vacíaMul (float * A , // matriz de entrada Afloat * B , // matriz de entrada Bfloat * C , // matriz de salida Cint tamaño ) // tamaño de las matrices{// N x N x N iteracionespara ( int fila = 0 ; fila < tamaño ; fila ++ ) {para ( int col = 0 ; col < size ; col ++ ) {int id = fila * tamaño + columna ;suma flotante = 0.0 ;para ( int m = 0 ; m < tamaño ; m ++ ) {suma += ( A [ fila * tamaño + m ] * B [ m * tamaño + columna ]);}C [ id ] = suma ;}}}

// Multiplicación de matrices MxM en OpenCL__kernel void matrixMul (__global float * A , // matriz de entrada A__global float * B , // matriz de entrada B__global float * C , // matriz de salida C__global int size ) // tamaño de las matrices{size_t id = get_global_id ( 0 ); // cada hilo trabaja en un elementosize_t fila = id / tamaño ;tamaño_t col = id % tamaño ;suma flotante = 0.0 ;// N iteracionespara ( int m = 0 ; m < tamaño ; m ++ ) {suma += ( A [ fila * tamaño + m ] * B [ m * tamaño + columna ]);}C [ id ] = suma ;}

^{La convolución es una operación muy utilizada en el procesamiento digital de señales} (DSP). Para calcular la convolución bidimensional de dos señales m × m , se requieren m² multiplicaciones y m × ( m – 1 ) sumas para cada elemento de salida. Es decir, la complejidad temporal total es Θ(n⁴ ⁾ para toda la señal de salida . Como muestra el siguiente ejemplo de OpenCL, con la aceleración de GPGPU, el tiempo total de cálculo se reduce efectivamente a Θ(n² ⁾ ya que todos los elementos de salida son independientes de los datos.

Ecuación de convolución 2D:

${\displaystyle y(n_{1},n_{2})=x(n_{1},n_{2})**h(n_{1},n_{2})=\sum _{k_{1}=0}^{m-1}\sum _{k_{2}=0}^{m-1}x(k_{1},k_{2})h(n_{1}-k_{1},n_{2}-k_{2})}$

// Implementación de convolución 2D en OpenCL__kernel void convolución (__global float * x , // señal de entrada x__global float * h , // filtrar h__global float * y , // señal de salida y__global int size ) // tamaño de ROS de la señal de entrada y del filtro{size_t id = get_global_id ( 0 ); // cada hilo trabaja en un elementosize_t fila = tamaño + tamaño - 1 ; // número de filas de la señal de salidasize_t col = size + size - 1 ; // número de columnas de la señal de salidatamaño_t n1 = id / fila ;tamaño_t n2 = id % col ;suma flotante = 0.0 ;// N x N iteracionespara ( int k1 = 0 ; k1 < size ; k1 ++ ) {para ( int k2 = 0 ; k2 < size ; k2 ++ ) {suma += ( x [ k1 * fila + k2 ] * h [( n1 * fila - k1 ) + ( n2 - k2 )]);}}C [ id ] = suma ;}

Cabe señalar que, si bien el ejemplo mostrado anteriormente es una convolución 2D, se puede adoptar un enfoque similar para un sistema de mayor dimensión. En general, para una convolución sD, una implementación GPGPU tiene una complejidad temporal Θ(n ^s ) , mientras que una implementación CPU tiene una complejidad temporal Θ(n ^2s ) .

Ecuación de convolución MD:

${\displaystyle y(n_{1},n_{2},...,n_{s})=x(n_{1},n_{2},...,n_{s})**h(n_{1},n_{2},...,n_{s})=\sum _{k_{1}=0}^{m_{1}-1}\sum _{k_{2}=0}^{m_{2}-1}...\sum _{k_{s}=0}^{m_{s}-1}x(k_{1},k_{2},...,k_{s})h(n_{1}-k_{1},n_{2}-k_{2},...,n_{s}-k_{s})}$

Además de la convolución, la transformada discreta de Fourier (DTFT) es otra técnica que se utiliza con frecuencia en el análisis de sistemas.

${\displaystyle X(\Omega _{1},\Omega _{2},...,\Omega _{s})=\sum _{n_{1}=0}^{m_{1}-1}\sum _{n_{2}=0}^{m_{2}-1}...\sum _{n_{s}=0}^{m_{s}-1}x(n_{1},n_{2},...,n_{s})e^{-j(\Omega _{1}n_{1}+\Omega _{1}n_{1}+...+\Omega _{s}n_{s})}}$

En la práctica, para implementar una MD DTFT, podemos realizar M veces una DFTF 1-D y transponer la matriz con respecto a cada dimensión. Con una operación DTFT 1-D, GPGPU puede reducir conceptualmente la complejidad de Θ(n² ⁾ a Θ(n ) , como se ilustra en el siguiente ejemplo de implementación OpenCL . Es decir, una MD DTFT con la complejidad de GPGPU se puede calcular en una GPU con una complejidad de Θ(n² ⁾ . Si bien algunas GPGPU también están equipadas con aceleradores FFT de hardware internos, esta implementación también podría optimizarse invocando las API o bibliotecas FFT proporcionadas por el fabricante de la GPU. ^{[ 8 ]}

// DTFT en OpenCL__kernel void convolución (__global float * x_re ,__global float * x_im ,__global float * X_re ,__global float * X_im ,__tamaño entero global ){size_t id = get_global_id ( 0 ); // cada hilo trabaja en un elementoX_re [ id ] = 0.0 ;X_im [ id ] = 0.0 ;para ( int i = 0 ; i < size ; i ++ ) {X_re += ( x_re [ id ] * cos ( 2 * 3.1415 * id / size ) - x_im [ id ] * sin ( 2 * 3.1415 * id / size ));X_im += ( x_re [ id ] * sin ( 2 * 3.1415 * id / size ) + x_im [ id ] * cos ( 2 * 3.1415 * id / size ));}}

El diseño de filtros digitales multidimensionales representa un gran desafío, especialmente para los filtros IIR . Generalmente, se recurre a computadoras para resolver ecuaciones en diferencias y obtener un conjunto de soluciones aproximadas. Si bien la computación en GPGPU se está popularizando, se han propuesto diversos algoritmos adaptativos para diseñar filtros FIR y/o IIR multidimensionales mediante GPGPU. ^{[ 9 ] [ 10 ] [ 11 ]}

Los sistemas de radar suelen necesitar reconstruir numerosas muestras de datos 3D o 4D en tiempo real. Tradicionalmente, sobre todo en el ámbito militar, esto requiere el apoyo de supercomputadoras. Actualmente, las GPGPU también se emplean para reemplazar a las supercomputadoras en el procesamiento de señales de radar. Por ejemplo, para procesar señales de radar de apertura sintética (SAR) , generalmente se requieren cálculos de FFT multidimensionales . ^{[ 12 ] [ 13 ] [ 14 ]} Las GPGPU pueden utilizarse para realizar rápidamente FFT y/o iFFT en este tipo de aplicaciones.

Muchos vehículos autónomos aplican técnicas de reconocimiento de imágenes 3D para el autocontrol. Evidentemente, para adaptarse al entorno exterior que cambia rápidamente, los procesos de reconocimiento y toma de decisiones deben realizarse en tiempo real. Las GPGPU son dispositivos excelentes para lograr este objetivo. ^{[ 15 ]}

Para obtener un diagnóstico preciso, las señales médicas bidimensionales o tridimensionales, como las obtenidas mediante ultrasonido , rayos X , resonancia magnética y tomografía computarizada , suelen requerir una alta frecuencia de muestreo y una elevada resolución de imagen para su reconstrucción. Se ha demostrado que, gracias a la potencia de cálculo superior de las GPGPU, se pueden obtener imágenes médicas de mejor calidad ^{[ 16 ] [ 17 ].}