def _project_images(reference_sources, estimated_source, flen, G=None):
    """Least-squares projection of estimated source on the subspace spanned by
    delayed versions of reference sources, with delays between 0 and flen-1.
    Passing G as all zeros will populate the G matrix and return it so it can
    be passed into the next call to avoid recomputing G (this will only works
    if not computing permutations).
    """
    nsrc = reference_sources.shape[0]
    nsampl = reference_sources.shape[1]
    nchan = reference_sources.shape[2]
    reference_sources = np.reshape(np.transpose(reference_sources, (2, 0, 1)),
                                   (nchan*nsrc, nsampl), order='F')

    # computing coefficients of least squares problem via FFT ##
    # zero padding and FFT of input data
    reference_sources = np.hstack((reference_sources,
                                   np.zeros((nchan*nsrc, flen - 1))))
    estimated_source = \
        np.hstack((estimated_source.transpose(), np.zeros((nchan, flen - 1))))
    n_fft = int(2**np.ceil(np.log2(nsampl + flen - 1.)))
    sf = scipy.fftpack.fft(reference_sources, n=n_fft, axis=1)
    sef = scipy.fftpack.fft(estimated_source, n=n_fft)

    # inner products between delayed versions of reference_sources
    if G is None:
        saveg = False
        G = np.zeros((nchan * nsrc * flen, nchan * nsrc * flen))
        for i in range(nchan * nsrc):
            for j in range(i+1):
                ssf = sf[i] * np.conj(sf[j])
                ssf = np.real(scipy.fftpack.ifft(ssf))
                ss = toeplitz(np.hstack((ssf[0], ssf[-1:-flen:-1])),
                              r=ssf[:flen])
                G[i * flen: (i+1) * flen, j * flen: (j+1) * flen] = ss
                G[j * flen: (j+1) * flen, i * flen: (i+1) * flen] = ss.T
    else:  # avoid recomputing G (only works if no permutation is desired)
        saveg = True  # return G
        if np.all(G == 0):  # only compute G if passed as 0
            G = np.zeros((nchan * nsrc * flen, nchan * nsrc * flen))
            for i in range(nchan * nsrc):
                for j in range(i+1):
                    ssf = sf[i] * np.conj(sf[j])
                    ssf = np.real(scipy.fftpack.ifft(ssf))
                    ss = toeplitz(np.hstack((ssf[0], ssf[-1:-flen:-1])),
                                  r=ssf[:flen])
                    G[i * flen: (i+1) * flen, j * flen: (j+1) * flen] = ss
                    G[j * flen: (j+1) * flen, i * flen: (i+1) * flen] = ss.T

    # inner products between estimated_source and delayed versions of
    # reference_sources
    D = np.zeros((nchan * nsrc * flen, nchan))
    for k in range(nchan * nsrc):
        for i in range(nchan):
            ssef = sf[k] * np.conj(sef[i])
            ssef = np.real(scipy.fftpack.ifft(ssef))
            D[k * flen: (k+1) * flen, i] = \
                np.hstack((ssef[0], ssef[-1:-flen:-1])).transpose()

    # Computing projection
    # Distortion filters
    try:
        C = np.linalg.solve(G, D).reshape(flen, nchan*nsrc, nchan, order='F')
    except np.linalg.linalg.LinAlgError:
        C = np.linalg.lstsq(G, D)[0].reshape(flen, nchan*nsrc, nchan,
                                             order='F')
    # Filtering
    sproj = np.zeros((nchan, nsampl + flen - 1))
    for k in range(nchan * nsrc):
        for i in range(nchan):
            sproj[i] += fftconvolve(C[:, k, i].transpose(),
                                    reference_sources[k])[:nsampl + flen - 1]
    # return G only if it was passed in
    if saveg:
        return sproj, G
    else:
        return sproj
def _project_images ( reference_sources , estimated_source , flen , G = None ) : nsrc = reference_sources . shape [ 0 ] nsampl = reference_sources . shape [ 1 ] nchan = reference_sources . shape [ 2 ] reference_sources = np . reshape ( np . transpose ( reference_sources , ( 2 , 0 , 1 ) ) , ( nchan * nsrc , nsampl ) , order = 'F' ) reference_sources = np . hstack ( ( reference_sources , np . zeros ( ( nchan * nsrc , flen - 1 ) ) ) ) estimated_source = np . hstack ( ( estimated_source . transpose ( ) , np . zeros ( ( nchan , flen - 1 ) ) ) ) n_fft = int ( 2 ** np . ceil ( np . log2 ( nsampl + flen - 1. ) ) ) sf = scipy . fftpack . fft ( reference_sources , n = n_fft , axis = 1 ) sef = scipy . fftpack . fft ( estimated_source , n = n_fft ) if G is None : saveg = False G = np . zeros ( ( nchan * nsrc * flen , nchan * nsrc * flen ) ) for i in range ( nchan * nsrc ) : for j in range ( i + 1 ) : ssf = sf [ i ] * np . conj ( sf [ j ] ) ssf = np . real ( scipy . fftpack . ifft ( ssf ) ) ss = toeplitz ( np . hstack ( ( ssf [ 0 ] , ssf [ - 1 : - flen : - 1 ] ) ) , r = ssf [ : flen ] ) G [ i * flen : ( i + 1 ) * flen , j * flen : ( j + 1 ) * flen ] = ss G [ j * flen : ( j + 1 ) * flen , i * flen : ( i + 1 ) * flen ] = ss . T else : saveg = True if np . all ( G == 0 ) : G = np . zeros ( ( nchan * nsrc * flen , nchan * nsrc * flen ) ) for i in range ( nchan * nsrc ) : for j in range ( i + 1 ) : ssf = sf [ i ] * np . conj ( sf [ j ] ) ssf = np . real ( scipy . fftpack . ifft ( ssf ) ) ss = toeplitz ( np . hstack ( ( ssf [ 0 ] , ssf [ - 1 : - flen : - 1 ] ) ) , r = ssf [ : flen ] ) G [ i * flen : ( i + 1 ) * flen , j * flen : ( j + 1 ) * flen ] = ss G [ j * flen : ( j + 1 ) * flen , i * flen : ( i + 1 ) * flen ] = ss . T D = np . zeros ( ( nchan * nsrc * flen , nchan ) ) for k in range ( nchan * nsrc ) : for i in range ( nchan ) : ssef = sf [ k ] * np . conj ( sef [ i ] ) ssef = np . real ( scipy . fftpack . ifft ( ssef ) ) D [ k * flen : ( k + 1 ) * flen , i ] = np . hstack ( ( ssef [ 0 ] , ssef [ - 1 : - flen : - 1 ] ) ) . transpose ( ) try : C = np . linalg . solve ( G , D ) . reshape ( flen , nchan * nsrc , nchan , order = 'F' ) except np . linalg . linalg . LinAlgError : C = np . linalg . lstsq ( G , D ) [ 0 ] . reshape ( flen , nchan * nsrc , nchan , order = 'F' ) sproj = np . zeros ( ( nchan , nsampl + flen - 1 ) ) for k in range ( nchan * nsrc ) : for i in range ( nchan ) : sproj [ i ] += fftconvolve ( C [ : , k , i ] . transpose ( ) , reference_sources [ k ] ) [ : nsampl + flen - 1 ] if saveg : return sproj , G else : return sproj
Least-squares projection of estimated source on the subspace spanned by
    delayed versions of reference sources, with delays between 0 and flen-1.
    Passing G as all zeros will populate the G matrix and return it so it can
    be passed into the next call to avoid recomputing G (this will only works
    if not computing permutations).