code: def _compute_reference_correlations(reference_sources, filters_len):
      
  
      # reshape references as nsrc X nchan X nsampl
      (nsrc, nsampl, nchan) = reference_sources.shape
      reference_sources = np.moveaxis(reference_sources, (1), (2))
  
      # zero padding and FFT of references
      reference_sources = _zeropad(reference_sources, filters_len - 1, axis=2)
      n_fft = int(2**np.ceil(np.log2(nsampl + filters_len - 1.)))
      sf = scipy.fftpack.fft(reference_sources, n=n_fft, axis=2)
  
      # compute intercorrelation between sources
      G = np.zeros((nsrc, nsrc, nchan, nchan, filters_len, filters_len))
      for ((i, c1), (j, c2)) in itertools.combinations_with_replacement(
          itertools.product(
              list(range(nsrc)), list(range(nchan))
          ),
          2
      ):
  
          ssf = sf[j, c2] * np.conj(sf[i, c1])
          ssf = np.real(scipy.fftpack.ifft(ssf))
          ss = toeplitz(
              np.hstack((ssf[0], ssf[-1:-filters_len:-1])),
              r=ssf[:filters_len]
          )
          G[j, i, c2, c1] = ss
          G[i, j, c1, c2] = ss.T
      return G, sf
docs: Compute the inner products between delayed versions of reference_sources
      reference is nsrc X nsamp X nchan.
      Returns
      * G, matrix : nsrc X nsrc X nchan X nchan X filters_len X filters_len
      * sf, reference spectra: nsrc X nchan X filters_len