Multidimensional DSP with GPU acceleration

Digital signal processing (DSP) is a ubiquitous methodology in scientific and engineering computations. However, practically, DSP problems are often not only 1-D. For instance, image data are 2-D signals and radar signals are 3-D signals. While the number of dimension increases, the time and/or storage complexity of processing digital signals grow dramatically. Therefore, solving DSP problems in real-time is extremely difficult in reality.

Modern general purpose graphics processing units (GPGPUs) are considered as having an excellent throughput on vector operations and numeric manipulations through a high degree of parallel computations. While processing digital signals, particularly multidimensional signals, often involves a series of vector operations on massive amount of independent data samples, GPGPUs are now widely employed to accelerate multidimensional DSP, such as image processing, video codecs, radar signal analysis, sonar signal processing, and ultrasound scanning. Conceptually, using GPGPU devices to perform multidimensional DSP is able to dramatically reduce the computation complexity compared with central processing units (CPUs), digital signal processors (DSPs), or other FPGA accelerators.

Motivation

Processing multidimensional signals is a common problem in scientific researches and/or engineering computations. Typically, a DSP problem's computation complexity grows exponentially while the number of dimension increases. Notwithstanding, with a high degree of time and storage complexity, it is extremely difficult to process multidimensional signals in real-time. Although there are many fast algorithms (e.g. FFT) have been proposed in 1-D DSP problems, they are still not efficient enough to be adapted in high dimensional DSP problems. Therefore, it is still hard to obtain the computation results with digital signal processors (DSPs). Hence, a better solution of software algorithm or hardware architecture to accelerate multidimensional DSP computations is strongly required.

Existing Approaches

Practically, to accelerate multidimensional DSP, some common approaches have been proposed and developed in the past decades.

Lower Sampling Rate

A makeshift to achieve real-time requirement in multidimensional DSP applications is using a lower sampling rate can efficiently reduce the number of samples to be processed at one time and thereby decreasing the total processing time. However, this can lead to the aliasing problem in the sampling theorem and make a poor quality of outputs. In some applications, such as military radars and medical images, we are eager to have highly precise and accurate results. In such cases, using a lower sampling rate to reduce amount of computation in multidimensional DSP domain is not always allowable.

Digital Signal Processors (DSPs)

Digital signal processors are designed specifically to process vector operations. They have been widely used in DSP computations for decades. However, most digital signal processors are only capable of manipulating couple operations in parallel. This kind of designs is sufficient to accelerate audio processing (1-D signals) and image processing (2-D signals). However, with a large amount of data samples in multidimensional signals, this is still not powerful enough to retrieve computation results in real-time.

Supercomputer Assistance

In order to accelerate multidimensional DSP computations, using dedicated supercomputers or cluster computers is required in some circumstances, e.g., weather forecasting and military radars. Nevertheless, using supercomputers designated to simply perform DSP operations takes considerable money cost and energy consumption. It is not practical and suitable for all multidimensional DSP applications.

GPU Acceleration

GPUs are originally devised to accelerate image processing and video stream rendering. Moreover, since modern GPUs' have good ability to perform numeric computations in parallel with a relatively low cost and better energy efficiency, GPUs are becoming a popular alternative to replace supercomputers performing multidimensional DSP.

GPGPU Computations

GPGPU/SIMD computation model.

Modern GPU designs are mainly based on SIMD computation paradigm. This type of GPU devices is so-called general-purpose GPUs (GPGPUs).

GPGPUs are able to perform an operation on multiple independent data concurrently with their vector or SIMD functional units. A modern GPGPU can spawn thousands of concurrent threads and process all threads in a batch manner. With this nature, GPGPUs can be employed as DSP accelerators easily while many DSP problems can be solved by divide-and-conquer algorithms. A large scale and complex DSP problem can be divided into bunch of small numeric problems and be processed altogether at one time so that the overall time complexity can be reduced significantly. For example, multiplying two M × M matrices can be processed by M × M concurrent threads on a GPGPU device without any output data dependency. Therefore, theoretically, by means of GPGPU acceleration, we can gain up to M × M speedup compared with a traditional CPU or digital signal processor.

GPU Programming Languages

Currently, there are several existing programming languages or interfaces which support GPGPU programming.

CUDA

CUDA is the standard interface to program NVIDIA GPUs. NVIDIA also provides many CUDA libraries to support DSP acceleration on NVIDIA GPU devices.

OpenCL

OpenCL is an industrial standard which was originally proposed by Apple Inc. and is maintained and developed by Khronos Group now. OpenCL provides C++ like APIs for programming different devices universally, including GPGPUs.

OpenCL program execution flow

C++ Amp

C++ Amp is a programming model proposed by Microsoft. C++ Amp is a C++ based library designed for programming SIMD processors

OpenAcc

OpenAcc is a programming standard for parallel computing developed by Cray, CAPS, NVIDIA and PGI. OpenAcc targets to program for CPU and GPU heterogeneous systems with C, C++, and Fortran extensions.

Examples of GPU Programming for Multidimensional DSP

m × m Matrix Multiplication

Suppose A and B are two m × m matrices and we would like to compute C = A × B.

$\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}}$

$\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}$

To compute each element in C takes m multiplications and (m – 1) additions. Therefore, with a CPU implementation, the time complexity to achieve this computation is Θ(n³) in the following C example. However, we have known that elements in C are independent to each other. Hence, the computation can be fully parallelized by SIMD processors, such as GPGPU devices. With a GPGPU implementation, the time complexity significantly reduces to Θ(n) by unrolling the for-loop showing in the following OpenCL example.

// MxM matrix multiplication in C
void matrixMul(
    float *A,  // input matrix A
    float *B,  // input matrix B
    float *C,  // output matrix C
    int size)  // size of the matrices
{
    // N x N x N iterations
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            int id = row * size + col;
            flot sum = 0.0;

            for (int m = 0; m < size; m++) {
                sum += (A[row * size + m] * B[m * size + col]);
            }

            C[id] = sum;
        }
    }
}

// MxM matrix multiplication in OpenCL
__kernel void matrixMul(
    __global float *A,  // input matrix A
    __global float *B,  // input matrix B
    __global float *C,  // output matrix C
    __global int size)  // size of the matrices
{
    size_t id = get_global_id(0);   // each thread works on an element
    size_t row = id / size;
    size_t col = id % size;
    float sum = 0.0;

    // N iterations
    for (int m = 0; m < size; m++) {
        sum += (A[row * size + m] * B[m * size + col]);
    }

    C[id] = sum;
}

Multidimensional Convolution (M-D Convolution)

Convolution is a frequently used operation in DSP. To compute 2-D convolution of two m × m signals, it requires m² multiplications and m × (m – 1) additions for an output element. That is, the overall time complexity is Θ(n⁴) for the entire output signal. As the following OpenCL example shows, with GPGPU acceleration, the total computation time effectively decreases to Θ(n²) since all output elements are data independent.

$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})$

// 2-D convolution implementation in OpenCL
__kernel void convolution(
    __global float *x,  // input signal x
    __global float *h,  // filter h
    __global float *y,  // output signal y
    __global int size)  // size of ROS of the input signal and filter
{
    size_t id = get_global_id(0);   // each thread works on an element
    size_t row = size + size - 1;   // number of rows of the output signal
    size_t col = size + size - 1;   // number of columns of the output signal
    size_t n1 = id / row;
    size_t n2 = id % col;
    float sum = 0.0;

    // N x N iterations
    for (int k1 = 0; k1 < size; k1++) {
        for (int k2 = 0; k2 < size; k2++) {
            sum += (x[k1 * row + k2] * h[(n1 * row - k1) + (n2 - k2)]);
        }
    }

    C[id] = sum;
}

Note that, although the example demonstrated above is a 2-D convolution, a similar approach can be adopted for a higher dimension system. Overall, for a s-D convolution, a GPGPU implementation has time complexity Θ(n^s), whereas a CPU implementation has time complexity Θ(n^2s).

$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})$

Multidimensional Discrete Time Fourier Transform (M-D DTFT)

In addition to convolution, discrete time Fourier transform (DTFT) is another technique which is often used in system analysis.

$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})}$

Practically, to implement M-D DTFT, if assuming the system is separable, we can perform M times 1-D DFTF and matrix transpose with respect to each dimension. With a 1-D FFT operation, GPGPU can conceptually reduce the complexity from Θ(n²) to Θ(n) by the following example of OpenCL implementation. That is, to perform M-D DTFT, the complexity of GPGPU can be achieved by Θ(n²). While some GPGPUs' are also equipped hardware FFT accelerator internally, this implementation might be also optimized by invoking the FFT APIs or libraries provided by GPU manufactures.

// DTFT in OpenCL
__kernel void convolution(
    __global float *x_re,
    __global float *x_im,
    __global float *X_re,
    __global float *X_im,
    __global int size)
{
    size_t id = get_global_id(0);   // each thread works on an element
    X_re[id] = 0.0;
    X_im[id] = 0.0;

    for (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));
    }
}

Real Applications

Digital Filter Design

Designing a digital filter in multidimensional area is a big challenge, especially IIR filters. Typically it relies on computers to solve difference equations and obtain a set of approximated solutions. While GPGPU computation is becoming popular, several adaptive algorithms have been proposed to design multidimensional FIR and/or IIR filters by means of GPGPUs.

Radar Signal Reconstruction and Analysis

Radar systems usually require to reconstruct a numerous amount of 3-D or 4-D data samples in real-time. Traditionally, particularly in military, this needs supercomputers' support. Nowadays, GPGPUs are also employed to replace supercomputers to process radar signals. For example, to process synthetic aperture radar (SAR) signals, it usually involves multidimensional FFT computations. GPGPUs can be used to rapidly perform FFT and/or iFFT in this kind of applications.

Self-Driving Car

Many self-driving cars applies 3-D image recognition techniques to auto control the vehicles. Clearly, to accommodate the fast changing exterior environment, the recognition and decision processes must be done in real-time. GPGPUs are excellent devices to achieve the goal.

Medical Image Processing

In order to have accurate diagnosis, 2-D or 3-D medical signals, such as ultrasound, X-ray, MRI, and CT, often requires very high sampling rate and image resolutions to reconstruct images. By applying GPGPUs' superior computation power, it has been proved that we can acquire better quality on medical images