# How to perform fitting with convolution of two functions in Mathematica?

I have an experimental data as shown bellow.

data = {{1582.41939, 0}, {1582.44618, 2}, {1582.47297, 3}, {1582.49976, 2}, {1582.52655, 6}, {1582.55334, 9}, {1582.58013, 13}, {1582.60692, 21}, {1582.63371, 24}, {1582.6605, 27}, {1582.68729, 27}, {1582.71408, 29}, {1582.74087, 30}, {1582.76766, 36}, {1582.79445, 25}, {1582.82124, 54}, {1582.84803, 38}, {1582.87482, 34}, {1582.90161, 57}, {1582.9284, 74}, {1582.95519, 106}, {1582.98198, 130}, {1583.00877, 178}, {1583.03556, 213}, {1583.06235, 265}, {1583.08914, 273}, {1583.11593, 312}, {1583.14272, 287}, {1583.16951, 295}, {1583.1963, 295}, {1583.22309, 265}, {1583.24988, 255}, {1583.27667, 292}, {1583.30346, 331}, {1583.33025, 455}, {1583.35704, 566}, {1583.38383, 744}, {1583.41062, 952}, {1583.43741, 1180}, {1583.4642, 1376}, {1583.49099, 1483}, {1583.51778, 1493}, {1583.54457, 1524}, {1583.57136, 1485}, {1583.59815, 1307}, {1583.62494, 1133}, {1583.65173, 1000}, {1583.67852, 957}, {1583.70531, 921}, {1583.7321, 1063}, {1583.75889, 1303}, {1583.78568, 1626}, {1583.81247, 1999}, {1583.83926, 2459}, {1583.86605, 2675}, {1583.89284, 2909}, {1583.91963, 2917}, {1583.94642, 2965}, {1583.97321, 2498}, {1584., 2267}, {1584.02679, 1823}, {1584.05358, 1560}, {1584.08037, 1146}, {1584.10716, 911}, {1584.13395, 860}, {1584.16074, 886}, {1584.18753, 951}, {1584.21432, 1130}, {1584.24111, 1313}, {1584.2679, 1376}, {1584.29469, 1586}, {1584.32148, 1574}, {1584.34827, 1602}, {1584.37506, 1432}, {1584.40185, 1226}, {1584.42864, 1021}, {1584.45543, 826}, {1584.48222, 592}, {1584.50901, 466}, {1584.5358, 366}, {1584.56259, 303}, {1584.58938, 315}, {1584.61617, 311}, {1584.64296, 353}, {1584.66975, 425}, {1584.69654, 413}, {1584.72333, 467}, {1584.75012, 426}, {1584.77691, 411}, {1584.8037, 366}, {1584.83049, 318}, {1584.85728, 244}, {1584.88407, 203}, {1584.91086, 171}, {1584.93765, 107}, {1584.96444, 67}, {1584.99123, 67}, {1585.01802, 70}, {1585.04481, 79}, {1585.0716, 85}, {1585.09839, 92}, {1585.12518, 84}, {1585.15197, 118}, {1585.17876, 100}, {1585.20555, 84}, {1585.23234, 48}, {1585.25913, 49}, {1585.28592, 46}, {1585.31271, 29}, {1585.3395, 14}, {1585.36629, 20}, {1585.39308, 6}, {1585.41987, 8}, {1585.44666, 10}, {1585.47345, 7}, {1585.50024, 9}, {1585.52703, 8}, {1585.55382, 6}, {1585.58061, 1}};

    figexp = ListPlot[data, Joined -> False,   PlotRange -> {Automatic, {0, 3300}}, Mesh -> All,   ImageSize -> {200, 160},   PlotStyle -> {Black, PointSize[Medium]},     AspectRatio -> 0.8, FrameTicks -> {{{0, 1000, 2000, 3000}, None}, {{1583, 1584, 1585},      None}}, FrameLabel -> {"wavelength", "Intensity"},   PlotLegends -> Placed[{"Exp data"}, Above], Frame -> True]


The theoretical model for the data is the convolution of two functions : f (x) and g (x).

    f (x) = y0 + A/(w Sqrt[π/2]) Exp[-2*((x - xc)/w)^2] (Cos[(4 π c)/x t + ϕ] + 1), with  c =  2.99792458*10^5 and t = 10.0019 as two constant; other parameters  are  {y0, A, w, xc, ϕ}.

g (x) = A2/(w2 Sqrt[π/2])*Exp[-2*((x - xc2)/w2)^2], with the parameters  of  {A2, w2, xc2}.


f (x) is the spectral distribution function, which is the product of a Guassion function and a Cos function. g (x) is the filter function, which is the standard Guassion function. I check these two functions as follow.

f[x_, y0_, A_, w_, xc_, ϕ_, t_, c_] :=  y0 + A/(w Sqrt[π/2]) Exp[-2*((x - xc)/w)^2] (Cos[(4 π c)/x t + ϕ] + 1);
g[x_, A2_, xc2_, w2_] := A2/(w2 Sqrt[π/2])*Exp[-2*((x - xc2)/w2)^2];

fig1 = Plot[{f[x, 0, 0.02, 0.7, 1584, 0, 10.0019, 2.99792458*10^5], g[x, 0.0055, 1584, 0.1]},
{x, 1582.5, 1585.5}, PlotStyle -> {Red, Blue}, PlotRange -> {{1582.5, 1585.5}, All}, PlotRange -> All, ImageSize -> {200, 160}, AspectRatio -> 0.8,  FrameTicks -> {{Automatic, None}, {{1583, 1584, 1585}, None}}, Frame -> True, Axes -> False, PlotLegends -> Placed[{"f(x)", "g(x)"}, Above], FrameLabel -> {"wavelength", "Intensity"}]


To perform the fitting, I need to set a theoretical model first. However, it is difficult to obtain an analytical expression for the convolution of f (x) and g (x).

f[x_] := y0 +  A/(w Sqrt[π/2])Exp[-2*((x - xc)/w)^2] (Cos[(4 π c)/x t + ϕ] + 1);
g[x_] := A2/(w2 Sqrt[π/2])*Exp[-2*((x - xc2)/w2)^2];
Assuming[{c > 0, t > 0, y0 > 0, A > 0, xc > 0, w > 0, ϕ > 0,  A2 > 0, xc2 > 0, w2 > 0}, Convolve[f[y], g[y], y, x]]


There is no effective output for the above codes.

I also considered to perform a numerical fitting, following the method given by Origin.

dx = 0.01;
resl = 0.15;
vX = Table[i, {i, 1582, 1586, dx}];
L = Length[vX];
vX2 = PadRight[vX, 2 L - 1];
Length[vX2];
c = 2.99792458*10^5;

f[x_, y0_, A_, w_, xc_, ϕ_] := y0 + A/(w Sqrt[π/2])Exp[-2*((x - xc)/w)^2] (Cos[(4 π c)/x 10.0019 + ϕ] + 1);
f2[x_] := f[x, 0, 17, 0.67, 1583.97, 5.3];
g[x_, w2_, xc2_] := 1/(w2 Sqrt[π/2]) Exp[-2*((x - xc2)/w2)^2];
g2[x_] := g[x, resl, 1583.97];
Plot[{f2[x], g2[x]}, {x, 1582, 1586}, PlotRange -> All, ImageSize -> {200, 150}, Frame -> True,
GridLines -> {{1584, 1583.93, 1583.86}, None}, PlotStyle -> { Red, Blue}, PlotLegends -> Placed[{"f(x)", "g(x)"}, Above]]

vF = Table[f2[i], {i, 1582, 1586, dx}];
vG = Table[g2[i], {i, 1582, 1586, dx}];
ListLinePlot[{vF, vG}, PlotRange -> All, Frame -> True, GridLines -> {{1584, 1583.93, 1583.86}, None}, PlotStyle -> { Red, Blue}];

vH = ListConvolve[vG, vF, {1, 1}];
ListLinePlot[{vF, vG, vH}, PlotRange -> All, PlotStyle -> {Blue, Red, Black}];

vF2 = PadRight[vF, 2 L - 1];
vG2 = PadRight[vG, 2 L - 1];
vH2 = ListConvolve[vG2, vF2, {1, 1}];
vH3 = Take[vH2, {Floor[L/2], L + Floor[L/2] - 1}];

vH3x = Transpose[{vX, vH3}];
figH = ListLinePlot[vH3x, PlotRange -> All, PlotStyle -> Green, ImageSize -> {200, 150}, Frame -> True, GridLines -> {{1584, 1583.93, 1583.86}, None}, PlotLegends -> Placed[{"H(x)"}, Above]];
Show[{figexp, figH}]


Using the codes above, I can plot a figure similar to the experimental data. But this is not a fitting. My questions is : How to perform this fitting with convolution of two functions in Mathematica? Any help or suggestion are highly appreciated.

• You might Fourier transform your data and functions. By the convolution theorem, your convolution becomes a product. NonlinearFit might then be able to estimate the parameters. Jul 17, 2020 at 5:57
• @mikado, Thanks a lot for your good suggestion. But it is very difficult to obtain the Fourier tranfrom of the function f (x). Jul 18, 2020 at 9:41
• Perhaps if xc is not too close to zero, you could linearise the argument of the cosine? Jul 18, 2020 at 10:10
• @mikado, Thank you. Yes, xc is far from zero. But how to linearise the argument? Jul 19, 2020 at 9:59
• in f(x) is t a static variable or is this a function of f(x,t) ? Jul 19, 2020 at 10:46

This is a bit slow but you could construct the fit yourself as a minimization of square residuals.

This gives me: {A -> 18.9346, w -> 0.768869, xc -> 1583.96, ϕ -> -0.632702}

dx = 0.01;
resl = 0.15;
vX = Table[i, {i, 1582, 1586, dx}];
L = Length[vX];
vX2 = PadRight[vX, 2 L - 1];
Length[vX2];
c = 2.99792458*10^5;

f[x_, y0_, A_, w_, xc_, ϕ_] :=
y0 + A/(w Sqrt[π/2]) Exp[-2*((x - xc)/w)^2] (Cos[(4 π c)/x 10.0019 + ϕ] + 1);
g[x_, w2_, xc2_] := 1/(w2 Sqrt[π/2]) Exp[-2*((x - xc2)/w2)^2];

conv[A_?NumericQ, w_?NumericQ, xc_?NumericQ, ϕ_?NumericQ] :=
Module[{f2, g2, vF, vG},
f2 = Function[{x}, f[x, 0, A, w, xc, ϕ]];
g2 = Function[{x}, g[x, resl, xc]];
vF = PadRight[Table[f2[i], {i, 1582, 1586, dx}], 2 L - 1];
vG = PadRight[Table[g2[i], {i, 1582, 1586, dx}], 2 L - 1];
Return[
Interpolation[
Transpose[{vX,
Take[ListConvolve[vG, vF, {1, 1}], {Floor[L/2], L + Floor[L/2] - 1}]}]
]
]
]
sqresiduals[A_?NumericQ, w_?NumericQ, xc_?NumericQ, ϕ_?NumericQ] :=
With[{convintp = conv[A, w, xc, ϕ]},
Total[(#[[2]] - convintp[#[[1]]])^2 & /@ data]
]
result = NMinimize[{sqresiduals[A, w, xc, ϕ],
0 < A < 0.5, 0 < w < 2, 1580 < xc < 1586, -20 < ϕ < 20}, {A, w, xc, ϕ},
MaxIterations -> 10]
With[{convintp = conv @@ ({A, w, xc, ϕ} /. Last[result])},
ListLinePlot[{data, {#, convintp[#]} & /@ data[[All, 1]]}]
]


• Thank you so much for this good answer. Is it possible to estimate the standard error for the value of each parameter? e.g., the standard error for w = 0.768869? Jul 21, 2020 at 5:23
• The total square error is in the first result of the NMinimize, i.e First[result]`, so the mean squared error for the whole fit is just the square root of this. However, I don't know how to get the standard error per-parameter. Jul 21, 2020 at 13:47
• Thanks a lot for your suggestion! Jul 22, 2020 at 0:55