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 '20 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 '20 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 '20 at 10:10
  • $\begingroup$ @mikado, Thank you. Yes, xc is far from zero. But how to linearise the argument? $\endgroup$ – user14634 Jul 19 '20 at 9:59
  • 1
    $\begingroup$ in f(x) is t a static variable or is this a function of f(x,t) ? $\endgroup$ – morbo Jul 19 '20 at 10:46
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 '20 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 '20 at 13:47
  • $\begingroup$ Thanks a lot for your suggestion! $\endgroup$ – user14634 Jul 22 '20 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.