Fit a gaussian function to the difference spectra to be used for estimation
   of the GABA peak.
   
   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_gaussian(spectra, f_ppm, lb=2.6, ub=3.6):
   """
   Fit a gaussian function to the difference spectra to be used for estimation
   of the GABA peak.
   
   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 = 5
   fit_func = ut.gaussian
   # Set the bounds for the optimization
   bounds = [(lb,ub), # peak location
             (0,None), # sigma
             (0,None), # amp
             (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])
      # Use the signal for a rough estimate of the parameters for
      # initialization :
      max_idx = np.argmax(signal[ii])
      max_sig = np.max(signal[ii])
      initial_f0 = f_ppm[idx][max_idx]
      half_max_idx = np.argmin(np.abs(signal[ii] - max_sig/2))
      # We estimate sigma as the hwhm:
      initial_sigma = np.abs(initial_f0 - f_ppm[idx][half_max_idx])
      initial_off = np.min(signal[ii])
      initial_drift = 0
      initial_amp = max_sig
      
      initial = (initial_f0,
                 initial_sigma,
                 initial_amp,
                 initial_off,
                 initial_drift)
      
      params[ii], _ = lsq.leastsqbound(mopt.err_func,
                                       initial,
                                       args=(f_ppm[idx],
                                             np.real(signal[ii]),
                                             fit_func), bounds=bounds)

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