# Deconvolution using Fourier transforms

I have a 2D signal in the form of a function $$g(x_m,y_m)$$ given as

\begin{aligned} g(x_m,y_m) = \ & \int^{\infty}_{-\infty} \mathrm{d}x_o \ \mathrm{d}y_o \ \frac{1}{\varepsilon^2} P(x_o - x_m,y_o -y_m) \\ & \quad \times \int^{\infty}_{-\infty} \mathrm{d}x_i \ \mathrm{d}y_i \ R(x_o - x_i,y_o -y_i)L(x_i,y_i) \end{aligned} \tag{1}

The integrals in equation $$(1)$$ can be seen as a convolution of $$P,R$$ and $$L$$ as $$g(x_m,y_m)= ((P/\varepsilon^2)*R*L)(x_m,y_m) \tag{2}$$

I want to find $$L$$ when $$g$$,$$R$$ and $$P$$ are given.

I tried to use Fourier transforms to find $$L$$ as:

$$L = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}[g]}{\mathcal{F}[R] \cdot \mathcal{F}[P/\varepsilon^2]} \right] \tag{3}$$

R[x_, y_] := 1.04012 (0.610433 E^(-443.943 (-0.00244797 + x)^2 - 443.943 (-0.0322528 + y)^2) + 0.069016 E^(-48.5944 (-0.0105835 + x)^2 - 48.5944 (0.00545908 + y)^2) +  0.665139 E^(-426.198 (0.00315836 + x)^2 - 426.198 (0.0248491 + y)^2));

g[x_, y_] := 3.06909 E^(-18.585 x^2 + x (13.6144 - 27.7795 y) + (12.2542 - 18.1432 y) y) +  0.402245 E^(-7.37814 x^2 + x (9.61481 - 5.57202 y) + (6.46554 - 6.35048 y) y) + 120.468 E^(-0.0245919 x^2 + x (-0.00325668 + 0.00197362 y) + (-0.00103919 - 0.0281421 y) y) + 0.773818 E^(-3.79704 x^2 + (-15.0351 - 16.3606 y) y + x (1.75472 + 2.86666 y)) + 0.0833316 E^(-38.7396 x^2 + (-10.6949 - 14.1731 y) y + x (32.7513 + 18.7984 y));

epsilon = 0.048;
P[x_, y_] := E^(-(2/epsilon)^2 (x^2 + y^2));

FTR = FourierTransform[R[x, y], {x, y}, {u, v}, FourierParameters -> {0, -2 \[Pi]}];
FTg = FourierTransform[g[x, y], {x, y}, {u, v}, FourierParameters -> {0, -2 \[Pi]}];
FTP = FourierTransform[P[x, y], {x, y}, {u, v}, FourierParameters -> {0, -2 \[Pi]}];

InverseFourierTransform[FTg/((1/epsilon^2)FTP*FTR),{u, v}, {xi, yi}, FourierParameters -> {1,-2*Pi}]


However, the InverseFourierTransform takes a lot of time and does not return any result. How do I find $$L(x_i,y_i)$$?, am I doing something wrong here?.

UPDATE:

I made a naive attempt by using ListDeconvolve. Just as a test example, I neglect $$L$$ here and take the convolution equation as $$g(x_m,y_m)= ((P/\varepsilon^2)*R)(x_m,y_m)$$

To find $$R$$, for a given $$g$$ and $$P$$, I use $$R = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}[g]}{\mathcal{F}[P/\varepsilon^2]} \right] \ . \tag{4}$$

In the code below I compare the deconvolution results between the direct FourierTransform method and the ListDeconvolve method:

(* Define the equation g(xm,ym) to be deconvolved *)
g[x_, y_] := 1.04012 (0.610433 E^(-443.943 (-0.00244797 + x)^2 - 443.943 (-0.0322528 + y)^2) + 0.069016 E^(-48.5944 (-0.0105835 + x)^2 - 48.5944 (0.00545908 + y)^2) +  0.665139 E^(-426.198 (0.00315836 + x)^2 - 426.198 (0.0248491 + y)^2));

(* Define the kernel *)
epsilon = 0.048;
P[x_, y_] := E^(-(2/epsilon)^2 (x^2 + y^2));

(* Take the Fourier Transforms and carryout the operation given by equation 4 *)
FTg = FourierTransform[g[x, y], {x, y}, {u, v}, FourierParameters -> {0, -2*Pi}];
FTP = FourierTransform[P[x, y], {x, y}, {u, v}, FourierParameters -> {0, -2*Pi}];
R = InverseFourierTransform[FTg/FTP, {u, v}, {x, y}, FourierParameters -> {0, -2*Pi}];

(* Normalize the deconvolved function to its maximum value *)
NomalizedR = (First[NMaximize[R, {x, y}]])^-1 R;

(* Sample the function g(xm,ym) and deconvolve using ListDeconvolve *)
Sampledg = Table[g[x, y], {x, -0.2, 0.2, 0.005}, {y, -0.2, 0.2, 0.005}];
gData = Table[{x, y, g[x, y]}, {x, -0.2, 0.2, 0.005}, {y, -0.2, 0.2, 0.005}];
SampledP = Table[P[x, y], {x, -0.048046875/2, 0.048046875/2, 0.048046875/8}, {y, -0.048046875/2, 0.048046875/2, 0.048046875/8}];
DeconvolvedR = ListDeconvolve[SampledP, Sampledg];
DataDeconR = Transpose[{Flatten[gData, 1][[All, 1]], Flatten[gData, 1][[All, 2]], Flatten[DeconvolvedR]}];

(* Fit a nonlinear model to the ListDeconvolve result in order to compare the graphs *)
modelDeconR =  E^(-((x - c1[1])^2/b1[1]^2) - (y - c2[1])^2/b2[1]^2) a[1] +
E^(-((x - c1[2])^2/b1[2]^2) - (y - c2[2])^2/b2[2]^2) a[2] +
E^(-((x - c1[3])^2/b1[3]^2) - (y - c2[3])^2/b2[3]^2) a[3];

Param = {a[1], a[2], a[3], b1[1], b1[2], b1[3], b2[1], b2[2], b2[3], {c1[1], 0.01}, {c1[2], 0.012}, {c1[3], 0.01}, {c2[1], 0.00}, {c2[2], 0.08}, {c2[3], 0.007}};

nlm = NonlinearModelFit[DataDeconR, modelDeconR, Param, {x, y},  MaxIterations -> 5000];
FitParam = nlm["BestFitParameters"]

(* Normalize the fitted model to its maximum value *)
NomalizedDeconR = (First[NMaximize[modelDeconR /. FitParam, {x, y}]])^-1 modelDeconR /.FitParam;

(* Compare the results in a 1D plot *)
LogPlot[{g[x, 0], NomalizedR /. {y -> 0}, NomalizedDeconR /. {y -> 0}}, {x, -0.25, 0.25}, PlotRange -> All]


The graph shows a little difference between the FourierTransform method and ListDeconvolve method. How far correct is this?. Can I extend the same process to equation $$3$$ as well?.