docs: def _two_func_initializer(freqs, signal):
     """
     This is a helper function for heuristic estimation of the initial parameters
     used in fitting dual peak functions
  
     _do_two_lorentzian_fit
     _do_two_gaussian_fit
     """
     # Use the signal for a rough estimate of the parameters for initialization:
     r_signal = np.real(signal)
     # The local maxima have a zero-crossing in their derivative, so we start by
     # calculating the derivative:
     diff_sig = np.diff(r_signal)
     # We look for indices that have zero-crossings (in the right direction - we
     # are looking for local maxima, not minima!)
     local_max_idx = []
     for ii in range(len(diff_sig)-1):
        if diff_sig[ii]>0 and diff_sig[ii+1]<0:
          local_max_idx.append(ii)
  
     # Array-ify it before moving on:
     local_max_idx = np.array(local_max_idx)
     # Our guesses for the location of the interesting local maxima is the two
     # with the largest signals in them: 
     max_idx = local_max_idx[np.argsort(r_signal[local_max_idx])[::-1][:2]]
     # We sort again, so that we can try to get the first one to be the left peak:
     max_idx = np.sort(max_idx)
     if len(max_idx)==1:
        max_idx = [max_idx[0], max_idx[0]]
     # And thusly: 
     max_idx_1 = max_idx[0]
     max_idx_2 = max_idx[1]
     # A few of the rest just follow:
     max_sig_1 = r_signal[max_idx_1]
     max_sig_2 = r_signal[max_idx_2]
     initial_amp_1 = max_sig_1
     initial_amp_2 = max_sig_2
     initial_f0_1 = freqs[max_idx_1]
     initial_f0_2 = freqs[max_idx_2]
     half_max_idx_1 = np.argmin(np.abs(np.real(signal) - max_sig_1/2))
     initial_hwhm_1 = np.abs(initial_f0_1 - freqs[half_max_idx_1])
     half_max_idx_2 = np.argmin(np.abs(np.real(signal) - max_sig_2/2))
     initial_hwhm_2 = np.abs(initial_f0_2 - freqs[half_max_idx_2])
     # Everything should be treated as real, except for the phase!
     initial_ph_1 = np.angle(signal[max_idx_1])
     initial_ph_2 = np.angle(signal[max_idx_2])
     # We only fit one offset and one drift, for both functions together! 
     initial_off = np.min(np.real(signal))
     initial_drift = 0
  
     initial_a_1 = (np.sum(np.real(signal)[max_idx_1:max_idx_1 +
                                           np.abs(half_max_idx_1)*2]) ) * 2
  
     initial_a_2 = (np.sum(np.real(signal)[max_idx_2:max_idx_2 +
                                           np.abs(half_max_idx_2)*2]) ) * 2
  
     return (initial_f0_1,
             initial_f0_2,
             initial_amp_1,
             initial_amp_2,
             initial_a_1,
             initial_a_2,
             initial_hwhm_1,
             initial_hwhm_2,
             initial_ph_1,
             initial_ph_2,
             initial_off,
             initial_drift)