I wrote some code that allows to estimate the initial values for the fit from the data.
Also your function does not allow for certain values so these need to be constrained.
And made the function it fits both peaks. Hope it helps.
(*import data and remove string values*)
data = N@DeleteCases[
Import["C:\\Users\\mfroelin.LA018015\\Desktop\\40MHz_0001.ascii.\
csv", "Data"], "", 2];
(*Put data in list form*)
dataf = Flatten[MapIndexed[Flatten[{#2, #1}] &, data, {2}], 1];
(*plot subset of data*)
datPlot =
ListPlot3D[dataf[[;; ;; 1000]], PlotRange -> Full,
PlotStyle -> Blue];
(*define the function*)
fun = a0 +
a1 PDF[BinormalDistribution[{\[Mu]x1, \[Mu]y1}, {\[Sigma]x1, \
\[Sigma]y1}, \[Rho]1], {x, y}] +
a2 PDF[
BinormalDistribution[{\[Mu]x2, \[Mu]y2}, {\[Sigma]x2, \
\[Sigma]y2}, \[Rho]2], {x, y}];
(*find the initial condition*)
{dx, dy} = Dimensions[data];
im = Closing[
Image[Reverse@
UnitStep[GaussianFilter[data, 10] - 5 Mean[Flatten[data]]]], 3];
comps = ComponentMeasurements[
im, {"Centroid", "EquivalentDiskRadius"}];
(*peak location*)
{\[Mu]y1i, \[Mu]x1i} = First[1 /. comps];
{\[Mu]y2i, \[Mu]x2i} = First[2 /. comps];
(*peak width*)
\[Sigma]x1i = \[Sigma]y1i = Last[1 /. comps]/2;
\[Sigma]x2i = \[Sigma]y2i = Last[2 /. comps]/2;
(*baseline and amplitude*)
a0i = Mean[Flatten[data]];
a1i = a2i = 10^3 Max[data];
(*Define the function contraints*)
cons = Join[
Thread[{a0, a1, a2} > 0],
Thread[{\[Sigma]x1, \[Sigma]y1, \[Sigma]x2, \[Sigma]y2} >
0.1 Mean[{\[Sigma]x2i, \[Sigma]y2i, \[Sigma]x1i, \[Sigma]y1i}]],
Thread[1 < {\[Mu]x1, \[Mu]x2} < dx],
Thread[1 < {\[Mu]y1, \[Mu]y2} < dy],
Thread[-0.95 < {\[Rho]1, \[Rho]2} < 0.95]
]
(*fit the data on subset of points to speed up*)
model = NonlinearModelFit[dataf[[;; ;; 100]], {fun, cons},
{{a0, a0i},
{a1, a1i}, {\[Mu]x1, \[Mu]x1i}, {\[Mu]y1, \[Mu]y1i}, {\[Sigma]x1, \
\[Sigma]x2i}, {\[Sigma]y1, \[Sigma]y1i}, {\[Rho]1, 0},
{a2, a2i}, {\[Mu]x2, \[Mu]x2i}, {\[Mu]y2, \[Mu]y2i}, {\[Sigma]x2, \
\[Sigma]x2i}, {\[Sigma]y2, \[Sigma]y2i}, {\[Rho]2, 0}}, {x, y}]
(*show the fit result and model*)
sol = model["BestFitParameters"]
Simplify[fun /. sol]
(*plot the fit*)
{xmax, ymax} = Dimensions[data];
fitPlot =
Plot3D[model[x, y], {x, 1, xmax}, {y, 1, ymax}, PlotStyle -> Red,
PlotRange -> Full];
Show[datPlot, fitPlot]
