<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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