def _do_two_gaussian_fit(freqs, signal, bounds=None):
   """
   Helper function for the two gaussian fit
   """
   initial = _two_func_initializer(freqs, signal)
   # Edit out the ones we want in the order we want them: 
   initial = (initial[0], initial[1],
              initial[6], initial[7],
              initial[2], initial[3],
              initial[10], initial[11])

   # We want to preferntially weight the error on estimating the height of the
   # individual peaks, so we formulate an error-weighting function based on
   # these peaks, which is simply a two-gaussian bumpety-bump:
   w = (ut.gaussian(freqs, initial[0], 0.075, 1, 0, 0) +
        ut.gaussian(freqs, initial[1], 0.075, 1, 0, 0))

   # Further, we want to also optimize on the individual gaussians error, to
   # restrict the fit space a bit more. For this purpose, we will pass a list
   # of gaussians with indices into the parameter list, so that we can do
   # that (see mopt.err_func for the mechanics).
   func_list = [[ut.gaussian, [0,2,4,6,7],
                 ut.gaussian(freqs, initial[0], 0.075, 1, 0, 0)],
                [ut.gaussian, [1,3,5,6,7],
                 ut.gaussian(freqs, initial[1], 0.075, 1, 0, 0)]]

   params, _ = lsq.leastsqbound(mopt.err_func, initial,
                                args=(freqs, np.real(signal),
                                ut.two_gaussian, w, func_list),
                                bounds=bounds)

   return params
Helper function for the two gaussian fit