Fit a gaussian function to the difference spectra

   This is useful for estimation of the Glx peak, which tends to have two
   peaks.


   Parameters
   ----------
   spectra : array of shape (n_transients, n_points)
      Typically the difference of the on/off spectra in each transient.

   f_ppm : array

   lb, ub : floats
      In ppm, the range over which optimization is bounded

def fit_two_gaussian(spectra, f_ppm, lb=3.6, ub=3.9):
   """
   Fit a gaussian function to the difference spectra

   This is useful for estimation of the Glx peak, which tends to have two
   peaks.


   Parameters
   ----------
   spectra : array of shape (n_transients, n_points)
      Typically the difference of the on/off spectra in each transient.

   f_ppm : array

   lb, ub : floats
      In ppm, the range over which optimization is bounded

   """
   idx = ut.make_idx(f_ppm, lb, ub)
   # We are only going to look at the interval between lb and ub
   n_points = idx.stop - idx.start
   n_params = 8
   fit_func = ut.two_gaussian
   # Set the bounds for the optimization
   bounds = [(lb,ub), # peak 1 location
             (lb,ub), # peak 2 location
             (0,None), # sigma 1
             (0,None), # sigma 2
             (0,None), # amp 1
             (0,None), # amp 2
             (None, None), # offset
             (None, None),  # drift
             ]

   model = np.empty((spectra.shape[0], n_points))
   signal = np.empty((spectra.shape[0], n_points))
   params = np.empty((spectra.shape[0], n_params))
   for ii, xx in enumerate(spectra):
      # We fit to the real spectrum:
      signal[ii] = np.real(xx[idx])
      params[ii] = _do_two_gaussian_fit(f_ppm[idx], np.real(signal[ii]),
                                      bounds=bounds)

      model[ii] = fit_func(f_ppm[idx], *params[ii])

   return model, signal, params