Python scipy.signal.convolve() Examples

def convolve2d(in1, in2, mode='full'):
    """
        note only support H * W * N * 1 convolve 2d
    """
    in1 = in1.transpose(2, 3, 0, 1) # to N * C * H * W
    in2 = in2.transpose(2, 3, 0, 1)
    out_c, _, kh, kw = in2.shape
    n, _, h, w = in1.shape

    if mode == 'full':
        ph, pw = kh-1, kw-1
        out_h, out_w = h-kh+1+ph*2, w-kw+1+pw*2# TODO
    elif mode == 'valid':
        ph, pw = 0, 0
        out_h, out_w = h-kh+1, w-kw+1 # TODO
    else:
        raise NotImplementedError

    y = cp.empty((n, out_c, out_h, out_w), dtype=in1.dtype)

    col = im2col_gpu(in1, kh, kw, 1, 1, ph, pw)
    y = cp.tensordot(
            col, in2, ((1, 2, 3), (1, 2, 3))).astype(in1.dtype, copy=False)
    y = cp.rollaxis(y, 3, 1)
    return y.transpose(2, 3, 0, 1) 

def test_convolve_generalization():
        ag_convolve = autograd.scipy.signal.convolve
        A_35 = R(3, 5)
        A_34 = R(3, 4)
        A_342 = R(3, 4, 2)
        A_2543 = R(2, 5, 4, 3)
        A_24232 = R(2, 4, 2, 3, 2)

        for mode in ['valid', 'full']:
            assert npo.allclose(ag_convolve(A_35,      A_34, axes=([1], [0]), mode=mode)[1, 2],
                                sp_convolve(A_35[1,:], A_34[:, 2], mode))
            assert npo.allclose(ag_convolve(A_35, A_34, axes=([],[]), dot_axes=([0], [0]), mode=mode),
                                npo.tensordot(A_35, A_34, axes=([0], [0])))
            assert npo.allclose(ag_convolve(A_35, A_342, axes=([1],[2]),
                                            dot_axes=([0], [0]), mode=mode)[2],
                                sum([sp_convolve(A_35[i, :], A_342[i, 2, :], mode)
                                    for i in range(3)]))
            assert npo.allclose(ag_convolve(A_2543, A_24232, axes=([1, 2],[2, 4]),
                                            dot_axes=([0, 3], [0, 3]), mode=mode)[2],
                                sum([sum([sp_convolve(A_2543[i, :, :, j],
                                                    A_24232[i, 2, :, j, :], mode)
                                        for i in range(2)]) for j in range(3)])) 

def sample(self, x_instrument, wl_hi, rdn_hi):
        """Apply instrument sampling to a radiance spectrum, returning predicted measurement."""

        if self.calibration_fixed and all((self.wl_init - wl_hi) < wl_tol):
            return rdn_hi
        wl, fwhm = self.calibration(x_instrument)
        if rdn_hi.ndim == 1:
            return resample_spectrum(rdn_hi, wl_hi, wl, fwhm)
        else:
            resamp = []
            # The "fast resample" option approximates a complete resampling
            # by a convolution with a uniform FWHM.
            if self.fast_resample:
                for i, r in enumerate(rdn_hi):
                    ssrf = spectral_response_function(np.arange(-10, 11), 0, fwhm[0])
                    blur = convolve(r, ssrf, mode='same')
                    resamp.append(interp1d(wl_hi, blur)(wl))
            else:
                for i, r in enumerate(rdn_hi):
                    r2 = resample_spectrum(r, wl_hi, wl, fwhm)
                    resamp.append(r2)
            return np.array(resamp) 

def test_valid_mode2(self):
        # See gh-5897
        a = [1, 2, 3, 6, 5, 3]
        b = [2, 3, 4, 5, 3, 4, 2, 2, 1]
        expected = [70, 78, 73, 65]

        out = convolve(a, b, 'valid')
        assert_array_equal(out, expected)

        out = convolve(b, a, 'valid')
        assert_array_equal(out, expected)

        a = [1 + 5j, 2 - 1j, 3 + 0j]
        b = [2 - 3j, 1 + 0j]
        expected = [2 - 3j, 8 - 10j]

        out = convolve(a, b, 'valid')
        assert_array_equal(out, expected)

        out = convolve(b, a, 'valid')
        assert_array_equal(out, expected) 

def comp_induce_potential(self):
        """
            Compute the induce potential corresponding to the density change
            calculated in get_spatial_density
        """

        from scipy.signal import convolve

        Nx, Ny, Nz = self.mesh[0].size, self.mesh[1].size, self.mesh[2].size

        grid = np.zeros((Nx, Ny, Nz), dtype = np.float64)
        factor = self.dr[0]*self.dr[1]*self.dr[2]/(np.sqrt(2*np.pi)**3)

        libnao.comp_spatial_grid_pot(
            self.dr.ctypes.data_as(POINTER(c_double)), 
            self.mesh[0].ctypes.data_as(POINTER(c_double)),
            self.mesh[1].ctypes.data_as(POINTER(c_double)),
            self.mesh[2].ctypes.data_as(POINTER(c_double)),
            grid.ctypes.data_as(POINTER(c_double)),
            c_int(Nx), c_int(Ny), c_int(Nz))

        return convolve(grid, self.dn_spatial, mode="same", method="fft")*factor 

def convolve2d(in1, in2, mode='full'):
    """
        note only support H * W * N * 1 convolve 2d
    """
    in1 = in1.transpose(2, 3, 0, 1) # to N * C * H * W
    in2 = in2.transpose(2, 3, 0, 1)
    out_c, _, kh, kw = in2.shape
    n, _, h, w = in1.shape

    if mode == 'full':
        ph, pw = kh-1, kw-1
        out_h, out_w = h-kh+1+ph*2, w-kw+1+pw*2# TODO
    elif mode == 'valid':
        ph, pw = 0, 0
        out_h, out_w = h-kh+1, w-kw+1 # TODO
    else:
        raise NotImplementedError

    y = cp.empty((n, out_c, out_h, out_w), dtype=in1.dtype)

    col = im2col_gpu(in1, kh, kw, 1, 1, ph, pw)
    y = cp.tensordot(
            col, in2, ((1, 2, 3), (1, 2, 3))).astype(in1.dtype, copy=False)
    y = cp.rollaxis(y, 3, 1)
    return y.transpose(2, 3, 0, 1) 

def ref_depthwise_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[0]):
            __y = []
            for k in range(w.shape[1]):
                __y.append(signal.convolve(x[i, j], w[j, k], mode=padding))
            _y.append(np.stack(__y, axis=0))
        y.append(np.concatenate(_y, axis=0))
    y = np.array(y)
    return y 

def ref_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[1]):
            __y = []
            for k in range(w.shape[0]):
                __y.append(signal.convolve(x[i, k], w[k, j], mode=padding))
            _y.append(np.sum(np.stack(__y, axis=-1), axis=-1))
        y.append(_y)
    y = np.array(y)
    return y 

def ref_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[1]):
            __y = []
            for k in range(w.shape[0]):
                __y.append(signal.convolve(x[i, k], w[k, j], mode=padding))
            _y.append(np.sum(np.stack(__y, axis=-1), axis=-1))
        y.append(_y)
    y = np.array(y)
    return y 

def fftconvolve_old(in1, in2, in3=None, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    copied from scipy.signal.signaltools, but here used to try out inverse filter
    doesn't work or I can't get it to work

     2010-10-23:
    looks ok to me for 1d,
    from results below with padded data array (fftp)
    but it doesn't work for multidimensional inverse filter (fftn)
    original signal.fftconvolve also uses fftn

    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1+s2-1

    # Always use 2**n-sized FFT
    fsize = 2**np.ceil(np.log2(size))
    IN1 = fftn(in1,fsize)
    #IN1 *= fftn(in2,fsize) #JP: this looks like the only change I made
    IN1 /= fftn(in2,fsize)  # use inverse filter
    # note the inverse is elementwise not matrix inverse
    # is this correct, NO  doesn't seem to work for VARMA
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if product(s1,axis=0) > product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret,osize)
    elif mode == "valid":
        return _centered(ret,abs(s2-s1)+1) 

def temporally_smooth(landmarks):
        """ apply temporal filtering on the 2D points """
        logger.debug("Temporally Smooth")
        filter_half_length = 2
        temporal_filter = np.ones((1, 1, 2 * filter_half_length + 1))
        temporal_filter = temporal_filter / temporal_filter.sum()

        start_tileblock = np.tile(landmarks[:, :, 0][:, :, np.newaxis], [1, 1, filter_half_length])
        end_tileblock = np.tile(landmarks[:, :, -1][:, :, np.newaxis], [1, 1, filter_half_length])
        landmarks_padded = np.dstack((start_tileblock, landmarks, end_tileblock))

        retval = signal.convolve(landmarks_padded, temporal_filter, mode='valid', method='fft')
        logger.debug("Temporally Smoothed: %s", retval)
        return retval 

def dmeas_dinstrumentb(self, x_instrument, wl_hi, rdn_hi):
        """Jacobian of radiance with respect to the instrument parameters
        that are unknown and not retrieved, i.e., the inevitable persisting
        uncertainties in instrument spectral and radiometric calibration.

        Input: meas, a vector of size n_chan
        Returns: Kb_instrument, a matrix of size [n_measurements x nb_instrument]
        """

        # Uncertainty due to radiometric calibration
        meas = self.sample(x_instrument, wl_hi, rdn_hi)
        dmeas_dinstrument = np.hstack(
            (np.diagflat(meas), np.zeros((self.n_chan, 2))))

        # Uncertainty due to spectral calibration
        if self.bval[-2] > 1e-6:
            dmeas_dinstrument[:, -2] = self.sample(x_instrument, wl_hi,
                                                   np.hstack((np.diff(rdn_hi), np.array([0]))))

        # Uncertainty due to spectral stray light
        if self.bval[-1] > 1e-6:
            ssrf = spectral_response_function(np.arange(-10, 11), 0, 4)
            blur = convolve(meas, ssrf, mode='same')
            dmeas_dinstrument[:, -1] = blur - meas

        return dmeas_dinstrument 

def convolve_numpy_array(self, input, gaussian_filter):
        from scipy.signal import convolve
        return convolve(input, gaussian_filter, mode='same') 

def smooth_ppg_signal(signal, sample_rate=PPG_SAMPLE_RATE, numtaps=PPG_FIR_FILTER_TAP_NUM, cutoff=PPG_FILTER_CUTOFF):
    if numtaps % 2 == 0:
        numtaps += 1
    return convolve(signal, firwin(numtaps, [x*2/sample_rate for x in cutoff], pass_zero=False), mode='valid').tolist() 

def ref_depthwise_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[0]):
            __y = []
            for k in range(w.shape[1]):
                __y.append(signal.convolve(x[i, j], w[j, k], mode=padding))
            _y.append(np.stack(__y, axis=0))
        y.append(np.concatenate(_y, axis=0))
    y = np.array(y)
    return y 

def ref_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[1]):
            __y = []
            for k in range(w.shape[0]):
                __y.append(signal.convolve(x[i, k], w[k, j], mode=padding))
            _y.append(np.sum(np.stack(__y, axis=-1), axis=-1))
        y.append(_y)
    y = np.array(y)
    return y 

def ref_depthwise_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[0]):
            __y = []
            for k in range(w.shape[1]):
                __y.append(signal.convolve(x[i, j], w[j, k], mode=padding))
            _y.append(np.stack(__y, axis=0))
        y.append(np.concatenate(_y, axis=0))
    y = np.array(y)
    return y 

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_array_equal(convolve(small, big, 'valid'),
                           out_array[1:3, 1:3, 1:3])
        assert_array_equal(convolve(big, small, 'valid'),
                           out_array[1:3, 1:3, 1:3]) 

def ref_conv(x, w, padding, data_format):
    y = []
    for i in range(x.shape[0]):
        _y = []
        for j in range(w.shape[1]):
            __y = []
            for k in range(w.shape[0]):
                __y.append(signal.convolve(x[i, k], w[k, j], mode=padding))
            _y.append(np.sum(np.stack(__y, axis=-1), axis=-1))
        y.append(_y)
    y = np.array(y)
    return y 

def cascade_filters(b1,a1,b2,a2):
    """
    Cascade two IIR digital filters into a single (b,a) coefficient set.
    
    To cascade two digital filters (system functions) given their numerator
    and denominator coefficients you simply convolve the coefficient arrays.
    
    Parameters
    ----------
    b1 : ndarray of numerator coefficients for filter 1
    a1 : ndarray of denominator coefficients for filter 1
    b2 : ndarray of numerator coefficients for filter 2
    a2 : ndarray of denominator coefficients for filter 2
    
    Returns
    -------
    b : ndarray of numerator coefficients for the cascade
    a : ndarray of denominator coefficients for the cascade
    
    Examples
    --------
    >>> from scipy import signal
    >>> b1,a1 = signal.butter(3, 0.1)
    >>> b2,a2 = signal.butter(3, 0.15)
    >>> b,a = cascade_filters(b1,a1,b2,a2)
    """
    return signal.convolve(b1,b2), signal.convolve(a1,a2) 

def CIC(M, K):
    """
    A functional form implementation of a cascade of integrator comb (CIC) filters.

    Parameters
    ----------
    M : Effective number of taps per section (typically the decimation factor).
    K : The number of CIC sections cascaded (larger K gives the filter a wider image rejection bandwidth).

    Returns
    -------
    b : FIR filter coefficients for a simple direct form implementation using the filter() function.

    Notes
    -----
    Commonly used in multirate signal processing digital down-converters and digital up-converters. A true CIC filter
    requires no multiplies, only add and subtract operations. The functional form created here is a simple FIR requiring
    real coefficient multiplies via filter().

    Mark Wickert July 2013
    """

    if K == 1:
        b = np.ones(M)
    else:
        h = np.ones(M)
        b = h
        for i in range(1, K):
            b = signal.convolve(b, h)  # cascade by convolving impulse responses

    # Make filter have unity gain at DC
    return b / np.sum(b) 

def save_wav(wav, path, hparams):
	wav = wav / np.abs(wav).max() * 0.999
	f1 = 0.5 * 32767 / max(0.01, np.max(np.abs(wav)))
	f2 = np.sign(wav) * np.power(np.abs(wav), 0.95)
	wav = f1 * f2
	wav = signal.convolve(wav, signal.firwin(hparams.num_freq, [hparams.fmin, hparams.fmax], pass_zero=False, fs=hparams.sample_rate))
	#proposed by @dsmiller
	wavfile.write(path, hparams.sample_rate, wav.astype(np.int16)) 

def lsf2poly(L):
    order = len(L)
    Q = L[::2]
    P = L[1::2]
    poles_P = np.r_[np.exp(1j*P),np.exp(-1j*P)]
    poles_Q = np.r_[np.exp(1j*Q),np.exp(-1j*Q)]
    
    P = np.poly(poles_P)
    Q = np.poly(poles_Q)
    
    P = convolve(P, np.array([1.0, -1.0]))
    Q = convolve(Q, np.array([1.0, 1.0]))
    
    a = 0.5*(P+Q)
    return a[:-1] 

def test_poly_vs_filtfilt(self):
        # Check that up=1.0 gives same answer as filtfilt + slicing
        random_state = np.random.RandomState(17)
        try_types = (int, np.float32, np.complex64, float, complex)
        size = 10000
        down_factors = [2, 11, 79]

        for dtype in try_types:
            x = random_state.randn(size).astype(dtype)
            if dtype in (np.complex64, np.complex128):
                x += 1j * random_state.randn(size)

            # resample_poly assumes zeros outside of signl, whereas filtfilt
            # can only constant-pad. Make them equivalent:
            x[0] = 0
            x[-1] = 0

            for down in down_factors:
                h = signal.firwin(31, 1. / down, window='hamming')
                yf = filtfilt(h, 1.0, x, padtype='constant')[::down]

                # Need to pass convolved version of filter to resample_poly,
                # since filtfilt does forward and backward, but resample_poly
                # only goes forward
                hc = convolve(h, h[::-1])
                y = signal.resample_poly(x, 1, down, window=hc)
                assert_allclose(yf, y, atol=1e-7, rtol=1e-7) 

def test_random_data(self):
        np.random.seed(1234)
        a = np.random.rand(1233) + 1j * np.random.rand(1233)
        b = np.random.rand(1321) + 1j * np.random.rand(1321)
        c = fftconvolve(a, b, 'full')
        d = np.convolve(a, b, 'full')
        assert_(np.allclose(c, d, rtol=1e-10)) 

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 test_mismatched_dims(self):
        # Input arrays should have the same number of dimensions
        assert_raises(ValueError, convolve, [1], 2, method='direct')
        assert_raises(ValueError, convolve, 1, [2], method='direct')
        assert_raises(ValueError, convolve, [1], 2, method='fft')
        assert_raises(ValueError, convolve, 1, [2], method='fft')
        assert_raises(ValueError, convolve, [1], [[2]])
        assert_raises(ValueError, convolve, [3], 2) 

def test_convolve_method_large_input(self):
        # This is really a test that convolving two large integers goes to the
        # direct method even if they're in the fft method.
        for n in [10, 20, 50, 51, 52, 53, 54, 60, 62]:
            z = np.array([2**n], dtype=np.int64)
            fft = convolve(z, z, method='fft')
            direct = convolve(z, z, method='direct')

            # this is the case when integer precision gets to us
            # issue #6076 has more detail, hopefully more tests after resolved
            if n < 50:
                assert_equal(fft, direct)
                assert_equal(fft, 2**(2*n))
                assert_equal(direct, 2**(2*n)) 

def test_invalid_params(self):
        a = [3, 4, 5]
        b = [1, 2, 3]
        assert_raises(ValueError, convolve, a, b, mode='spam')
        assert_raises(ValueError, convolve, a, b, mode='eggs', method='fft')
        assert_raises(ValueError, convolve, a, b, mode='ham', method='direct')
        assert_raises(ValueError, convolve, a, b, mode='full', method='bacon')
        assert_raises(ValueError, convolve, a, b, mode='same', method='bacon') 

def test_same_mode(self):
        a = [1, 2, 3, 3, 1, 2]
        b = [1, 4, 3, 4, 5, 6, 7, 4, 3, 2, 1, 1, 3]
        c = convolve(a, b, 'same')
        d = array([57, 61, 63, 57, 45, 36])
        assert_array_equal(c, d)