Given my previous comments, I wonder if your data was generated by a trapezoidal wave rather than a square wave.

data={{-6.25, 0.00219649}, {-5.75, 0.0291773}, {-5.25, 0.0428364}, {-4.75, 0.0418075}, 
  {-4.25, 0.0390189}, {-3.75, 0.0142977}, {-3.25, 0.000728449}, {-2.75, 0.000577218}, 
  {-2.25, 0.00314383}, {-1.75, 0.0275447}, {-1.25, 0.0428128}, {-0.75, 0.0433703}, 
  {-0.25, 0.0433703}, {0.25, 0.0403935}, {0.75, 0.0160277}, {1.25, 0.000686316}, 
  {1.75, 0.00036222}, {2.25, 0.00288114}, {2.75, 0.0256044}, {3.25, 0.0416728}, 
  {3.75, 0.0429419}, {4.25, 0.0396677}, {4.75, 0.0162309}};

nlm = NonlinearModelFit[data,
   (a/\[Pi])*(ArcSin[Sin[(\[Pi]/m) x + l]] + ArcCos[Cos[(\[Pi]/m) x + l]]) - a/2 + c,
   {{a, 0.04}, {m, 2}, {l, 3}, {c, 0.02}}, x];
Show[ListPlot[data], Plot[nlm[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]]

The equation was lifted from Amrtanshu Raj 's answer at https://uk.mathworks.com/matlabcentral/answers/684163-generating-periodic-trapezoidal-waves-with-ramps. (Yes, a MATLAB site.)

I'm curious as to how you created the simulated data. I ask because if you sampled a square wave from $x=-6.3$ to $x=4.8$ with a constant error variance over a cycles of a square wave, it would seem near impossible to get any of those 6 points in the middle of the figure (i.e., those points between $y=0.01$ and $y=0.03$).

In other words, those 6 points could be considered outliers and maybe given very little weight.

Frankly I would find the phase separately by brute force, using reasonable candidate values (timesinbetween below). SqWave is just a pathologically nasty function to fit automatically. I would also, probably, rewrite SqWave in some other form (like SquareWave?), but here the method below doesn't really care about that:

Thanks, David. Using a Sin fit looks like a good way to define x0 and estimate starting values for a & b. There are several other techniques you included that I can make use of. I'm puzzled that the amplitude (a) turns out to be negative (not sure what that means), and I'm a bit surprised/disappointed that the fit value doesn't match the approximate max & min values of the data better than it does, because I would like to use this technique to estimate the gauge modulation ((max-min)/(max+min)).

Again, thanks much for your help - I really appreciate it.

Quite so, David, and a subjective data analysis was precisely what I was hoping to avoid by doing this fit. Thanks again. - Jim

NonlinearModelFit works better for this: