Multidimensional discrete convolution: Difference between revisions

Revision as of 23:07, 23 February 2019

Revision as of 23:15, 14 March 2019

Line 117:

===Row-Column Decomposition===

The row-column method can be applied when one of the signals in the convolution is separable. The method exploits the properties of separability in order to achieve a method of calculating the convolution of two multidimensional signals that is more computationally efficient than direct computation of each sample (given that one of the signals are separable).<ref>{{cite journalbook|last1=Sihvo|first1=Tero|title=International Symposium on Signals, Circuits and Systems, 2005. ISSCS 2005|volume=1|pages=99–102|last2=Niittylahti|first2=Jarkko|titlechapter=Row-Column Decomposition Based 2D Transform Optimization on Subword Parallel Processors|date=5 June 2005|chapter-url=https://ieeexplore.ieee.org/abstract/document/1509860/|doi=10.1109/ISSCS.2005.1509860|isbn=978-0-7803-9029-4}}</ref> The following shows the mathematical reasoning behind the row-column decomposition approach (typically <math>h(n_1,n_2)</math> is the separable signal):

<math>y(n_1,n_2)=\sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} h(k_1,k_2)x(n_1-k_1,n_2-k_2)</math>

Line 195:

==Overlap and Add==

Another method to perform multidimensional convolution is the '''overlap and add''' approach. This method helps reduce the computational complexity often associated with multidimensional convolutions due to the vast amounts of data inherent in modern-day digital systems.<ref>{{cite journal|last1 = Fernandez|first1 = Joseph|last2 = Kumar|first2 = Vijaya|title = Multidimensional Overlap-Add and Overlap-Save for Correlation and Convolution|journal = IEEE|date = |issue = Image Processing (ICIP)|pages = 509–513|url=https://ieeexplore.ieee.org/abstract/document/6738105/|doi = 10.1109/ICIP.2013.6738105|year = 2013|isbn = 978-1-4799-2341-0}}</ref> For sake of brevity, the two-dimensional case is used as an example, but the same concepts can be extended to multiple dimensions.

Consider a two-dimensional convolution using a direct computation:

Line 257:

The overlap and save method is very similar to the overlap and add methods with a few notable exceptions. The overlap-add method involves a linear convolution of discrete-time signals, whereas the overlap-save method involves the principle of circular convolution. In addition, the overlap and save method only uses a one-time zero padding of the impulse response, while the overlap-add method involves a zero-padding for every convolution on each input component. Instead of using zero padding to prevent time-___domain aliasing like its overlap-add counterpart, overlap-save simply discards all points of aliasing, and saves the previous data in one block to be copied into the convolution for the next block.

In one dimension, the performance and storage metric differences between the two methods is minimal. However, in the multidimensional convolution case, the overlap-save method is preferred over the overlap-add method in terms of speed and storage abilities.<ref>{{Cite journal|url = |title = High-Speed Multidimensional Convolution|last = Kim|first = Chang|date = May 1980|journal = IEEE Transactions on Pattern Analysis and Machine Intelligence|doi = 10.1109/tpami.1980.4767017 |last2 = Strintzis|first2 = Michael|volume=PAMI-2|issue = 3|pages=269–273}}</ref> Just as in the overlap and add case, the procedure invokes the two-dimensional case but can easily be extended to all multidimensional procedures.

===Breakdown of Procedure===

Line 302:

<math>l_{Z''} =</math> <math>l_{Y''} +</math><math>l_{X''}</math> <math>= (M+K-1)</math><math>\times</math><math>(N+L-1)</math>

This vector length is equivalent to the dimensions of the original matrix output <math>Z</math>, making converting back to a matrix a direct transformation. Thus, the vector, <math>Z''</math>, is converted back to matrix form, which produces the output of the two-dimensional discrete convolution.<ref name=":1">{{Cite journal|url = https://www.researchgate.net/profile/Mostafa_Naghizadeh/publication/274360447_Multidimensional_convolution_via_a_1D_convolution_algorithm/links/584ed2ff08aecb6bd8d02227/Multidimensional-convolution-via-a-1D-convolution-algorithm.pdf274360447|title = Multidimensional convolution via a 1D convolution algorithm|last = Naghizadeh|first = Mostafa|date = November 2009|journal = The Leading Edge|doi = |pmid = |access-date = |last2 = Sacchi|first2 = Mauricio}}</ref>

===Filtering on a Helix===

Line 329:

<math> h(n) = -1, 0, ... , 0, -1, 4, 1, 0, ..., 0, -1, 0, ...</math>

Notice in the one-dimensional filter that there are no leading zeroes as illustrated in the one-dimensional filtering strip after being unwound. The entire one-dimensional strip could have been convolved with; however, it is less computationally expensive to simply ignore the leading zeroes. In addition, none of these backside zero values will need to be stored in memory, preserving precious memory resources.<ref name=":2">{{Cite journal|url = |title = Multidimensional recursive filters via a helix|last = Claerbout|first = Jon|date = September 1998|journal = Geophysics|volume = 63|issue = 5|doi = 10.1190/1.1444449|pmid = |access-date = |page = 9|citeseerx = 10.1.1.76.1193|bibcode = 1998Geop...63.1532C}}</ref>

===Applications===

Helix transformations to implement recursive filters via convolution are used in various areas of signal processing. Although frequency ___domain Fourier analysis is effective when systems are stationary, with constant coefficients and periodically-sampled data, it becomes more difficult in unstable systems. The helix transform enables three-dimensional post-stack migration processes that can process data for three-dimensional variations in velocity.<ref name=":2" /> In addition, it can be applied to assist with the problem of implicit three-dimensional wavefield extrapolation.<ref>{{Cite journal|title = Exploring three-dimensional implicit wavefield extrapolation with the helix transform|last = Fomel|first = Sergey|date = 1997|journal = SEP reportReport|doi = |pmid = |access-date = |last2 = Claerbout|first2 = Jon|pages = 43–60|url=https://pdfs.semanticscholar.org/5993/a76fdb37d3b1a3dfba795d6ed596c608cec9.pdf#page=179}}</ref> Other applications include helpful algorithms in seismic data regularization, prediction error filters, and noise attenuation in geophysical digital systems.<ref name=":1" />

==Gaussian Convolution==

Line 347:

===Approximation by FIR Filter===

Gaussian convolution can be effectively approximated via implementation of a [[Finite impulse response]] (FIR) filter. The filter will be designed with truncated versions of the Gaussian. For a two-dimensional filter, the transfer function of such a filter would be defined as the following:<ref name=":0">{{cite journal|last1=Getreuer|first1=Pascal|title=A Survey of Gaussian Convolution Algorithms|journal=Image Processing Onon Line|date=2013|pages=286–310|doi=10.5201/ipol.2013.87|volume=3}}</ref>

<math>H(z_1,z_2)=\frac{1}{s(r_1,r_2)} \sum_{n_1=-r_1}^{r_1}\sum_{n_2=-r_2}^{r_2}G(n_1,n_2){z_1}^{-n_1}{z_2}^{-n_2}</math>

Line 363:

<math>H(z)=\frac{1}{2r+1} \frac{z^r-z^{-r-1}}{1-z^-1}</math>

Typically, recursive passes 3, 4, or 5 times are performed in order to obtain an accurate approximation.<ref name=":0" /> A suggested method for computing ''r'' is then given as the following:<ref>{{Cite journal|title = Efficient synthesis of Gaussian filters by cascaded uniform filters|last = Wells|first = W.M.|date = 1986|journal = IEEE Transactions on Pattern Analysis and Machine Intelligence|doi = 10.1109/TPAMI.1986.4767776|pmid =|volume=PAMI-8|issue = 2|pages=234–239}}</ref>

<math>\sigma^2=\frac{1}{12}K((2r+1)^2-1)</math> where ''K'' is the number of recursive passes through the filter.