This task is a straightforward generalization of Deconvolution/1D to higher dimensions. For example, the one dimensional case would be applicable to audio signals, whereas two dimensions would pertain to images. Define the discrete convolution in 





d




{\displaystyle {\mathit {d}}}

 dimensions of two functions





H
,
F
:


Z


d


→

R



{\displaystyle H,F:\mathbb {Z} ^{d}\rightarrow \mathbb {R} }

 
taking 





d




{\displaystyle {\mathit {d}}}

-tuples of integers to real numbers as the function





G
:


Z


d


→

R



{\displaystyle G:\mathbb {Z} ^{d}\rightarrow \mathbb {R} }


also taking 





d




{\displaystyle {\mathit {d}}}

-tuples of integers to reals and satisfying





G
(

n

0


,
…
,

n

d
−
1


)
=

∑


m

0


=
−
∞


∞


…

∑


m

d
−
1


=
−
∞


∞


F
(

m

0


,
…
,

m

d
−
1


)
H
(

n

0


−

m

0


,
…
,

n

d
−
1


−

m

d
−
1


)


{\displaystyle G(n_{0},\dots ,n_{d-1})=\sum _{m_{0}=-\infty }^{\infty }\dots \sum _{m_{d-1}=-\infty }^{\infty }F(m_{0},\dots ,m_{d-1})H(n_{0}-m_{0},\dots ,n_{d-1}-m_{d-1})}


for all 





d




{\displaystyle {\mathit {d}}}

-tuples of integers 



(

n

0


,
…
,

n

d
−
1


)
∈


Z


d




{\displaystyle (n_{0},\dots ,n_{d-1})\in \mathbb {Z} ^{d}}

. Assume






F




{\displaystyle {\mathit {F}}}

 and 





H




{\displaystyle {\mathit {H}}}

 (and therefore 





G




{\displaystyle {\mathit {G}}}

) are non-zero over only a finite domain bounded by the origin, hence possible to represent as finite multi-dimensional arrays or nested lists 





f




{\displaystyle {\mathit {f}}}

, 





h




{\displaystyle {\mathit {h}}}

, and 





g




{\displaystyle {\mathit {g}}}

.

For this task, implement a function (or method, procedure, subroutine, etc.) deconv to perform deconvolution (i.e., the inverse of convolution) by solving for 





h




{\displaystyle {\mathit {h}}}

 given 





f




{\displaystyle {\mathit {f}}}

 and 





g




{\displaystyle {\mathit {g}}}

. (See Deconvolution/1D for details.)

 The function should work for 





g




{\displaystyle {\mathit {g}}}

 of arbitrary length in each dimension (i.e., not hard coded or constant) and 





f




{\displaystyle {\mathit {f}}}

 of any length up to that of 





g




{\displaystyle {\mathit {g}}}

 in the corresponding dimension.
 The deconv function will need to be parameterized by the dimension 





d




{\displaystyle {\mathit {d}}}

 unless the dimension can be inferred from the data structures representing 





g




{\displaystyle {\mathit {g}}}

 and 





f




{\displaystyle {\mathit {f}}}

.
 There may be more equations than unknowns. If convenient, use a function from a library that finds the best fitting solution to an overdetermined system of linear equations (as in the Multiple regression task). Otherwise, prune the set of equations as needed and solve as in the Reduced row echelon form task.
 Debug your solution using this test data, of which a portion is shown below. Be sure to verify both that the deconvolution of 





g




{\displaystyle {\mathit {g}}}

 with 





f




{\displaystyle {\mathit {f}}}

 is 





h




{\displaystyle {\mathit {h}}}

 and that the deconvolution of 





g




{\displaystyle {\mathit {g}}}

 with 





h




{\displaystyle {\mathit {h}}}

 is 





f




{\displaystyle {\mathit {f}}}

. Display the results in a human readable form for the three dimensional case only.
dimension 1:

h: [-8, 2, -9, -2, 9, -8, -2]
f: [ 6, -9, -7, -5]
g: [-48, 84, -16, 95, 125, -70, 7, 29, 54, 10]

dimension 2:

h: [
      [-8, 1, -7, -2, -9, 4], 
      [4, 5, -5, 2, 7, -1], 
      [-6, -3, -3, -6, 9, 5]]
f: [
      [-5, 2, -2, -6, -7], 
      [9, 7, -6, 5, -7], 
      [1, -1, 9, 2, -7], 
      [5, 9, -9, 2, -5], 
      [-8, 5, -2, 8, 5]]
g: [
      [40, -21, 53, 42, 105, 1, 87, 60, 39, -28], 
      [-92, -64, 19, -167, -71, -47, 128, -109, 40, -21], 
      [58, 85, -93, 37, 101, -14, 5, 37, -76, -56], 
      [-90, -135, 60, -125, 68, 53, 223, 4, -36, -48], 
      [78, 16, 7, -199, 156, -162, 29, 28, -103, -10], 
      [-62, -89, 69, -61, 66, 193, -61, 71, -8, -30], 
      [48, -6, 21, -9, -150, -22, -56, 32, 85, 25]]

dimension 3:

h: [
      [[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], 
      [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]]
f: [
      [[-9, 5, -8], [3, 5, 1]], 
      [[-1, -7, 2], [-5, -6, 6]], 
      [[8, 5, 8],[-2, -6, -4]]]
g: [
      [
         [54, 42, 53, -42, 85, -72], 
         [45, -170, 94, -36, 48, 73], 
         [-39, 65, -112, -16, -78, -72], 
         [6, -11, -6, 62, 49, 8]], 
      [
         [-57, 49, -23, 52, -135, 66], 
         [-23, 127, -58, -5, -118, 64], 
         [87, -16, 121, 23, -41, -12], 
         [-19, 29, 35, -148, -11, 45]], 
      [
         [-55, -147, -146, -31, 55, 60], 
         [-88, -45, -28, 46, -26, -144], 
         [-12, -107, -34, 150, 249, 66], 
         [11, -15, -34, 27, -78, -50]], 
      [
         [56, 67, 108, 4, 2, -48], 
         [58, 67, 89, 32, 32, -8], 
         [-42, -31, -103, -30, -23, -8],
         [6, 4, -26, -10, 26, 12]]]

import sequtils, typetraits
 
type Size = uint64
 
type M[T: SomeNumber] = object
  dims: seq[Size]
  subsizes: seq[Size]
  data: seq[T]
 
####################################################################################################
# Miscellaneous.
 
func dotProduct[T: SomeNumber](a, b: openArray[T]): T =
  assert a.len == b.len
  for i in 0..a.high:
    result += a[i] * b[i]
 
 
####################################################################################################
# Operations on M objects.
 
func setDimensions(m: var M; dimensions: varargs[Size]) =
 
  for dim in dimensions:
    if dim == 0:
      raise newException(IndexDefect, "wrong dimension: 0")
 
  m.dims = @dimensions
  m.subsizes = m.dims
  for i in 0..dimensions.high:
    m.subsizes[i] = m.dims[(i+1)..^1].foldl(a * b, Size(1))
 
  let dlength = m.dims[0] * m.subsizes[0]
  if Size(m.data.len) != dlength:
    m.data.setLen(dlength)
 
#---------------------------------------------------------------------------------------------------
 
func initM(m: var M; dimensions: varargs[Size]) =
  m.setDimensions(dimensions)
 
#---------------------------------------------------------------------------------------------------
 
func set1DArray(m: var M; t: varargs[m.T]) =
 
  let minLen = min(m.data.len, t.len)
  m.data.setLen(minLen)
  m.data[0..<minLen] = t[0..<minLen]
 
#---------------------------------------------------------------------------------------------------
 
func seqToIdx(m: M; s: Size): seq[Size] =
 
  var acc = s
  for subsize in m.subsizes:
    result.add(acc div subsize)
    acc = acc mod subsize
 
#---------------------------------------------------------------------------------------------------
 
template size(m: M): Size = Size(m.data.len)
 
#---------------------------------------------------------------------------------------------------
 
func checkBounds(m: M; indexes: varargs[Size]): bool =
 
  if indexes.len > m.dims.len:
    return false
 
  for i, dim in indexes:
    if dim >= m.dims[i]:
      return false
 
  result = true
 
#---------------------------------------------------------------------------------------------------
 
func `[]`(m: M; indexes: varargs[Size]): m.T =
 
  if not m.checkBounds(indexes):
    raise newException(IndexDefect, "index out of range: " & $indexes)
 
  m.data[dotProduct(indexes, m.subsizes)]
 
#---------------------------------------------------------------------------------------------------
 
func `[]=`(m: M; indexes: varargs[int]; val: m.T) =
 
  if not m.checkBounds(indexes):
    raise newException(IndexDefect, "index out of range: " & $indexes)
 
  m.data[dotProduct(indexes, m.subsizes)] = val
 
#---------------------------------------------------------------------------------------------------
 
func `==`(a, b: M): bool = a.dims == b.dims and a.data == b.data
 
#---------------------------------------------------------------------------------------------------
 
func `$`(m: M): string = $m.data
 
 
####################################################################################################
# Convolution/deconvolution.
 
func convolute(h, f: M): M =
  ## Result is "g".
 
  var dims = h.dims
  for i in 0..dims.high:
    dims[i] += f.dims[i] - 1
  result.initM(dims)
 
  let bound = result.size
  for i in 0..<h.size:
    let hIndexes = h.seqToIdx(i)
 
    for j in 0..<f.size:
      let fIndexes = f.seqToIdx(j)
      for k in 0..dims.high:
        dims[k] = hIndexes[k] + fIndexes[k]
      let idx1d = dotProduct(dims, result.subsizes)
      if idx1d < bound:
        result.data[idx1d] += h.data[i] * f.data[j]
      else:
        break   # Bound reached.
 
#---------------------------------------------------------------------------------------------------
 
func deconvolute(g, f: M): M =
  ## Result is "h".
 
  var dims = g.dims
  for i, d in dims:
    if d + 1 <= f.dims[i]:
      raise newException(IndexDefect, "a dimension is zero or negative")
 
  for i in 0..dims.high:
    dims[i] -= f.dims[i] - 1
  result.initM(dims)
 
  for i in 0..<result.size:
    let iIndexes = result.seqToIdx(i)
    result.data[i] = g[iIndexes]
 
    for j in 0..<i:
      let jIndexes = result.seqToIdx(j)
      for k in 0..dims.high:
        dims[k] = iIndexes[k] - jIndexes[k]
      if f.checkBounds(dims):
        result.data[i] -= result.data[j] * f[dims]
 
    when result.T is SomeInteger:
      result.data[i] = result.data[i] div f.data[0]
    else:
      result.data[i] /= f.data[0]
 
 
####################################################################################################
# Transformation of a sequence into an M object.
 
func fold[T](a: seq[T]; d: var seq[Size]): auto =
 
  if d.len == 0:
    d.add(Size(a.len))
 
  when a.elementType is seq:
    if a.len == 0:
      raise newException(ValueError, "empty dimension")
    d.add(Size(a[0].len))
    for elem in a:
      if elem.len != a[0].len:
        raise newException(ValueError, "not rectangular")
    result = fold(a.foldl(a & b), d)
 
  else:
    if Size(a.len) != d.foldl(a * b):
      raise newException(ValueError, "not same size")
    result = a
 
#---------------------------------------------------------------------------------------------------
 
func arrtoM[T](a: T): auto =
 
  var dims: seq[Size]
  let d = fold(a, dims)
  var res: M[d.elementType]
  res.initM(dims)
  res.set1DArray(d)
  return res
 
 
#———————————————————————————————————————————————————————————————————————————————————————————————————
 
const H = @[ @[ @[-6, -8, -5,  9], @[-7,  9, -6, -8], @[ 2, -7,  9,  8] ],
             @[ @[ 7,  4,  4, -6], @[ 9,  9,  4, -4], @[-3,  7, -2, -3] ] ]
 
const F = @[ @[ @[-9,  5, -8], @[ 3,  5,  1] ],
             @[ @[-1, -7,  2], @[-5, -6,  6] ],
             @[ @[ 8,  5,  8], @[-2, -6, -4] ] ]
 
let h = arrToM(H)
let f = arrToM(F)
 
let g = h.convolute(f)
 
echo "g == f convolute h ? ", g == f.convolute(h)
echo "h == g deconv f    ? ", h == g.deconvolute(f)
echo "f == g deconv h    ? ", f == g.deconvolute(h)
echo "         f = ", f
echo "g deconv f = ", g.deconvolute(h)