docs: def fit_gaba(self, reject_outliers=3.0, fit_lb=2.8, fit_ub=3.4,
                   phase_correct=True, fit_func=None):
          """
          Fit either a single Gaussian, or a two-Gaussian to the GABA 3 PPM
          peak.
  
          Parameters
          ----------
          reject_outliers : float
              Z-score criterion for rejection of outliers, based on their model
              parameter
  
          fit_lb, fit_ub : float
              Frequency bounds (in ppm) for the region of the spectrum to be
              fit.
  
          phase_correct : bool
              Where to perform zero-order phase correction based on the fit of
              the creatine peaks in the sum spectra
  
          fit_func : None or callable (default None).
              If this is set to `False`, an automatic selection will take place,
              choosing between a two-Gaussian and a single Gaussian, based on a
              split-half cross-validation procedure. Otherwise, the requested
              callable function will be fit. Needs to conform to the conventions
              of `fit_gaussian`/`fit_two_gaussian` and
              `ut.gaussian`/`ut.two_gaussian`.
  
          """
          # We need to fit the creatine, so that we know which transients to
          # exclude in fitting this peak:
          if not hasattr(self, 'creatine_params'):
              self.fit_creatine()
  
          fit_spectra = np.ones(self.diff_spectra.shape) * np.nan
          # Silence warnings:
          with warnings.catch_warnings():
              warnings.simplefilter("ignore")
              fit_spectra =\
                  self.diff_spectra[self._cr_transients].copy()
  
          if phase_correct:
              for ii, this_spec in enumerate(fit_spectra):
                  # Silence warnings:
                  with warnings.catch_warnings():
                      warnings.simplefilter("ignore")
                      fit_spectra[ii] = ut.phase_correct_zero(this_spec,
                              self.creatine_params[self._cr_transients][ii, 3])
  
          if fit_func is None:
              # Cross-validate!
              fitter, self.gaba_func, self.gaba_model_err, self.gaba_signal_err=\
                  self._xval_choose_funcs(fit_spectra,
                                          reject_outliers,
                                          fit_lb, fit_ub)
          # Otherwise, you had better supply a couple of callables that can be
          # used to fit these spectra!
          else:
              fitter = fit_func[0]
              self.gaba_func = fit_func[1]
              self.gaba_model_err, self.gaba_signal_err = \
                  self._xval_model_error(fit_spectra, reject_outliers,
                                         fit_lb, fit_ub, fitter, self.gaba_func)
          # Either way, we end up fitting to everything in the end: 
          choose_transients, model, signal, params, this_idx = self._fit_helper(
                                           fit_spectra, reject_outliers,
                                           fit_lb, fit_ub, fitter)
  
          self._gaba_transients = choose_transients
          self.gaba_model = model
          self.gaba_signal = signal
          self.gaba_params = params
          self.gaba_idx = this_idx
          mean_params = stats.nanmean(params, 0)
  
          self.gaba_auc =  self._calc_auc(self.gaba_func, params, self.gaba_idx)