1
$\begingroup$

I have some experimental data to analyse. The function $f(t)$ that describes it is however convolved with a gaussian, due to finite resolution of our experiment. $f(t)$ is too complicated to be convolved in Mathematica, and I was therefore thinking about deconvolving the data before fitting to $f(t)$, but I have been unable to do so. I have tried with ListDeconvolve and ListConvolve so far, without success. Here is the code:

step = 0.050;
d = 0;

f[t_] := (1 - E^- ((3 Abs[t])/2)) + d;
fList = Table[{step*i, f[step*i]}, {i, -50, 50}];

fConv[x_] := 
  1/Sqrt[2 \[Pi] c^2]
 Convolve[(1 - E^- ((3 Abs[t])/(2))), E^(-(t^2/(2 c^2))), 
 t, x] /. {c -> 0.5};
fConvList = Table[{step*i, fConv[step*i]}, {i, -50, 50}];

ConvList = 
  ListConvolve[GaussianMatrix[{{22}}], Table[f[step*i], {i, -50, 50}],
23];
ConvData = Table[{dd*(i - 51), %[[i]]}, {i, 1, 100}];

ListLinePlot[{fList, fConvList, ConvData}, PlotRange -> All, 
 PlotStyle -> {Normal, Dashed, Dotted}]

The problem is that if I try to get back to the original data, that is list, using ListDeconvolve, I get an absurd answer.

DeConvList = 
  ListDeconvolve[GaussianMatrix[{{22}}], ConvList, 
   Method -> "TotalVariation"];
DeConvData = Table[{step*(i - 51), DeConvList[[i]]}, {i, 1, 100}];

ListDeconvolve[GaussianMatrix[{{22}}], Table[f[step*i], {i, -25, 25}],
   Method -> "TotalVariation"];
DeConvfConv = Table[{step*(i - 51), %[[i - 51]]}, {i, 1, 100}];

ListLinePlot[{fList, fConvList, DeConvfConv, ConvData, DeConvData}, 
 PlotRange -> {0, 1.2}, 
 PlotStyle -> {Normal, Dashed, Dashed, Dotted, Dotted}]

not what I want

What is the matter?

$\endgroup$
3
$\begingroup$

Rather than deconvoluting the data before fitting, I'd suggest calculating non-convoluted target values, numerically convolving them with numerical values of your known response function, then minimizing the squared difference between these two data sets.

The process described below is a sample from code that I use to fit time-dependent fluorescence emission decays to determine the lifetime of the fluorescent species.

First of all, let's make up some data to play with:

{tmin, tmax, pitch} = {0, 25, 0.25};

(* The actual signal that we would like to fit *)
signal = With[{tau = 5}, Table[10000 Exp[-t/tau], {t, tmin, tmax, pitch}]];

(* The instrumental response function due to finite resolution, delays, etc *)
ifun = Table[10000 Exp[-(t - 2)^2/4], {t, tmin, tmax, pitch}];

(* The actual measured data: convolution between  the physical phenomenon and the *)
(* response function; add some random noise in for good measure *)
response = 10000 (Rescale@ListConvolve[ifun, signal, {-1, -1}] + 
              RandomReal[0.1, Length[signal]]);

Let's visualize the data we generated:

ListLinePlot[
 {signal, ifun, response},
 DataRange -> {tmin, tmax}, PlotRange -> All, 
 PlotStyle -> {Black, Gray, {Thick, Red}},
 PlotLegends -> {"actual signal", "instrumental response function", "measured response"}
]

data, ifun, response

Out of these three, you should know or measure your response function (and indeed you do), and you measure the response; by fitting you want to get at the parameters that best describe the actual signal.


Then I design a sum-of-squares function describing the difference between the model and the data; I will use this function as an argument to one of the Mathematica minimizers, e.g. NMinimize to obtain an estimation of the best fitting parameters.

Clear[tominimize]
tominimize[a_?NumericQ, tau_?NumericQ] := Module[
   {calculated, convoluted},
   calculated = a Exp[-#/tau] &@Range[tmin, tmax, pitch];
   convoluted = a Rescale@ListConvolve[ifun, calculated, {-1, -1}];
   Total[(convoluted - response)^2]
 ]

This function calculates a set of points using the non-convoluted model; then uses fast numerical convolution (ListConvolve) to convolve these points with the known response function, and finally compares these results to the actual experimental data, reporting the sum of the squared differences to the minimizer. Since this function will only evaluate correctly when a and tau have numerical values, NumericQ is there to prevent premature symbolic evaluation.


Using that function we can find a fit to our experimental data:

fitparams = Last@NMinimize[{tominimize[a, tau], 3 < tau < 7}, {a, tau}]

(* Out: {a -> 10548., tau -> 5.63476} *)

This is a pretty good fit to the actual values, which had been selected to be $a=10000$ and $\tau=5$.

Using these parameters we can compare the measured response to the fit results:

fittedresponse = 
  a Rescale@
     ListConvolve[ifun, a Exp[-#/tau] &@Range[tmin, tmax, pitch] /. 
       fitparams, {-1, -1}] /. fitparams;


ListPlot[
 {response, fittedresponse}, DataRange -> {tmin, tmax},
 PlotStyle -> {Directive[PointSize[0.01], Black], Directive[Thick, Darker@Green]},
 Joined -> {False, True}
]

response and fit

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.