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]]]

-- demo\rosetta\Deconvolution.exw
with javascript_semantics
function m_size(sequence m)
--
-- returns the size of a matrix as a list of lengths
--
    sequence res = {}
    object me = m
    while sequence(me) do
        res &= length(me)
        me = me[1]
    end while
    return res
end function
 
function make_coordset(sequence size)
--
-- returns all points in the matrix, in zero-based indexes, 
-- eg {{0,0,0}..{3,3,5}} for a 4x4x6 matrix [96 in total]
--
    sequence res = {}
    integer count = product(size)
    for i=0 to count-1 do
        sequence coords = {}
        integer j = i
        for s=length(size) to 1 by -1 do
            integer dimension = size[s]
            coords &= mod(j,dimension)
            j = floor(j/dimension)
        end for
        coords = reverse(coords)
        res = append(res,coords)
    end for
    return res
end function
 
function row(sequence g, f, gs, gc, fs, hs)
--
--# Assembles a row, which is one of the simultaneous equations that needs
--# to be solved by reducing the whole set to reduced row echelon form. Note
--# that each row describes the equation for a single cell of the 'g' function.
--#
--# Arguments:
--# g   The "result" matrix of the convolution being undone.
--# h   The known "input" matrix of the convolution being undone.
--# gs  The size descriptor of 'g', passed as argument for efficiency.
--# gc  The coordinate in 'g' that we are generating the equation for.
--# fs  The size descriptor of 'f', passed as argument for efficiency.
--# hs  The size descriptor of 'h' (the unknown "input" matrix), passed
--#     as argument for efficiency.
--
    sequence row = {},
             coords = make_coordset(hs)
    for i=1 to length(coords) do
        sequence hc = coords[i]
        object fn = f
        for k=1 to length(gc) do
            integer d = gc[k]-hc[k]
            if d<0 or d>=fs[k] then
                fn = 0
                exit
            end if
            fn = fn[d+1]
        end for
        row = append(row,fn)
    end for
    object gn = g
    for i=1 to length(gc) do
        gn = gn[gc[i]+1]
    end for
    row = append(row,gn)
    return row
end function
 
function toRREF(sequence m)
--
-- [renamed] copy of Reduced_row_echelon_form.htm#Phix
-- plus one small tweak, as noted below, exit->return,
-- not that said seems to make any actual difference.
--
integer lead = 1,
        rows = length(m),
        cols = length(m[1])
    m = deep_copy(m)
    for r=1 to rows do
        if lead>=cols then exit end if
        integer i = r
        while m[i][lead]=0 do
            i += 1
            if i=rows then
                i = r
                lead += 1
--              if lead=cols then exit end if
                if lead=cols then return m end if
            end if
        end while
        sequence mi = deep_copy(m[r]),
                 mr = sq_div(m[i],m[i][lead])
        m[i] = mi
        m[r] = mr
        for j=1 to rows do
            if j!=r then
                m[j] = sq_sub(m[j],sq_mul(m[j][lead],m[r]))
            end if
        end for
        lead += 1
    end for
    return m
end function
 
function lset(sequence h, sequence idx, object v)
-- helper routine: store v somewhere deep inside h
    h = deep_copy(h)
    integer i1 = idx[1]+1
    if length(idx)=1 then
        h[i1] = v
    else
        h[i1] = lset(h[i1],idx[2..$],v)
    end if
    return h
end function
 
function deconvolve(sequence g, f)
--
--# Deconvolve a pair of matrixes. Solves for 'h' such that 'g = f convolve h'.
--#
--# Arguments:
--# g     The matrix of data to be deconvolved.
--# f     The matrix describing the convolution to be removed.
--  
    -- Compute the sizes of the various matrixes involved.
    sequence gsize = m_size(g),
             fsize = m_size(f),
             hsize = sq_add(sq_sub(gsize,fsize),1)
 
    -- Prepare the set of simultaneous equations to solve
    sequence toSolve = {},
             coords = make_coordset(gsize)
    for i=1 to length(coords) do
        toSolve = append(toSolve,row(g,f,gsize,coords[i],fsize,hsize))
    end for
 
    -- Solve the equations
    sequence solved = toRREF(toSolve)
 
    -- Create a result matrix of the right size
    object h = 0
    for i=length(hsize) to 1 by -1 do
        h = repeat(h,hsize[i])
    end for
 
    -- Fill the results from the equations into the result matrix
    coords = make_coordset(hsize)
    for i=1 to length(coords) do
        h = lset(h,coords[i],solved[i][$])
    end for
    return h
end function
 
constant f1 = { 6, -9, -7, -5},
         g1 = {-48, 84, -16, 95, 125, -70, 7, 29, 54, 10},
         h1 = {-8, 2, -9, -2, 9, -8, -2}
 
if deconvolve(g1, f1)!=h1 then ?9/0 end if
if deconvolve(g1, h1)!=f1 then ?9/0 end if
 
constant f2 = {{-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}},
         g2 = {{ 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}},
         h2 = {{-8, 1,-7,-2,-9, 4}, 
               { 4, 5,-5, 2, 7,-1}, 
               {-6,-3,-3,-6, 9, 5}}
 
if deconvolve(g2, f2)!=h2 then ?9/0 end if
if deconvolve(g2, h2)!=f2 then ?9/0 end if
 
constant f3 = {{{-9,  5, -8}, { 3,  5,  1}},
               {{-1, -7,  2}, {-5, -6,  6}},
               {{ 8,  5,  8}, {-2, -6, -4}}},
         g3 = {{{ 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}}},
         h3 = {{{ -6, -8, -5,  9},
                { -7,  9, -6, -8},
                {  2, -7,  9,  8}},
               {{  7,  4,  4, -6},
                {  9,  9,  4, -4},
                { -3,  7, -2, -3}}}
 
if deconvolve(g3, f3)!=h3 then ?9/0 end if
if deconvolve(g3, h3)!=f3 then ?9/0 end if
 
ppOpt({pp_Nest,2,pp_IntFmt,"%3d"})
pp(deconvolve(g3, f3))
pp(deconvolve(g3, h3))
