3
$\begingroup$

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]

enter image description here

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

enter image description here

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

enter image description here

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.

$\endgroup$
13
  • 1
    $\begingroup$ 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. $\endgroup$
    – mikado
    Jul 17, 2020 at 5:57
  • $\begingroup$ @mikado, Thanks a lot for your good suggestion. But it is very difficult to obtain the Fourier tranfrom of the function f (x). $\endgroup$
    – user14634
    Jul 18, 2020 at 9:41
  • 1
    $\begingroup$ Perhaps if xc is not too close to zero, you could linearise the argument of the cosine? $\endgroup$
    – mikado
    Jul 18, 2020 at 10:10
  • $\begingroup$ @mikado, Thank you. Yes, xc is far from zero. But how to linearise the argument? $\endgroup$
    – user14634
    Jul 19, 2020 at 9:59
  • 1
    $\begingroup$ in f(x) is t a static variable or is this a function of f(x,t) ? $\endgroup$ Jul 19, 2020 at 10:46

1 Answer 1

4
+50
$\begingroup$

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

fit

$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – user14634
    Jul 21, 2020 at 5:23
  • $\begingroup$ 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. $\endgroup$
    – flinty
    Jul 21, 2020 at 13:47
  • $\begingroup$ Thanks a lot for your suggestion! $\endgroup$
    – user14634
    Jul 22, 2020 at 0:55

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.