GPUアクセラレーションを備えた多次元DSP

多次元DSPアプリケーションでリアルタイム要件を満たすための暫定的な方法として、サンプリングレートを低くすることが挙げられます。これにより、一度に処理するサンプル数を効率的に削減し、総処理時間を短縮できます。しかし、サンプリング定理によるエイリアシング問題や出力品質の低下につながる可能性があります。軍事レーダーや医療画像などの一部のアプリケーションでは、非常に高精度で正確な結果が求められます。このような場合、多次元DSP領域での計算量を削減するためにサンプリングレートを低くすることは、必ずしも許容されるとは限りません。

デジタル信号プロセッサは、ベクトル演算を処理するために特化して設計されています。これらは数十年にわたりDSP演算で広く使用されてきました。しかし、ほとんどのデジタル信号プロセッサは、並列処理できる演算数が限られています。このような設計は、音声処理（1次元信号）や画像処理（2次元信号）の高速化には十分ですが、多次元信号の膨大なデータサンプルを扱う場合、リアルタイムで演算結果を取得するには依然として力不足です。

多次元DSP演算を高速化するために、気象予報や軍事レーダーなど、特定の状況では専用のスーパーコンピュータやクラスタコンピュータの使用が必要となる。しかしながら、DSP演算のみを行うためのスーパーコンピュータを使用するには、相当な費用とエネルギー消費がかかる。また、すべての多次元DSPアプリケーションにとって実用的かつ適切な方法とは言えない。

GPUはもともと画像処理やビデオストリームのレンダリングを高速化するために考案されました。さらに、最新のGPUは比較的低コストでエネルギー効率も高く、並列数値計算を効率的に実行できるため、多次元DSPを実行するスーパーコンピュータの代替として人気が高まっています。^{[ 1 ]}

CUDAはNvidia GPUをプログラミングするための標準インターフェースです。NvidiaはNvidia GPUデバイス上でDSPアクセラレーションをサポートするために多くのCUDAライブラリも提供しています。^{[ 4 ]}

OpenCLは、元々はApple Inc.によって提案され、現在はKhronos Groupによって維持および開発されている業界標準です。^{[ 5 ]} OpenCLは、 GPGPUを含むさまざまなデバイスを普遍的にプログラミングするためのC++ライクなAPIを提供します。

次の図は、GPUデバイス上でOpenCLプログラムを実行する際のフローを示しています。まずCPUがOpenCLデバイス（この場合はGPU）を検出し、次にジャストインタイムコンパイラを呼び出してOpenCLソースコードをターゲットバイナリに変換します。その後、CPUは計算を実行するためにデータをGPUに送信します。GPUがデータを処理している間、CPUは自身のタスクを処理できます。

C++ AMPは、 Microsoftが提案したプログラミング モデルです。C++ AMP は、SIMD プロセッサのプログラミング用に設計されたC++ベースのライブラリです^{[ 6 ]}

OpenACCは、 Cray、CAPS、NVIDIA、PGIによって開発された並列コンピューティングのプログラミング標準です。 ^{[ 7 ]} OpenAccは、 C、C++、Fortran拡張機能を備えたCPUとGPUの異種システム向けのプログラミングを対象としています。

AとBが2つのm × m行列であると仮定し、 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}}$

Cの各要素を計算するには、m 回の乗算と( m – 1 )回の加算が必要です。したがって、CPU 実装では、この計算を実行するための時間計算量は、次の C の例でΘ(n ³ )となります。しかし、 Cの要素は互いに独立していることがわかっています。したがって、GPGPU デバイスなどの SIMD プロセッサによって計算を完全に並列化できます。GPGPU 実装では、次の OpenCL の例に示すように for ループを展開することで、時間計算量はΘ(n)に大幅に削減されます。

// C言語におけるMxM行列の乗算void matrixMul (float * A 、// 入力行列 Afloat * B 、// 入力行列 Bfloat * C 、// 出力行列 Cint size ) // 行列のサイズ{// N x N x N 回の反復for ( int row = 0 ; row < size ; row ++ ) {for ( int col = 0 ; col < size ; col ++ ) {int id = row * size + col ;float sum = 0.0 ;for ( int m = 0 ; m < size ; m ++ ) {sum += ( A [ row * size + m ] * B [ m * size + col ]);}C [ id ] =合計;}}}

// OpenCLにおけるMxM行列乗算__kernel void matrixMul (__global float * A , // 入力行列 A__global float * B , // 入力行列 B__global float * C 、// 出力行列 C__global int size ) // 行列のサイズ{size_t id = get_global_id ( 0 ); // 各スレッドは要素に対して作業を行うsize_t row = id / size ;size_t col = id % size ;float sum = 0.0 ;// N回の反復for ( int m = 0 ; m < size ; m ++ ) {sum += ( A [ row * size + m ] * B [ m * size + col ]);}C [ id ] =合計;}

畳み込みはDSPで頻繁に使用される演算です。2つのm × m信号の2次元畳み込みを計算するには、出力要素ごとにm ^× ( m -1)回の乗算とm ×( m - 1 )回の加算が必要です。つまり、出力信号全体の時間計算量はΘ(n4 ⁾です。次のOpenCLの例が示すように、GPGPUアクセラレーションを使用すると、すべての出力要素がデータに依存しないため、全体の計算時間は実質的にΘ(n2 ⁾に短縮されます。

2次元畳み込み方程式：

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

// OpenCLによる2次元畳み込みの実装__kernel void convolution (__global float * x 、// 入力信号 x__global float * h , // h をフィルタリング__global float * y , // 出力信号 y__global int size ) // 入力信号とフィルタの ROS のサイズ{size_t id = get_global_id ( 0 ); // 各スレッドは要素に対して作業を行うsize_t row = size + size - 1 ; // 出力信号の行数size_t col = size + size - 1 ; // 出力信号の列数size_t n1 = id / row ;size_t n2 = id % col ;float sum = 0.0 ;// N x N 回繰り返し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 ] =合計;}

上記の例では2次元畳み込みを示しましたが、同様のアプローチはより高次元のシステムにも適用できます。全体として、sD畳み込みの場合、GPGPU実装の時間計算量はΘ(n ^s )であるのに対し、CPU実装の時間計算量はΘ(n ^2s )となります。

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

畳み込み演算に加えて、離散時間フーリエ変換（DTFT）もシステム解析でよく用いられる手法の一つである。

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

実際には、MD DTFT を実装するには、M 回 1-D DFTF を実行し、各次元に関して行列の転置を行うことができます。1-D DTFT 演算では、GPGPU は概念的に複雑さを Θ(n ² )からΘ(n )に削減できます。これは、次の OpenCL 実装例で示されています。つまり、GPGPU の複雑さの MD DTFT は、GPU 上でΘ(n ² )の複雑さで計算できます。一部の GPGPU にはハードウェア FFT アクセラレータが内部に搭載されていますが、この実装は、GPU メーカーが提供する FFT API またはライブラリを呼び出すことによって最適化することもできます。^{[ 8 ]}

// OpenCLにおけるDTFT__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 ); // 各スレッドは要素に対して作業を行う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 ));}}

多次元デジタルフィルタの設計は、特にIIRフィルタの場合、大きな課題です。通常、差分方程式を解いて近似解のセットを取得するためにコンピュータに依存します。GPGPUコンピューティングが普及するにつれて、GPGPUを使用して多次元FIRフィルタやIIRフィルタを設計するためのいくつかの適応アルゴリズムが提案されています。 ^{[ 9 ] [ 10 ] [ 11 ]}

レーダーシステムは通常、多数の 3 次元または 4 次元のデータサンプルをリアルタイムで再構築する必要があります。従来、特に軍事分野では、これにはスーパーコンピュータのサポートが必要でした。今日では、GPGPU もスーパーコンピュータの代わりにレーダー信号を処理するために使用されています。たとえば、合成開口レーダー (SAR)信号を処理するには、通常、多次元FFT計算が必要です。^{[ 12 ] [ 13 ] [ 14 ]} GPGPU は、このようなアプリケーションで FFT および/または iFFT を高速に実行するために使用できます。

多くの自動運転車は、車両の自動制御に3D画像認識技術を適用しています。急速に変化する外部環境に対応するためには、認識と意思決定のプロセスをリアルタイムで行う必要があることは明らかです。GPGPUは、この目標を達成するための優れたデバイスです。^{[ 15 ]}

正確な診断を行うためには、超音波、X線、MRI、CTなどの2次元または3次元の医療信号では、画像を再構成するために非常に高いサンプリングレートと画像解像度が必要となることがよくあります。GPGPUの優れた計算能力を適用することで、より高品質の医療画像を取得できることが示されています^{[ 16 ] [ 17 ] 。}