docs: 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()