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.

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