docs: def fit_two_lorentzian(spectra, f_ppm, lb=2.6, ub=3.6):
     """
     Fit a lorentzian function to the sum spectra to be used for estimation of
     the creatine and choline peaks.
     
     Parameters
     ----------
     spectra : array of shape (n_transients, n_points)
        Typically the sum of the on/off spectra in each transient.
  
     f_ppm : array
  
     lb, ub : floats
        In ppm, the range over which optimization is bounded
     
     """
     # We are only going to look at the interval between lb and ub
     idx = ut.make_idx(f_ppm, lb, ub)
     n_points = np.abs(idx.stop - idx.start) 
     n_params = 10 # Lotsa params!
     # Set the bounds for the optimization
     bounds = [(lb,ub), #peak1 
               (lb,ub), #peak2 
               (0,None), #area1 
               (0,None), #area2 
               (0,ub-lb), #hwhm1 
               (0,ub-lb), #hwhm2
               (-np.pi/2, np.pi/2), #phase
               (-np.pi/2, np.pi/2), #phase
               (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_lorentzian_fit(f_ppm[idx], np.real(signal[ii]),
                                        bounds=bounds)
        
        model[ii] = ut.two_lorentzian(f_ppm[idx], *params[ii])
     
     return model, signal, params