Python scipy.signal.correlate() Examples

def filter_and_smooth_angular_velocity(angular_velocity,
                                       low_pass_kernel_size, clip_percentile, plot=False):
  """Reduce the noise in a velocity signal."""

  max_value = np.percentile(angular_velocity, clip_percentile)
  print("Clipping angular velocity norms to {} rad/s ...".format(max_value))
  angular_velocity_clipped = np.clip(angular_velocity, -max_value, max_value)
  print("Done clipping angular velocity norms...")

  low_pass_kernel = np.ones((low_pass_kernel_size, 1)) / low_pass_kernel_size
  print("Smoothing with kernel size {} samples...".format(low_pass_kernel_size))

  angular_velocity_smoothed = signal.correlate(angular_velocity_clipped,
                                               low_pass_kernel, 'same')

  print("Done smoothing angular velocity norms...")

  if plot:
    plot_angular_velocities("Angular Velocities", angular_velocity,
                            angular_velocity_smoothed, True)

  return angular_velocity_smoothed.copy() 

def offset_compensation(x1, x2, ns, ndec, nlag=2000):
    '''Find and correct a constant time offset between two signals using a 
    cross-correlation
    
    Parameters:
        s1, s2:     Arrays containing the input signals
        ns:         Number of samples to use for cross-correlation
        ndec:       Decimation factor prior to cross-correlation
        nlag:       Number of lag bins for cross-correlation
    Returns:
        x2s:        The signal x2 time-shifted so that it aligns with x1. Edges 
                    are padded with zeros.         
    '''
    s1 = x1[0:int(ns)]
    s2 = x2[0:int(ns)]

    # cross-correlate to find the offset
    os = find_channel_offset(s1, s2, ndec, nlag)

    if(os == 0):
        return x2
    else:
        return shift(x2, os) 

def gabor_feature_single_job(a, filters, fm_i, label, cluster_center_number, save_flag):
    # convolution
    start_time = time.time()
    # b=SN.correlate(a,filters[i]) # too slow
    b = signal.correlate(a, filters[fm_i], mode='same')
    end_time = time.time()
    print('feature %d done (%f s)' % (fm_i, end_time - start_time))

    # show Gabor filter output
    if save_flag:
        img = (b[:, :, int(a.shape[2] / 2)]).copy()
        plt.imsave('./result/gabor_output(%d).png' % fm_i, img, cmap='gray')  # save fig

    # generate feature vector
    start_time = time.time()
    result = generate_feature_vector(b=b, label=label, cluster_center_number=cluster_center_number)
    end_time = time.time()
    print('feature vector %d done (%f s)' % (fm_i, end_time - start_time))
    return fm_i, result 

def calculateBackgroundSignal(self, mat, vmat, nuc_cov):
        offset=self.start-mat.start-vmat.w
        if offset<0:
            raise Exception("Insufficient flanking region on \
                    mat to calculate signal")
        self.vmat = vmat
        self.bias_mat = mat
        self.cov = CoverageTrack(self.chrom, self.start, self.end)
        self.cov.calculateCoverage(self.bias_mat, vmat.lower,
                                   vmat.upper, vmat.w*2+1)
        self.nuc_cov = nuc_cov.vals
        self.vals = signal.correlate(self.bias_mat.get(vmat.lower,vmat.upper,
                                                         self.bias_mat.start + offset,
                                                         self.bias_mat.end - offset),
                                       vmat.mat,mode = 'valid')[0]
        self.vals = self.vals * self.nuc_cov/ self.cov.vals 

def test_rank0(self, dt):
        a = np.array(np.random.randn()).astype(dt)
        a += 1j * np.array(np.random.randn()).astype(dt)
        b = np.array(np.random.randn()).astype(dt)
        b += 1j * np.array(np.random.randn()).astype(dt)

        y_r = (correlate(a.real, b.real)
               + correlate(a.imag, b.imag)).astype(dt)
        y_r += 1j * (-correlate(a.real, b.imag) + correlate(a.imag, b.real))

        y = correlate(a, b, 'full')
        assert_array_almost_equal(y, y_r, decimal=self.decimal(dt) - 1)
        assert_equal(y.dtype, dt)

        assert_equal(correlate([1], [2j]), correlate(1, 2j))
        assert_equal(correlate([2j], [3j]), correlate(2j, 3j))
        assert_equal(correlate([3j], [4]), correlate(3j, 4)) 

def test_consistency_correlate_funcs(self):
        # Compare np.correlate, signal.correlate, signal.correlate2d
        a = np.arange(5)
        b = np.array([3.2, 1.4, 3])
        for mode in ['full', 'valid', 'same']:
            assert_almost_equal(np.correlate(a, b, mode=mode),
                                signal.correlate(a, b, mode=mode))
            assert_almost_equal(np.squeeze(signal.correlate2d([a], [b],
                                                              mode=mode)),
                                signal.correlate(a, b, mode=mode))

            # See gh-5897
            if mode == 'valid':
                assert_almost_equal(np.correlate(b, a, mode=mode),
                                    signal.correlate(b, a, mode=mode))
                assert_almost_equal(np.squeeze(signal.correlate2d([b], [a],
                                                                  mode=mode)),
                                    signal.correlate(b, a, mode=mode)) 

def test_pdos_1d():
    pad=lambda x: pad_zeros(x, nadd=len(x)-1)
    n=500; w=welch(n)
    # 1 second signal
    t=np.linspace(0,1,n); dt=t[1]-t[0]
    # sum of sin()s with random freq and phase shift, 10 frequencies from
    # f=0...100 Hz
    v=np.array([np.sin(2*np.pi*f*t + rand()*2*np.pi) for f in rand(10)*100]).sum(0)
    f=np.fft.fftfreq(2*n-1, dt)[:n]

    c1=mirror(ifft(abs(fft(pad(v)))**2.0)[:n].real)
    c2=correlate(v,v,'full')
    c3=mirror(acorr(v,norm=False))
    assert np.allclose(c1, c2)
    assert np.allclose(c1, c3)

    p1=(abs(fft(pad(v)))**2.0)[:n]
    p2=(abs(fft(mirror(acorr(v,norm=False)))))[:n]
    assert np.allclose(p1, p2)

    p1=(abs(fft(pad(v*w)))**2.0)[:n]
    p2=(abs(fft(mirror(acorr(v*w,norm=False)))))[:n]
    assert np.allclose(p1, p2) 

def calculateSignal(self, mat, vmat):
        offset=self.start-mat.start-vmat.w
        if offset<0:
            raise Exception("Insufficient flanking region on \
                    mat to calculate signal")
        self.vals = signal.correlate(mat.get(vmat.lower,vmat.upper,
                                              mat.start + offset, mat.end - offset),
                                       vmat.mat,mode = 'valid')[0] 

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

def _rmatvec(self, x):
        x = np.reshape(x, self.dims)
        y = correlate(x, self.h, mode='same', method=self.method)
        y = y.ravel()
        return y 

def _setup_rank1(self, dt, mode):
        np.random.seed(9)
        a = np.random.randn(10).astype(dt)
        a += 1j * np.random.randn(10).astype(dt)
        b = np.random.randn(8).astype(dt)
        b += 1j * np.random.randn(8).astype(dt)

        y_r = (correlate(a.real, b.real, mode=mode) +
               correlate(a.imag, b.imag, mode=mode)).astype(dt)
        y_r += 1j * (-correlate(a.real, b.imag, mode=mode) +
                     correlate(a.imag, b.real, mode=mode))
        return a, b, y_r 

def test_rank1_valid(self, dt):
        a, b, y_r = self._setup_rank1(dt, 'valid')
        y = correlate(a, b, 'valid')
        assert_array_almost_equal(y, y_r, decimal=self.decimal(dt))
        assert_equal(y.dtype, dt)

        # See gh-5897
        y = correlate(b, a, 'valid')
        assert_array_almost_equal(y, y_r[::-1].conj(), decimal=self.decimal(dt))
        assert_equal(y.dtype, dt) 

def test_rank1_same(self, dt):
        a, b, y_r = self._setup_rank1(dt, 'same')
        y = correlate(a, b, 'same')
        assert_array_almost_equal(y, y_r, decimal=self.decimal(dt))
        assert_equal(y.dtype, dt) 

def test_rank1_full(self, dt):
        a, b, y_r = self._setup_rank1(dt, 'full')
        y = correlate(a, b, 'full')
        assert_array_almost_equal(y, y_r, decimal=self.decimal(dt))
        assert_equal(y.dtype, dt) 

def test_swap_full(self, dt):
        d = np.array([0.+0.j, 1.+1.j, 2.+2.j], dtype=dt)
        k = np.array([1.+3.j, 2.+4.j, 3.+5.j, 4.+6.j], dtype=dt)
        y = correlate(d, k)
        assert_equal(y, [0.+0.j, 10.-2.j, 28.-6.j, 22.-6.j, 16.-6.j, 8.-4.j]) 

def test_rank3(self, dt):
        a = np.random.randn(10, 8, 6).astype(dt)
        a += 1j * np.random.randn(10, 8, 6).astype(dt)
        b = np.random.randn(8, 6, 4).astype(dt)
        b += 1j * np.random.randn(8, 6, 4).astype(dt)

        y_r = (correlate(a.real, b.real)
               + correlate(a.imag, b.imag)).astype(dt)
        y_r += 1j * (-correlate(a.real, b.imag) + correlate(a.imag, b.real))

        y = correlate(a, b, 'full')
        assert_array_almost_equal(y, y_r, decimal=self.decimal(dt) - 1)
        assert_equal(y.dtype, dt) 

def computeBias(self, fasta, chromDict, pwm):
        """compute bias track based on sequence and pwm"""
        self.slop(chromDict, up = pwm.up, down = pwm.down)
        sequence = seq.get_sequence(self, fasta)
        seqmat = seq.seq_to_mat(sequence, pwm.nucleotides)
        self.vals = signal.correlate(seqmat,np.log(pwm.mat),mode='valid')[0]
        self.start += pwm.up
        self.end -= pwm.down 

def test_mismatched_dims(self):
        # Input arrays should have the same number of dimensions
        assert_raises(ValueError, correlate, [1], 2, method='direct')
        assert_raises(ValueError, correlate, 1, [2], method='direct')
        assert_raises(ValueError, correlate, [1], 2, method='fft')
        assert_raises(ValueError, correlate, 1, [2], method='fft')
        assert_raises(ValueError, correlate, [1], [[2]])
        assert_raises(ValueError, correlate, [3], 2) 

def xcorr(s1, s2, nlead, nlag):
    ''' cross-correlated s1 and s2'''
    return signal.correlate(s1, np.pad(s2, (nlag, nlead), mode='constant'),
        mode='valid') 

def find_channel_offset(s1, s2, nd, nl):
    '''use cross-correlation to find channel offset in samples'''
    B1 = signal.decimate(s1, nd)
    B2 = np.pad(signal.decimate(s2, nd), (nl, nl), 'constant')
    xc = np.abs(signal.correlate(B1, B2, mode='valid'))
    return (np.argmax(xc) - nl)*nd 

def calculate_time_offset_from_signals(times_A, signal_A,
                                       times_B, signal_B,
                                       plot=False, block=True):
  """ Calculates the time offset between signal A and signal B. """
  convoluted_signals = signal.correlate(signal_B, signal_A)
  dt_A = np.mean(np.diff(times_A))
  offset_indices = np.arange(-len(signal_A) + 1, len(signal_B))
  max_index = np.argmax(convoluted_signals)
  offset_index = offset_indices[max_index]
  time_offset = dt_A * offset_index + times_B[0] - times_A[0]
  if plot:
    plot_results(times_A, times_B, signal_A, signal_B, convoluted_signals,
                 time_offset, block=block)
  return time_offset 

def test__auto_correlation():
    n_time_samples, n_tapers = 100, 3
    test_data = np.random.rand(n_tapers, n_time_samples)
    rxx = _auto_correlation(test_data)[:, :n_time_samples]

    for taper_ind in np.arange(n_tapers):
        expected_correlation = correlate(
            test_data[taper_ind, :], test_data[taper_ind, :])[
                n_time_samples - 1:]
        assert np.allclose(rxx[taper_ind], expected_correlation) 

def __xcorr_normalization(size, dtype):
    arr = np.ones(shape=(size,), dtype=dtype)
    return correlate(arr, arr, mode="full") 

def _convolve_data_adjoint(output, filt, data_shape,
                           mode='full', strides=None,
                           multi_channel=False):
    D, b, B, m, n, s, c_i, c_o, p = _get_convolve_params(
        data_shape, filt.shape,
        mode, strides, multi_channel)

    # Normalize shapes.
    output = output.reshape((B, c_o) + p)
    filt = filt.reshape((c_o, c_i) + n)
    data = np.zeros((B, c_i) + m, dtype=output.dtype)
    slc = tuple(slice(None, None, s_d) for s_d in s)
    if mode == 'full':
        output_kj = np.zeros([m_d + n_d - 1 for m_d, n_d in zip(m, n)],
                             dtype=output.dtype)
        adjoint_mode = 'valid'
    elif mode == 'valid':
        output_kj = np.zeros([max(m_d, n_d) - min(m_d, n_d) + 1
                              for m_d, n_d in zip(m, n)],
                             dtype=output.dtype)
        if all(m_d >= n_d for m_d, n_d in zip(m, n)):
            adjoint_mode = 'full'
        else:
            adjoint_mode = 'valid'

    for k in range(B):
        for j in range(c_o):
            for i in range(c_i):
                output_kj[slc] = output[k, j]
                data[k, i] += signal.correlate(
                    output_kj, filt[j, i],
                    mode=adjoint_mode)

    # Reshape.
    data = data.reshape(data_shape)
    return data 

def _convolve_filter_adjoint(output, data, filt_shape,
                             mode='full', strides=None,
                             multi_channel=False):
    D, b, B, m, n, s, c_i, c_o, p = _get_convolve_params(
        data.shape, filt_shape,
        mode, strides, multi_channel)

    # Normalize shapes.
    data = data.reshape((B, c_i) + m)
    output = output.reshape((B, c_o) + p)
    slc = tuple(slice(None, None, s_d) for s_d in s)
    if mode == 'full':
        output_kj = np.zeros([m_d + n_d - 1 for m_d, n_d in zip(m, n)],
                             dtype=output.dtype)
        adjoint_mode = 'valid'
    elif mode == 'valid':
        output_kj = np.zeros([max(m_d, n_d) - min(m_d, n_d) + 1
                              for m_d, n_d in zip(m, n)],
                             dtype=output.dtype)
        if all(m_d >= n_d for m_d, n_d in zip(m, n)):
            adjoint_mode = 'valid'
        else:
            adjoint_mode = 'full'

    filt = np.zeros((c_o, c_i) + n, dtype=output.dtype)
    for k in range(B):
        for j in range(c_o):
            for i in range(c_i):
                output_kj[slc] = output[k, j]
                filt[j, i] += signal.correlate(
                    output_kj, data[k, i], mode=adjoint_mode)

    # Reshape.
    filt = filt.reshape(filt_shape)
    return filt 

def _setup_rank1(self, mode):
        np.random.seed(9)
        a = np.random.randn(10).astype(self.dt)
        a += 1j * np.random.randn(10).astype(self.dt)
        b = np.random.randn(8).astype(self.dt)
        b += 1j * np.random.randn(8).astype(self.dt)

        y_r = (correlate(a.real, b.real, mode=mode) +
               correlate(a.imag, b.imag, mode=mode)).astype(self.dt)
        y_r += 1j * (-correlate(a.real, b.imag, mode=mode) +
                correlate(a.imag, b.real, mode=mode))
        return a, b, y_r 

def estimate_time_shift(x, y):
    """ Computes the cross-correlation between time series x and y, grabs the
        index of where it's a maximum.  This yields the time difference in
        samples between x and y.
    """
    if DEBUG: print("computing cross-correlation")
    corr = signal.correlate(y, x, mode='same', method='fft')
    if DEBUG: print("finished computing cross-correlation")

    nx, ny = len(x), len(y)
    t_samples = np.arange(nx)
    ct_samples = t_samples - nx//2  # try to center time shift (x axis) on zero
    cmax_ind = np.argmax(corr)      # where is the max of the cross-correlation?
    dt = ct_samples[cmax_ind]       # grab the time shift value corresponding to the max c-corr

    if DEBUG:
        print("cmax_ind, nx//2, ny//2, dt =",cmax_ind, nx//2, ny//2, dt)
        fig, (ax_x, ax_y, ax_corr) = plt.subplots(3, 1)
        ax_x.get_shared_x_axes().join(ax_x, ax_y)
        ax_x.plot(t_samples, x)
        ax_y.plot(t_samples, y)
        ax_corr.plot(ct_samples, corr)
        plt.show()
        
    return dt


#  for use in filtering filenames 

def test_rank1_valid(self):
        a, b, y_r = self._setup_rank1()
        y = correlate(a, b, 'valid')
        assert_array_almost_equal(y, y_r[1:4])
        self.assertTrue(y.dtype == self.dt) 

def test_rank1_same(self):
        a, b, y_r = self._setup_rank1()
        y = correlate(a, b, 'same')
        assert_array_almost_equal(y, y_r[:-1])
        self.assertTrue(y.dtype == self.dt) 

def test_rank1_full(self):
        a, b, y_r = self._setup_rank1()
        y = correlate(a, b, 'full')
        assert_array_almost_equal(y, y_r)
        self.assertTrue(y.dtype == self.dt)