def coil_combine(data, w_idx=[1,2,3], coil_dim=2, sampling_rate=5000.):
    """
    Combine data across coils based on the amplitude of the water peak,
    according to:

    .. math::
        
        X = \sum_{i}{w_i S_i}

   Where X is the resulting combined signal, $S_i$ are the individual coil
   signals and $w_i$ are calculated as:

   .. math::
   
        w_i = mean(S_i) / var (S_i)
        
    following [Hall2013]_. In addition, we apply a phase-correction, so that
    all the phases of the signals from each coil are 0

    Parameters
    ----------
    data : float array
       The data as it comes from the scanner, with shape (transients, echos,
    coils, time points) 
    
    w_idx : list
       The indices to the non-water-suppressed transients. Per default we take
        the 2nd-4th transients. We dump the first one, because it seems to be
        quite different than the rest of them...

    coil_dim : int
        The dimension on which the coils are represented. Default: 2

    sampling rate : float
        The sampling rate in Hz. Default : 5000.
  
    References
    ----------

    .. [Hall2013] Emma L. Hall, Mary C. Stephenson, Darren Price, Peter
       G. Morris (2013). Methodology for improved detection of low
       concentration metabolites in MRS: Optimised combination of signals from 
       multi-element coil arrays. Neuroimage 86: 35-42.  

    .. [Wald1997] Wald, L. and Wright, S. (1997). Theory and application of
       array coils in MR spectroscopy. NMR in Biomedicine, 10: 394-410.

    .. [Keeler2005] Keeler, J (2005). Understanding NMR spectroscopy, second
       edition. Wiley (West Sussex, UK).
    
    """
    w_data, w_supp_data = separate_signals(data, w_idx)

    fft_w = np.fft.fftshift(fft.fft(w_data))
    fft_w_supp = np.fft.fftshift(fft.fft(w_supp_data))
    freqs_w = np.linspace(-sampling_rate/2.0,
                          sampling_rate/2.0,
                          w_data.shape[-1])
    
    # To determine phase and amplitude, fit a Lorentzian line-shape to each
    # coils data in each trial: 
    # No bounds except for on the phase:
    bounds = [(None,None),
              (0,None),
              (0,None),
              (-np.pi, np.pi),
              (None,None),
              (None, None)]

    n_params = len(bounds)
    params = np.zeros(fft_w.shape[:-1] + (n_params,))

    # Let's fit a Lorentzian line-shape to each one of these:
    for repeat in range(w_data.shape[0]):
       for echo in range(w_data.shape[1]):
          for coil in range(w_data.shape[2]):
             sig = fft_w[repeat, echo, coil]
             # Use the private function to do this:
             params[repeat, echo, coil] = _do_lorentzian_fit(freqs_w,
                                                             sig, bounds)


    # The area parameter stands for the magnitude:
    area_w = params[..., 1]

    # In each coil, we derive S/(N^2):
    s = np.mean(area_w.reshape(-1, area_w.shape[-1]), 0)
    n = np.var(area_w.reshape(-1, area_w.shape[-1]), 0)
    amp_weight = s/n 
    # Normalize to sum to 1: 
    amp_weight = amp_weight / np.sum(amp_weight)    
    
    # Next, we make sure that all the coils have the same phase. We will use
    # the phase of the Lorentzian to align the phases:
    phase_param = params[..., 3]
    zero_phi_w = np.mean(phase_param.reshape(-1, phase_param.shape[-1]),0)

    # This recalculates the weight with the phase alignment (see page 397 in
    # Wald paper):
    weight = amp_weight * np.exp(-1j * zero_phi_w) 

    # Multiply each one of the signals by its coil-weights and average across
    # coils:
    na = np.newaxis  # Short-hand

    # Collapse across coils for the combination in both the water 
    weighted_w_data = np.mean(np.fft.ifft(np.fft.fftshift(
       weight[na, na, :, na] * fft_w)), coil_dim)
    weighted_w_supp_data = np.mean(np.fft.ifft(np.fft.fftshift(
       weight[na, na, : ,na] * fft_w_supp)) , coil_dim)

    # Normalize each series by the sqrt(rms):
    def normalize_this(x):
       return  x * (x.shape[-1] / (np.sum(np.abs(x))))

    weighted_w_data = normalize_this(weighted_w_data)
    weighted_w_supp_data = normalize_this(weighted_w_supp_data)
    # Squeeze in case that some extraneous dimensions were introduced (can
    # happen for SV data, for example)
    return weighted_w_data.squeeze(), weighted_w_supp_data.squeeze()