View source for Parallel multidimensional digital signal processing

{{Orphan|date=December 2016}}

'''[[Parallel algorithm|Parallel]]''' '''[[Multidimensional signal processing|multidimensional]]''' '''[[Digital electronics|digital]]''' '''[[signal processing]]''' (mD-DSP) is defined as the application of [[parallel programming]] and [[multiprocessing]] to  [[digital signal processing]] techniques to process [[Digital signal (signal processing)|digital signals]] that have more than a single dimension.  The use of mD-DSP is fundamental to many application areas such as [[Digital data|digital]] [[Digital image processing|image]] and [[Video processing|video]]<ref>Keimel, Christian, Martin Rothbucher, Hao Shen, and Klaus Diepold. "[https://mediatum.ub.tum.de/doc/1082908/file.pdf Video is a cube]." IEEE Signal Processing Magazine 28, no. 6 (2011): 41–49.</ref> [[Signal processing|processing]], [[medical imaging]], [[geophysical signal analysis]], [[sonar]], [[radar]], [[lidar]], [[array processing]], [[computer vision]], [[computational photography]], and [[augmented reality|augmented]] and [[virtual reality]]. However, as the number of [[dimension]]s of a [[Signal (electrical engineering)|signal]] increases the [[Computational complexity theory|computational complexity]] to operate on the signal increases rapidly. This relationship between the number of dimensions and the amount of [[Computational complexity of mathematical operations|complexity]],  related to both [[Time complexity|time]] and [[DSPACE|space]], as studied in the field of [[Analysis of algorithms|algorithm analysis]], is analogues to the concept of the [[curse of dimensionality]].<ref>Hanka, Rudolf, and Thomas P. Harte. "[https://www.researchgate.net/profile/Rudolf_Hanka/publication/2736962_Curse_of_Dimensionality_Classifying_Large_Multi-Dimensional_Images_with_Neural_Networks/links/54ae7d8d0cf2213c5fe578b1/Curse-of-Dimensionality-Classifying-Large-Multi-Dimensional-Images-with-Neural-Networks.pdf Curse of dimensionality: classifying large multi-dimensional images with neural networks]." In Computer Intensive Methods in Control and Signal Processing, pp. 249–260. Birkhäuser Boston, 1997.</ref> This large complexity generally results in an extremely long [[Execution (computing)|execution run-time]] of a given mD-DSP application rendering its usage to become impractical for many applications; especially for [[Real-time computing|real-time]] applications.<ref>Saponara, Sergio, Antonio Plaza, Marco Diani, Matthias F. Carlsohn, and Giovanni Corsini. "[https://link.springer.com/article/10.1007/s11554-014-0395-0 Special issue on algorithms and architectures for real-time multi-dimensional image processing]." J. Real-Time Image Processing 9, no. 3 (2014): 393–396.</ref> This long run-time is the primary motivation of applying parallel algorithmic techniques to mD-DSP problems.

== Motivation, problem statement, and basic concepts  ==
Due to the end of [[frequency scaling]]<ref>Mims, Christopher. "Why CPUs Aren't Getting Any Faster." MIT Technology Review. October 22, 2012. Accessed December 09, 2016. https://www.technologyreview.com/s/421186/why-cpus-arent-getting-any-faster/.</ref><ref>"Introduction to Parallel Programming With CUDA | Udacity." Introduction to Parallel Programming With CUDA | Udacity. Accessed December 07, 2016. https://www.youtube.com/watch?v=ET9wxEuqp_Y.</ref> of [[Microprocessor|processors]], which is largely attributed to the effect of [[Dennard scaling]]<ref>Bohr, M., 2007. [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.109.4141&rep=rep1&type=pdf#page=11 A 30 year retrospective on Dennard's MOSFET scaling paper]. IEEE Solid-State Circuits Society Newsletter, 12(1), pp. 11–13.</ref> around the year 2005, a common trend of processor manufacturers was to continue to exploit [[Moore's Law|Moore's law]]<ref>Moore, Gordon. "[https://www.supplychain247.com/images/pdfs/gordon_moore_1965_article.pdf Cramming More Components onto Integrated Circuits]," Electronics Magazine Vol. 38, No. 8 (April 19, 1965)</ref> by increasing the number of processors on a single chip, which are termed [[multi-core processor]]s as opposed to [[uniprocessor system|uniprocessors]].<ref name=H&P>Hennessy, John L., and David A. Patterson. Computer architecture: a quantitative approach. Elsevier, 2011.</ref>

mD-DSP algorithms exhibit a large amount of complexity, as described in the previous section, which makes efficient implementation difficult in regard to run-time and power consumption. This article primarily addresses basic parallel concepts used to alleviate run-time of common mD-DSP applications. The concept of [[parallel computing]] can be applied to mD-DSP applications to exploit the fact that if a problem can be expressed in a [[Parallel algorithm|parallel]] [[algorithm]]ic form, then [[parallel programming]] and [[multiprocessing]] can be used in an attempt to increase the [[High-throughput computing|computational throughput]] of the mD-DSP procedure on a given [[Computer hardware|hardware]] platform. An increase in computational throughput can result in a decreased run-time, i.e. a [[speedup]] of a specific mD-DSP algorithm.

In addition to increasing computational throughput, a generally considered equally important goal is to maximally utilize the [[memory bandwidth]] of a given computing memory architecture. The combination of the computational throughput and memory bandwidth usage can be achieved through the concept of operational intensity, which is summarized in what is referred to as the [[roofline model]].<ref>Williams, S., Waterman, A. and Patterson, D., 2009. Roofline: an insightful visual performance model for multicore architectures. Communications of the ACM, 52(4), pp. 65–76.</ref> The concepts of operational intensity and the roofline model in general have recently become popular methods of quantifying the performance of mD-DSP algorithms.<ref>Ofenbeck, Georg, et al. "[http://spiral.ece.cmu.edu:8080/pub-spiral/pubfile/ispass-2013_177.pdf Applying the roofline model]." Performance Analysis of Systems and Software (ISPASS), 2014 IEEE International Symposium on. IEEE, 2014.</ref><ref>Bernabé, Gregorio, Javier Cuenca, Luis Pedro García, and Domingo Giménez. "[https://link.springer.com/article/10.1007/s11227-014-1302-y Improving an autotuning engine for 3D Fast Wavelet Transform on manycore systems]." The Journal of Supercomputing 70, no. 2 (2014): 830–844.</ref>

Increasing throughput can be beneficial to [[Scalability#Weak versus strong scaling|strong scaling]]<ref>"Introduction to Parallel Programming With CUDA | Udacity." Introduction to Parallel Programming With CUDA | Udacity. Accessed December 07, 2016. https://www.youtube.com/watch?v=539pseJ3EUc.</ref><ref>Humphries, Benjamin, Hansen Zhang, Jiayi Sheng, Raphael Landaverde, and Martin C. Herbordt. "[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4652940/ 3D FFTs on a Single FPGA]." In Field-Programmable Custom Computing Machines (FCCM), 2014 IEEE 22nd Annual International Symposium on, pp. 68–71. IEEE, 2014.</ref> of a given mD-DSP algorithm. Another possible benefit of increasing operational intensity is to allow for an increase in weak scaling, which allows the mD-DSP procedure to operate on increased [[data]] sizes or larger [[data set]]s, which is important for application areas such as [[data mining]] and the training of [[Deep learning|deep neural networks]]<ref>Catanzaro, Bryan. "Scaling Deep Learning." 2015. Accessed December 10, 2016. http://www.highperformancegraphics.org/wp-content/uploads/2016/CatanzaroHPG2016.pdf.</ref> using [[big data]].

The goal of parallizing an algorithm is not always to decrease the traditional concept of complexity of the algorithm because the term complexity as used in this context typically refers to the [[Random-access machine|RAM]] [[Abstract machine|abstract computer model]], which by definition is serial. [[Parallel programming model|Parallel abstract computer models]] such as [[Parallel random-access machine|PRAM]] have been proposed to describe [[Analysis of parallel algorithms|complexity for parallel algorithms]] such as mD signal processing algorithms.<ref>Saybasili, A. Beliz, Alexandros Tzannes, Bernard R. Brooks, and Uzi Vishkin. "[http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.188.796&rep=rep1&type=pdf Highly Parallel Multi-Dimensional Fast Fourier Transform on Fine-and Coarse-Grained Many-Core Approaches]." In Proceedings of the 21st IASTED International Conference, vol. 668, no. 018, p. 107. 2009.</ref>

Another factor that is important to the performance of mD-DSP algorithm implementations is the resulting energy consumption and power dissipation.<ref>Ukidave, Y., Ziabari, A.K., Mistry, P., Schirner, G. and Kaeli, D., 2013, April. [http://www.northeastern.edu/esl/sites/default/files/Ukidave_ISPASS13.pdf Quantifying the energy efficiency of FFT on heterogeneous platforms]. In Performance Analysis of Systems and Software (ISPASS), 2013 IEEE International Symposium on (pp. 235–244). IEEE.</ref>

== Existing approaches ==

=== Parallel implementations of multidimensional discrete fourier transforms ===
As a simple example of an mD-DSP algorithm that is commonly decomposed into a parallel form, let’s consider the parallelization of the discrete Fourier transform, which is generally implemented using a form of the Fast Fourier Transform ([[Fast Fourier transform|FFT]]). There are hundreds of available software libraries that offer optimized FFT algorithms,<ref>"2DECOMP&FFT – Review on Fast Fourier Transform Software." 2DECOMP&FFT – Review on Fast Fourier Transform Software. Accessed December 07, 2016. http://www.2decomp.org/fft.html.</ref> and many of which offer parallelized versions of mD-FFT algorithms<ref>SIAM Journal on Scientific Computing 2012, Vol. 34, No. 4, pp. C192–C209.</ref><ref>Sandia National Laboratories: Exceptional Service in the National Interest. Accessed December 07, 2016. http://www.sandia.gov/~sjplimp/docs/fft/README.html.</ref><ref>"FFTE: A Fast Fourier Transform Package." FFTE: A Fast Fourier Transform Package. Accessed December 07, 2016. http://www.ffte.jp/.</ref> with the most popular being the parallel versions of the [[FFTW|FFTw]]<ref>Parallel FFTw. Accessed December 07, 2016. http://www.fftw.org/parallel/parallel-fftw.html.</ref> library.

The most straightforward method of paralyzing the DFT is to utilize the row-column decomposition method. The following derivation is a close paraphrasing from the classical text ''Multidimensional Digital Signal Processing''.<ref name=D&M>Dudgeon, Dan E., and Russell M. Mersereau. Multidimensional Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1984.</ref> The row-column decomposition can be applied to an arbitrary number of dimensions, but for illustrative purposes, the 2D row-column decomposition of the DFT will be described first. The 2D DFT is defined as
:<math>X \left(  k_1, k_2 \right) = \sum_{n_1=0}^{N_1-1} \sum_{n_2=0}^{N_2-1} x \left(  n_1, n_2 \right) W_{N_2}^{n_2 k_2} W_{N_1}^{n_1 k_1}</math>
where term <math>W_{N} \stackrel{\text{def}}{=}\ \exp \left( -j \frac{2\pi}{N} \right)</math> is commonly referred to as the [[twiddle factor]] of the DFT in the signal processing literature.

The DFT equation can be re-written in the following form

:<math>X(k_1, k_2) = \sum_{n_1=0}^{N_1-1} \left[ \sum_{n_2=0}^{N_2-1} \left( x(n_1,n_2) W_{N_2}^{n_2 k_2} \right) \right] W_{N_1}^{n_1 k_1}</math>

where the quantity inside the brackets is a 2D sequence which we will denote as <math>G \left( n_1, k_2 \right)</math>. We can then express the above equation as the pair of relations

:<math>G(n_1,k_2) = \sum_{n_2=0}^{N_2-1} x(n_1,n_2) W_{N_2}^{n_2 k_2} </math>

:<math>X(k_1,k_2) = \sum_{n_1=0}^{N_1-1} G(n_1,k_2) W_{N_1}^{n_1 k_1} </math>

Each column of <math> G </math> is the 1D DFT of the corresponding column of <math> x </math>. Each row of <math> X </math> is the 1D DFT of the corresponding row of <math> G </math>. Expressing the 2D-DFT in the above form allows us to see that we can compute a 2D DFT by decomposing it into row and column DFTs. The DFT of each column of <math> x </math> can first be computed where the results of which are placed into an intermediate array. Then we can compute the DFT of each row of the intermediate array.

This row-column decomposition process can easily be extended to compute an mD DFT. First, the 1D DFT is computed with respect to one of the independent variables, say <math> n_M </math>, for each value of the remaining variables. Next, 1D DFTs are computed with respect to the variable <math> n_{M-1} </math> for all values of the <math> ( M-1) </math>-tuple <math> ( n_1, n_2,\ldots, n_{M-2},k_M) </math>. We continue in this fashion until all 1D DFTs have been evaluated with respect to all the spatial variables.<ref name="D&M"/>

The row-column decomposition of the DFT is parallelized in its most simplistic manner by noting that the row and column computations are independent of each other and therefore can be performed on separate processors in parallel. The parallel 1D DFT computations on each processor can then utilize the FFT algorithm for further optimization. One large advantage of this specific method of parallelizing an mD DFT is that each of the 1D FFTs being performed in parallel on separate processors can then be performed in a [[Concurrent computing|concurrent]] fashion on [[Shared memory]] [[Multi-core processor|multithreaded SIMD processors]]<ref name=P&H>Patterson, David A., and John L. Hennessy. Computer organization and design: the hardware/software interface. Newnes, 2013.</ref>
.<ref name="H&P"/>

A specifically convenient hardware platform that has the ability to simultaneous perform both parallel and concurrent DFT implementation techniques that is highly amenable to are GPUs due to common GPUs having both a separate set of multithreaded SIMD processors (which are referred to as "streaming multiprocessors" in the CUDA programming language, and "compute units" in the OpenCL language) and individual SIMD lanes (commonly referred to loosely as a "core", or more specifically a CUDA "thread processor" or as an OpenCL "processing element") within each multithreaded SIMD processor.

A disadvantage to this technique of applying a separate FFT on each shared memory multiprocessor is the required interleaving of the data among the shared memory. One of the most popular libraries that utilizes this basic form of concurrent FFT computation is in the shared memory version of the FFTw library.<ref>"Parallel Versions of FFTW." Parallel FFTW. Accessed December 07, 2016. http://www.fftw.org/parallel/parallel-fftw.html.</ref>

=== Parallel implementations of multidimensional FIR filter structures  ===
The section will describe a method of implementing an mD digital finite impulse response (FIR) filter in a completely parallel realization. The proposed method for a completely parallel realization of a general FIR filter is achieved through the use of a combination of parallel sections consisting of cascaded 1D digital filters.<ref name=Deng>Deng, Tian-Bo. "[https://www.sciencedirect.com/science/article/pii/S105120049790292X Completely Parallel Realizations for Multidimensional Digital Filters]." Digital Signal Processing 7, no. 3 (1997): 188–198.</ref>

Consider the general desired ideal finite extent mD FIR digital filter in the complex <math>z</math>-domain, given as

:<math>H_d(z_1,z_2,\ldots,z_M) =  \sum_{n_1=0}^{N_1} \sum_{n_2=0}^{N_2} \cdots \sum_{n_M=0}^{N_M} a(n_1,n_2,\ldots,n_M) z_1^{-n_1} z_2^{-n_2} \cdots z_3^{-n_M}</math>

Placing the coefficients of <math>a(n_1,n_2,\ldots,n_M)</math> into an array <math>C=[c(l_1,l_2,\ldots,l_M)]</math> and performing some algebraic manipulation as described in,<ref name="Deng"/> we are able to arrive at an expression that allows us to decompose the filter into a parallel filterbank, given as

:<math>H_d(z_1,z_2,\ldots,z_M) \approx H(z_1,z_2,\ldots,z_M) =  \sum_{i=1}^{r} E_i(z_1) F_i(z_2) \cdots G_i(z_M)</math>

where
:<math>E_i(z_1) = \sum_{n_1=0}^{N_1}C_{i1}(n_1+1) z_1^{-n_1}</math>

:<math>F_i(z_2) = \sum_{n_2=0}^{N_3}C_{i2}(n_2+1) z_2^{-n_2}</math>

:<math>G_i(z_M) = \sum_{n_M=0}^{N_M}C_{iM}(n_M+1) z_M^{-n_M}</math>

Therefore, the original MD digital filter <math>H_d(z_1,z_2,\ldots,z_M)</math> is approximately decomposed into a parallel filterbank realization composed of a set of separable parallel filters <math>E(z_1), F(z_2), G(z_M)</math>, such that <math>H_d(z_1,z_2,\ldots,z_M) \approx H(z_1,z_2,\ldots,z_M)</math>. This proposed parallel FIR filter realization is represented by the block diagram as seen in Figure 1.
[[File:A completely parallel realization of an M-Dimensional FIR filter.png|thumb|1024px|center|Figure 1. A completely parallel realization of an M-Dimensional FIR filter.]]

The completely parallel realization as seen in figure 1 can be implemented in hardware by noting that block diagrams, and their corresponding [[Signal-flow graph]]s (SFGs) are a useful method of graphically representing any DSP algorithm that can be expressed as a linear constant coefficient difference equation. SFGs allow for easy transition from a difference equation into a hardware implementation by allowing one to visualize the difference equation in terms of digital logic components such as [[shift register]]s, and basic [[arithmetic logic unit|ALU]] digital circuit elements such as [[adder (electronics)|adder]]s and [[binary multiplier|multipliers]]. For this specific parallel realization, one could place each parallel section on a separate parallel processor to allow for each section to be implemented in a completely [[Task parallelism|task-parallel]] fashion.

Using the [[Fork–join model|fork–join]] [[parallel programming model]], a 'fork' may be applied at the first pickoff point in Figure 1, and the summing junction can be implemented during a synchronization with a 'join' operation. Implementing an mD FIR filter in this fashion lends itself well to the [[MapReduce]] [[Programming model|general programming model]]<ref>Wang, Lizhe, Dan Chen, Rajiv Ranjan, Samee U. Khan, Joanna KolOdziej, and Jun Wang. "[http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.402.1961&rep=rep1&type=pdf Parallel processing of massive EEG data with MapReduce]." In Parallel and Distributed Systems (ICPADS), 2012 IEEE 18th International Conference on, pp. 164–171. Ieee, 2012.</ref>

=== Implementations of multidimensional discrete convolution via shift registers on an FPGA ===
Convolution on mD signals lends itself well to [[Pipeline (computing)|pipelining]] due to the fact each of single output convolution operation is independent of every other one. Due to this data independence between each convolution operation between the filters impulse response and the signal a new set of data calculations may begin at the instant the first convolution operation is finished. A common method of performing mD convolution in a [[raster scan]] fashion (including dimensions greater than 2) on a traditional general purpose CPU, or even a GPU, is to cache the set of output data from each scan line of each independent dimension into the local cache. By utilizing the unique custom re-configurable architecture of a [[field-programmable gate array]] (FPGA) we can optimize this procedure dramatically by customizing the cache structure.<ref name=altera_efficient_conv>Singh, Deshanand. "Efficient Implementation of Convolutional Neural Networks Using OpenCL on FPGAs." Lecture, Embedded Vision Summit, Cali, Santa Clara, California, May 12, 2015. Accessed November 3, 2016. http://www.embedded-vision.com/platinum-members/embedded-vision-alliance/embedded-vision-training/downloads/pages/may-2015-summit-proceedings.</ref>

As in the illustrative example found in the presentation this description is derived from<ref name="altera_efficient_conv"/> we are going to restrict our discussion to two dimensional signals. In the example we perform a set of convolutional operations between a general 2D signal and a 3x3 filter kernel. As the sequence of convolution operations proceed along each raster line the filter kernel is slid across one dimension of the input signal and the data read from the memory is cached. The first pass loads three new lines of data into cache. The OpenCL code for this procedure is scene below.
<syntaxhighlight lang="c" line="1">
// Using a cache to hide poor memory access patterns on a traditional general purpose processor
for (int y = 1; y < yDim-1; ++y) {
  for (int x = 1; x < xDim-1; ++x) {
    // 3x3 Filter Kernel
    for (int y2 = -1; y2 < 1; ++y2) {
      for (int x2 = -1; x2 < 1; ++x2) {
        cache_temp[y][x] += 2D_input_signal[y+y2][x+x2]
                          * kernel[y2][x2];
        //...
      }
    }
  }
}  
</syntaxhighlight>

This caching technique is used to hide poor data to memory access pattern efficiency in terms of coalescing. However, with each successive loop only a single cache-line is updated. If we make the reasonable assumption of a 1 pixel per cycle performance point, applying this proposed caching technique to an FPGA results in a cache requirement of 9 reads and one write per cycle. Utilizing this caching technique on an FPGA results in inefficient performance in terms of both power consumption and a creation of a larger [[memory footprint]] than is required because there is a great deal of redundant reads into the cache.
With an FPGA we can customize the cache structure to give rise to a much more efficient result.

The proposed method to alleviate this poor performance with an FPGA implementation as proposed in the corresponding literature<ref name="altera_efficient_conv"/> is to customize the cache architecture through utilization of the re-configurable hardware resources of the FPGA. The important attribute to note here is that the FPGA creates the hardware based on the code that we write as opposed to writing code to run on a fixed architecture with a set of fixed [[Instruction set|instructions]].

A description of how to modify the implementation to optimize the cache architecture on an FPGA will now be discussed. Again, let's begin with an initial 2D signal and assume it is of size <math>W \times W</math>. We can then remove all of the lines that aren't in the neighborhood of the window. We next can linearize the 2D signal of this restricted segment of the 2D signal after removing lines that aren't in the neighborhood of the window. We can achieve this linearization via a simple [[Row- and column-major order|row-major]] data layout. After linearizing the 2D signal into a 1D array, under the assumption that we are not concerned with the boundary conditions of the convolution, we can discard any pixels that only contribute to the boundary computations – which is common in practice for many practical applications. Now, when we slide the window over each array element to perform the next computation, we have effectively created a shift register. The corresponding code for our shift register implementation of achieving this 2D filtering can be seen below.
<syntaxhighlight lang="c" line="1">
// Shift register in OpenCL
pixel_t sr[2*W+3];
while (keep_going) {
  // Shift data in
  #pragma unroll
  for (int i=1; i<2*W+3; ++i)
    sr[i] = sr[i-1];
  sr[0] = data_in

  // Tap output data
  data_out = {sr[  0], sr[    1], sr[    2],
              sr[  2], sr[  W+1], sr[  W+2],
              sr[2*W], sr[2*W+1], sr[2*W+2]};
  //...
} 
</syntaxhighlight>

By performing the aforementioned steps, the goal is to manage the data movement to match the FPGAs architectural strengths to achieve the highest performance. These architectural strengths allow custom implementation that can be based on the work being done as opposed to leveraging fixed operations, fixed data paths, and fixed data widths, as would be done on a general purpose processor.

==Notes==
{{reflist}}

[[Category:Digital signal processing]]