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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
 
double PI;
typedef double complex cplx;
 
void _fft(cplx buf[], cplx out[], int n, int step)
{
	if (step < n) {
		_fft(out, buf, n, step * 2);
		_fft(out + step, buf + step, n, step * 2);
 
		for (int i = 0; i < n; i += 2 * step) {
			cplx t = cexp(-I * PI * i / n) * out[i + step];
			buf[i / 2]     = out[i] + t;
			buf[(i + n)/2] = out[i] - t;
		}
	}
}
 
void fft(cplx buf[], int n)
{
	cplx out[n];
	for (int i = 0; i < n; i++) out[i] = buf[i];
	_fft(buf, out, n, 1);
}
 
/* pad array length to power of two */
cplx *pad_two(double g[], int len, int *ns)
{
	int n = 1;
	if (*ns) n = *ns;
	else while (n < len) n *= 2;
 
	cplx *buf = calloc(sizeof(cplx), n);
	for (int i = 0; i < len; i++) buf[i] = g[i];
	*ns = n;
	return buf;
}
 
void deconv(double g[], int lg, double f[], int lf, double out[], int row_len) {
	int ns = 0;
	cplx *g2 = pad_two(g, lg, &ns);
	cplx *f2 = pad_two(f, lf, &ns);
 
	fft(g2, ns);
	fft(f2, ns);
 
	cplx h[ns];
	for (int i = 0; i < ns; i++) h[i] = g2[i] / f2[i];
	fft(h, ns);
 
	for (int i = 0; i < ns; i++) {
		if (cabs(creal(h[i])) < 1e-10) 
			h[i] = 0;
	}
 
	for (int i = 0; i > lf - lg - row_len; i--)
		out[-i] = h[(i + ns) % ns]/32;
	free(g2);
	free(f2);
}
 
double* unpack2(void *m, int rows, int len, int to_len)
{
	double *buf = calloc(sizeof(double), rows * to_len);
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < len; j++)
			buf[i * to_len + j] = ((double(*)[len])m)[i][j];
	return buf;
}
 
void pack2(double * buf, int rows, int from_len, int to_len, void *out)
{
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < to_len; j++)
			((double(*)[to_len])out)[i][j] = buf[i * from_len + j] / 4;
}
 
void deconv2(void *g, int row_g, int col_g, void *f, int row_f, int col_f, void *out) {
	double *g2 = unpack2(g, row_g, col_g, col_g);
	double *f2 = unpack2(f, row_f, col_f, col_g);
 
	double ff[(row_g - row_f + 1) * col_g];
	deconv(g2, row_g * col_g, f2, row_f * col_g, ff, col_g);
	pack2(ff, row_g - row_f + 1, col_g, col_g - col_f + 1, out);
 
	free(g2);
	free(f2);
}
 
double* unpack3(void *m, int x, int y, int z, int to_y, int to_z)
{
	double *buf = calloc(sizeof(double), x * to_y * to_z);
	for (int i = 0; i < x; i++)
		for (int j = 0; j < y; j++) {
			for (int k = 0; k < z; k++)
				buf[(i * to_y + j) * to_z + k] =
					((double(*)[y][z])m)[i][j][k];
		}
	return buf;
}
 
void pack3(double * buf, int x, int y, int z, int to_y, int to_z, void *out)
{
	for (int i = 0; i < x; i++)
		for (int j = 0; j < to_y; j++)
			for (int k = 0; k < to_z; k++)
				((double(*)[to_y][to_z])out)[i][j][k] =
					buf[(i * y + j) * z + k] / 4;
}
 
void deconv3(void *g, int gx, int gy, int gz, void *f, int fx, int fy, int fz, void *out) {
	double *g2 = unpack3(g, gx, gy, gz, gy, gz);
	double *f2 = unpack3(f, fx, fy, fz, gy, gz);
 
	double ff[(gx - fx + 1) * gy * gz];
	deconv(g2, gx * gy * gz, f2, fx * gy * gz, ff, gy * gz);
	pack3(ff, gx - fx + 1, gy, gz, gy - fy + 1, gz - fz + 1, out);
 
	free(g2);
	free(f2);
}
 
int main()
{
	PI = atan2(1,1) * 4;
	double h[2][3][4] = {
		{{-6, -8, -5,  9}, {-7, 9, -6, -8}, { 2, -7,  9,  8}},
		{{ 7,  4,  4, -6}, { 9, 9,  4, -4}, {-3,  7, -2, -3}}
	};
	int hx = 2, hy = 3, hz = 4;
	double f[3][2][3] = {	{{-9,  5, -8}, { 3,  5,  1}},
				{{-1, -7,  2}, {-5, -6,  6}},
				{{ 8,  5,  8}, {-2, -6, -4}} };
	int fx = 3, fy = 2, fz = 3;
	double g[4][4][6] = {
		{	{ 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}
		}
	};
	int gx = 4, gy = 4, gz = 6;
 
	double h2[gx - fx + 1][gy - fy + 1][gz - fz + 1];
	deconv3(g, gx, gy, gz, f, fx, fy, fz, h2);
	printf("deconv3(g, f):\n");
	for (int i = 0; i < gx - fx + 1; i++) {
		for (int j = 0; j < gy - fy + 1; j++) {
			for (int k = 0; k < gz - fz + 1; k++)
				printf("%g ", h2[i][j][k]);
			printf("\n");
		}
		if (i < gx - fx) printf("\n");
	}
 
	double f2[gx - hx + 1][gy - hy + 1][gz - hz + 1];
	deconv3(g, gx, gy, gz, h, hx, hy, hz, f2);
	printf("\ndeconv3(g, h):\n");
	for (int i = 0; i < gx - hx + 1; i++) {
		for (int j = 0; j < gy - hy + 1; j++) {
			for (int k = 0; k < gz - hz + 1; k++)
				printf("%g ", f2[i][j][k]);
			printf("\n");
		}
		if (i < gx - hx) printf("\n");
	}
}
 
/* two-D case; since task doesn't require showing it, it's commented out */
/*
int main()
{
	PI = atan2(1,1) * 4;
	double h[][6] = { 	{-8, 1, -7, -2, -9, 4},
				{4, 5, -5, 2, 7, -1},
				{-6, -3, -3, -6, 9, 5} };
	int hr = 3, hc = 6;
 
	double f[][5] = {	{-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} };
	int fr = 5, fc = 5;
	double g[][10] = {
			{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}	};
	int gr = 7, gc = 10;
 
	double h2[gr - fr + 1][gc - fc + 1];
	deconv2(g, gr, gc, f, fr, fc, h2);
	for (int i = 0; i < gr - fr + 1; i++) {
		for (int j = 0; j < gc - fc + 1; j++)
			printf(" %g", h2[i][j]);
		printf("\n");
	}
 
	double f2[gr - hr + 1][gc - hc + 1];
	deconv2(g, gr, gc, h, hr, hc, f2);
	for (int i = 0; i < gr - hr + 1; i++) {
		for (int j = 0; j < gc - hc + 1; j++)
			printf(" %g", f2[i][j]);
		printf("\n");
	}
}*/