code: def fit_two_gaussian(spectra, f_ppm, lb=3.6, ub=3.9):
     
     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
docs: 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