Python signal.convolve函数代码示例

def degrades(V, W, H, rate, s, niter=20):
    """
    Deconvolution by Gradient Descent with Sparsification.
    """
    V, W, H = normalize(V), normalize(W), normalize(H)
    for iter in xrange(niter):
        convolved_pieces = np.vstack([signal.convolve(W[r], H[r], mode='full') for r in xrange(W.shape[0])])
        convolved = np.sum(convolved_pieces, axis=0)
        delta = V - convolved
        projected_H = fft_deconvolve(delta, W)
        H = magical_entropy_hammer(H, H + projected_H * rate, 0.01, axis=1)
        
        
        convolved_pieces = np.vstack([signal.convolve(W[r], H[r], mode='full') for r in xrange(W.shape[0])])
        convolved = np.sum(convolved_pieces, axis=0)
        delta = V - convolved
        projected_W = fft_deconvolve(delta, H)
        W = magical_entropy_hammer(W, W + projected_W * rate, 0.001, axis=1)
        
        print np.sum(np.abs(delta))
        print H
        
    pylab.clf()
    pylab.subplot(311)
    pylab.plot(W.T)
    pylab.subplot(312)
    pylab.plot(H.T)
    pylab.subplot(313)
    pylab.plot(V)
    pylab.plot(convolved)
    return W, H

def test3DConvolution():
    kernel = array([[[-1, 1], [1, 0]], [[1, 0], [0, 0]]])
    kernel_cube = array([
        zeros((3,3)),
        [[0, 0, 0], [0, -1, 1], [0, 1, 0]],
        [[0, 0, 0], [0,  1, 0], [0, 0, 0]]])

    A = reshape(linspace(1,9,9), (3,3), order = 'C')
    B = array((A, A, A))

    C = operators.convolve(kernel, (3, 3, 3), order = 'F')
    adjointTest(C)

    np.testing.assert_allclose(
        reshape(C._forward(ravel(B, order = 'F')), (3, 3, 3), order = 'F'),
        signal.convolve(B, kernel, 'same'))

    np.testing.assert_allclose(
        reshape(C._forward(ravel(B, order = 'F')), (3, 3, 3), order = 'F'),
        signal.convolve(B, kernel_cube, 'same'))

    np.testing.assert_allclose(
        reshape(C._adjoint(ravel(B, order = 'F')), (3, 3, 3), order = 'F'),
        signal.convolve(B,
            np.flipud(np.fliplr(kernel_cube[:,:,::-1])), 'same'))

def get_corr (im, pars):

    # transform 2-d image to 1-d total signal vectors
    totsig, totsigla = get_totsig(im,pars.la)

    ln1 = im.npix / 2
    ln2 = im.npix
    quad_len = ln1*ln1

    # compute staypuft signal for each pixel
    mask = totsig > 40.0
    p0 = np.fabs(pars.ampscale * mask * totsig)
    p1 = np.fabs(pars.ampscale * mask * totsigla)

    ekern = np.exp(-np.arange(ln1*pars.tx)/pars.hh)
    qkern = pars.a1*np.arange(ln1*pars.tx) + pars.a2

    e = convolve (p0, ekern, mode='full')
    q = convolve (p1, qkern, mode='full')
    b = e[0:quad_len] + q[0:quad_len]

    # transform the correction vector back into a 2-d image quad
    b = b[::-1]
    b = np.reshape (b, (ln1,ln1))
    b = np.transpose(b)

    # replicate the correction into all 4 full image quads
    im.data[0:ln1,0:ln1] = b
    im.data[0:ln1,ln1:ln2] = b
    im.data[ln1:ln2,0:ln1] = b
    im.data[ln1:ln2,ln1:ln2] = b

    return im

 def test_consistency_convolve_funcs(self):
     # Compare np.convolve, signal.convolve, signal.convolve2d
     a = np.arange(5)
     b = np.array([3.2, 1.4, 3])
     for mode in ["full", "valid", "same"]:
         assert_almost_equal(np.convolve(a, b, mode=mode), signal.convolve(a, b, mode=mode))
         assert_almost_equal(np.squeeze(signal.convolve2d([a], [b], mode=mode)), signal.convolve(a, b, mode=mode))

def basefreq(audiofile):
    """
    This function reads in the audio file and does the hann windowed fft of 
    the right input. It then smooths the output using a gaussian filter and
    then finds the peaks. It returns the peaks in the right audio channel since
    testing showed there was no significant difference in the two.
    """
    #read the data into an ndarray using scikits-audiolab        
    data, rate, enc = al.aiffread(audiofile)
    #split the left and right channel
    datar = data[:,1]
    datal = data[:,0]
    #take the fft of both of the channels with the hann window applied
    #the hann window reduces spectral leakage in the FFT     
    dftr = abs(fft.fft(datar*signal.hann(len(datar))))
    dftl = abs(fft.fft(datal*signal.hann(len(datal))))
    #compute the frequencies in the FFT
    freq = float(rate)/float(len(datar))
    freqs = np.arange(len(dftr)/2+99)*freq
    dftr = dftr[0:np.size(dftr)/2]
    dftl = dftl[0:np.size(dftr)/2]
    #smooth the fft with a gaussian
    c = signal.gaussian(100,20)
    dftr = signal.convolve(dftr,c)
    dftl = signal.convolve(dftl,c)
    #find the significant peaks in each channel
    peaksr = findpeaks(dftr,freqs)
    peaksl = findpeaks(dftl,freqs)
    #plot the output fft for testing
    #plt.plot(freqs,dftr)
    #plt.show()
    #print peaksr
    return peaksr

def compute_harris_response(image):
    """ compute the Harris corner detector response function 
        for each pixel in the image"""

    #derivatives
    imx, imy = filtertools.gauss_derivatives(image, 3)

    #kernel for blurring
    gauss = filtertools.gauss_kernel(3)

    #compute components of the structure tensor
    Wxx = signal.convolve(imx*imx, gauss, mode='same')
    Wxy = signal.convolve(imx*imy, gauss, mode='same')
    Wyy = signal.convolve(imy*imy, gauss, mode='same')

    #determinant and trace
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    if numpy.count_nonzero(Wtr) == 0:
        return
    print Wtr.shape
    print numpy.count_nonzero(Wtr)

    return Wdet / Wtr

def ddwt(x, num_scales):
    """
    Дискретное вейвлет-преобразование без прореживания
    :param x: входной сигнал
    :param num_scales: число уровней разложения
    :return:
    """
    h = np.array([1, 3, 3, 1], float) / 8
    g = np.array([2, -2], float)
    signal_len = len(x)

    detail = []
    approx = []
    ap = x.copy()
    detail.append(ap.copy())  # на нулевом уровне храним исходный сигнал
    approx.append([])

    for s in range(num_scales):
        dly = 2**s
        hif = convolve(ap, g, mode="full")[dly:dly+signal_len]
        detail.append(hif)
        approx.append(ap)
        if s < num_scales-1:
            ap = convolve(ap, h, mode="full")[dly:dly+signal_len]
            dly_lo = len(h)-1
            ap[:dly_lo] = ap[dly_lo] # хак
            # вместо прореживания сигналов (Маллат) расширяем характеристики
            # фильтров
            h = fexpand(h)
            g = fexpand(g)

    return approx, detail

def lsf_to_lpc(all_lsf):
    if len(all_lsf.shape) < 2:
        all_lsf = all_lsf[None]
    order = all_lsf.shape[1]
    all_lpc = np.zeros((len(all_lsf), order + 1))
    for i in range(len(all_lsf)):
        lsf = all_lsf[i]
        zeros = np.exp(1j * lsf)
        sum_zeros = zeros[::2]
        diff_zeros = zeros[1::2]
        sum_zeros = np.hstack((sum_zeros, np.conj(sum_zeros)))
        diff_zeros = np.hstack((diff_zeros, np.conj(diff_zeros)))
        sum_filt = np.poly(sum_zeros)
        diff_filt = np.poly(diff_zeros)

        if order % 2 != 0:
            deconv_diff = sg.convolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff = sg.convolve(diff_filt, [1, -1])
            deconv_sum = sg.convolve(sum_filt, [1, 1])

        lpc = .5 * (deconv_sum + deconv_diff)
        # Last coefficient is 0 and not returned
        all_lpc[i] = lpc[:-1]
    return np.squeeze(all_lpc)

    def convolve(self, f):
        from scipy.signal import convolve

        e_field = self.E_field_as_numpy()
        f_data = np.zeros((self.dim_x(), self.dim_y()))
        result = np.zeros((self.numberEnergies(), self.dim_x(), self.dim_y(), 2), dtype=np.complex128)

        print("Convolve I_X ", end="")
        for i_x, x_cooridinate in enumerate(self.absolute_x_coordinates()):
            if i_x%100 ==0:
                print(" ",i_x , end="")
            for i_y, y_cooridinate in enumerate(self.absolute_y_coordinates()):
                f_data[i_x, i_y] = f(x_cooridinate,y_cooridinate)

        for index_energy in range(self.numberEnergies()):
            for pol in (0,1):
                print("Convolving pol", pol)
                #r = convolve(f_data, f_data,mode='same')
                r = convolve(e_field[index_energy,:,:,pol].real, f_data,mode='same')
                r = r + 1j *convolve(e_field[index_energy,:,:,pol].imag, f_data,mode='same')

                for i_x, x_cooridinate in enumerate(self.absolute_x_coordinates()):
                    for i_y, y_cooridinate in enumerate(self.absolute_y_coordinates()):
                        result[index_energy, i_x , i_y , pol] = r[i_x,i_y]

        convolved_wavefront = NumpyWavefront(result, self.x_start(), self.x_end(), self.y_start(), self.y_end(), self.z())

        return convolved_wavefront

    def __init__(self, shapein, kernel, mode="full", fft=False, **kwargs):
        if fft:
            from scipy.signal import fftconvolve as convolve
            # check kernel shape parity
            if np.any(np.asarray(kernel.shape) % 2 != 1):
                raise ValueError("Kernels with non-even shapes are not handled for now.")
        else:
            from scipy.signal import convolve

        self.kernel = kernel
        self.mode = mode

        # reverse kernel
        s = (slice(None, None, -1), ) * kernel.ndim
        self.rkernel = kernel[s]
        # reverse mode
        if mode == 'full':
            self.rmode = 'valid'
        elif mode == 'valid':
            self.rmode = 'full'
        elif mode == 'same':
            self.rmode = 'same'
        # shapeout
        if mode == 'full':
            shapeout = [s + ks - 1 for s, ks in zip(shapein, kernel.shape)]
        if mode == 'valid':
            shapeout = [s - ks + 1 for s, ks in zip(shapein, kernel.shape)]
        if mode == 'same':
            shapeout = shapein

        matvec = lambda x: convolve(x, self.kernel, mode=self.mode)
        rmatvec = lambda x: convolve(x, self.rkernel, mode=self.rmode)
        NDOperator.__init__(self, shapein, shapeout, matvec, rmatvec, **kwargs)

    def test_input_swapping(self):
        small = arange(8).reshape(2, 2, 2)
        big = 1j * arange(27).reshape(3, 3, 3)
        big += arange(27)[::-1].reshape(3, 3, 3)

        out_array = array(
            [[[0 + 0j, 26 + 0j, 25 + 1j, 24 + 2j],
              [52 + 0j, 151 + 5j, 145 + 11j, 93 + 11j],
              [46 + 6j, 133 + 23j, 127 + 29j, 81 + 23j],
              [40 + 12j, 98 + 32j, 93 + 37j, 54 + 24j]],

             [[104 + 0j, 247 + 13j, 237 + 23j, 135 + 21j],
              [282 + 30j, 632 + 96j, 604 + 124j, 330 + 86j],
              [246 + 66j, 548 + 180j, 520 + 208j, 282 + 134j],
              [142 + 66j, 307 + 161j, 289 + 179j, 153 + 107j]],

             [[68 + 36j, 157 + 103j, 147 + 113j, 81 + 75j],
              [174 + 138j, 380 + 348j, 352 + 376j, 186 + 230j],
              [138 + 174j, 296 + 432j, 268 + 460j, 138 + 278j],
              [70 + 138j, 145 + 323j, 127 + 341j, 63 + 197j]],

             [[32 + 72j, 68 + 166j, 59 + 175j, 30 + 100j],
              [68 + 192j, 139 + 433j, 117 + 455j, 57 + 255j],
              [38 + 222j, 73 + 499j, 51 + 521j, 21 + 291j],
              [12 + 144j, 20 + 318j, 7 + 331j, 0 + 182j]]])

        assert_array_equal(convolve(small, big, 'full'), out_array)
        assert_array_equal(convolve(big, small, 'full'), out_array)
        assert_array_equal(convolve(small, big, 'same'),
                           out_array[1:3, 1:3, 1:3])
        assert_array_equal(convolve(big, small, 'same'),
                           out_array[0:3, 0:3, 0:3])
        assert_raises(ValueError, convolve, small, big, 'valid')
        assert_array_equal(convolve(big, small, 'valid'),
                           out_array[1:3, 1:3, 1:3])

def filtervertical(H, minleafDomain):
	"""
	This function applies the matched filter on the horizontal 2D histogram.
	Returns filter response for both directions, as stairs can face both ways.
	https://www.youtube.com/watch?v=S7qbelm_4Y8 --> explains matched filter good, couldn't make spectralpython matched filter work
	"""
	# from skimage.feature import canny
	# edges = canny(H)
	# plt.imshow(edges,cmap=plt.cm.gray)
	# plt.show()

	filt = create_matched_filter(minleafDomain)
	# Note that the convolution of the time-reversed wavelet is identical to cross-correlation of the wavelet with the wavelet (autocorrelation) in the input signal --> http://crewes.org/ForOurSponsors/ResearchReports/2002/2002-46.pdf
	
	filty=np.transpose(filt) # transpose to also get stairs in other direction
	fr1=signal.convolve(H,filt, mode='same')

	# plt.subplot(1,3,1)
	# plt.imshow(fr1,cmap='spectral',interpolation='none')
	# plt.title('vert matched filter')
	# plt.colorbar()

	fr2=signal.convolve(H,filty, mode='same')
	# fr3=signal.convolve2d(H,filty, mode='same') # should givve the same result as fr2
	# fr3 = fr3**2 # doesn't really work as there are very high 'outliers' that just take everything away

	# plt.subplot(1,3,2)
	# plt.imshow(fr2,cmap='spectral',interpolation='none')
	# plt.title('vert matched filter transpose')
	# plt.colorbar()

	return fr1, fr2

def gaussian( u, v, size) :
    """Smooths the velocity field with a Gaussian kernel.
    
    Parameters
    ----------
    u : 2d np.ndarray
        the u velocity component field
        
    v : 2d np.ndarray
        the v velocity component field
        
    size : int
        the half width of the kernel. Kernel
        has shape 2*size+1
        
    Returns
    -------
    uf : 2d np.ndarray
        the smoothed u velocity component field
        
    vf : 2d np.ndarray
        the smoothed v velocity component field    
        
    """
    g = _gaussian_kernel( size=size )
    uf = convolve( u, g, mode='same')
    vf = convolve( v, g, mode='same')
    return uf, vf

def computeTensor(im, sigmaG=1, factorSigma=4):

    # returns 3d array of the size of the image
    # yx stores xx, yx, and yy components of the tensor

    # get the luminance of the image, use [0.3, 0.6, 0.1]
    # use numpy's dot
    # blur the image
    imLum = lum(im)
    
    imLumBlurred = zeros( ( height(im), width(im) ) )

    ndimage.filters.gaussian_filter( imLum, sigmaG, 0, imLumBlurred )

    gradX = signal.convolve(imLumBlurred, Sobel, mode='same')
    gradY = signal.convolve(imLumBlurred, transpose(Sobel), mode='same')
    
    # construct 3 2d arrays of the elements of the tensor
    gradXX = gradX*gradX
    gradYY = gradY*gradY
    gradXY = gradX*gradY

    ndimage.filters.gaussian_filter( gradXX, sigmaG * factorSigma, 0, gradXX )
    ndimage.filters.gaussian_filter( gradXY, sigmaG * factorSigma, 0, gradXY )
    ndimage.filters.gaussian_filter( gradYY, sigmaG * factorSigma, 0, gradYY )


    # construct RGB image based on these vals
    out = constantIm(height(im), width(im), 0.0)

    out[:,:,0] = gradXX
    out[:,:,1] = gradXY
    out[:,:,2] = gradYY

    return out

def gaussian_differentiation_kernel(sigma, num_stds, order, delta, scale):
    """
    http://en.wikipedia.org/wiki/Scale_space#Gaussian_derivatives
    :param sigma:
    :param num_stds:
    :param order:
    :param delta:
    :param scale:
    :return:
    """
    delta = list(delta)
    g = gaussian_kernel(sigma, num_stds)
    d = np.ones((1,) * len(order))
    for i in range(len(order)):
        _d = differentiation_kernel(order[i]).astype(float)
        if len(_d) % 2 != 1:
            _d = (np.pad(_d, ((0, 1),), 'constant') + np.pad(_d, ((1, 0),), 'constant')) / 2.
            delta[i] *= 2
        _d /= delta[i] ** order[i]
        _d *= np.sqrt(scale[i]) ** order[i]
        shp = np.ones(len(order), int)
        shp[i] = _d.shape[0]
        _d.shape = shp
        d = spsig.convolve(d, _d)
    d = spsig.convolve(g, d)
    return d