docs: def _compute_projection_filters(G, sf, estimated_source):
      """Least-squares projection of estimated source on the subspace spanned by
      delayed versions of reference sources, with delays between 0 and
      filters_len-1
      """
      # epsilon
      eps = np.finfo(np.float).eps
  
      # shapes
      (nsampl, nchan) = estimated_source.shape
      # handles the case where we are calling this with only one source
      # G should be nsrc X nsrc X nchan X nchan X filters_len X filters_len
      # and sf should be nsrc X nchan X filters_len
      if len(G.shape) == 4:
          G = G[None, None, ...]
          sf = sf[None, ...]
      nsrc = G.shape[0]
      filters_len = G.shape[-1]
  
      # zero pad estimates and put chan in first dimension
      estimated_source = _zeropad(estimated_source.T, filters_len - 1, axis=1)
  
      # compute its FFT
      n_fft = int(2**np.ceil(np.log2(nsampl + filters_len - 1.)))
      sef = scipy.fftpack.fft(estimated_source, n=n_fft)
  
      # compute the cross-correlations between sources and estimates
      D = np.zeros((nsrc, nchan, filters_len, nchan))
      for (j, cj, c) in itertools.product(
          list(range(nsrc)), list(range(nchan)), list(range(nchan))
      ):
          ssef = sf[j, cj] * np.conj(sef[c])
          ssef = np.real(scipy.fftpack.ifft(ssef))
          D[j, cj, :, c] = np.hstack((ssef[0], ssef[-1:-filters_len:-1]))
  
      # reshape matrices to build the filters
      D = D.reshape(nsrc * nchan * filters_len, nchan)
      G = _reshape_G(G)
  
      # Distortion filters
      try:
          C = np.linalg.solve(G + eps*np.eye(G.shape[0]), D).reshape(
              nsrc, nchan, filters_len, nchan
          )
      except np.linalg.linalg.LinAlgError:
          C = np.linalg.lstsq(G, D)[0].reshape(
              nsrc, nchan, filters_len, nchan
          )
  
      # if we asked for one single reference source,
      # return just a nchan X filters_len matrix
      if nsrc == 1:
          C = C[0]
      return C