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

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