docs: 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