# Author: Travis Oliphant |
| # 1999 -- 2002 |
|
|
| from __future__ import division, print_function, absolute_import |
|
|
| import operator |
| import threading |
| import sys |
| import timeit |
|
|
| from . import sigtools, dlti |
| from ._upfirdn import upfirdn, _output_len |
| from scipy._lib.six import callable |
| from scipy._lib._version import NumpyVersion |
| from scipy import fftpack, linalg |
| from numpy import (allclose, angle, arange, argsort, array, asarray, |
| atleast_1d, atleast_2d, cast, dot, exp, expand_dims, |
| iscomplexobj, mean, ndarray, newaxis, ones, pi, |
| poly, polyadd, polyder, polydiv, polymul, polysub, polyval, |
| product, r_, ravel, real_if_close, reshape, |
| roots, sort, take, transpose, unique, where, zeros, |
| zeros_like) |
| import numpy as np |
| import math |
| from scipy.special import factorial |
| from .windows import get_window |
| from ._arraytools import axis_slice, axis_reverse, odd_ext, even_ext, const_ext |
| from .filter_design import cheby1, _validate_sos |
| from .fir_filter_design import firwin |
|
|
| if sys.version_info.major >= 3 and sys.version_info.minor >= 5: |
| from math import gcd |
| else: |
| from fractions import gcd |
|
|
|
|
| __all__ = ['correlate', 'fftconvolve', 'convolve', 'convolve2d', 'correlate2d', |
| 'order_filter', 'medfilt', 'medfilt2d', 'wiener', 'lfilter', |
| 'lfiltic', 'sosfilt', 'deconvolve', 'hilbert', 'hilbert2', |
| 'cmplx_sort', 'unique_roots', 'invres', 'invresz', 'residue', |
| 'residuez', 'resample', 'resample_poly', 'detrend', |
| 'lfilter_zi', 'sosfilt_zi', 'sosfiltfilt', 'choose_conv_method', |
| 'filtfilt', 'decimate', 'vectorstrength'] |
|
|
|
|
| _modedict = {'valid': 0, 'same': 1, 'full': 2} |
|
|
| _boundarydict = {'fill': 0, 'pad': 0, 'wrap': 2, 'circular': 2, 'symm': 1, |
| 'symmetric': 1, 'reflect': 4} |
|
|
|
|
| _rfft_mt_safe = (NumpyVersion(np.__version__) >= '1.9.0.dev-e24486e') |
|
|
| _rfft_lock = threading.Lock() |
|
|
|
|
| def _valfrommode(mode): |
| try: |
| val = _modedict[mode] |
| except KeyError: |
| if mode not in [0, 1, 2]: |
| raise ValueError("Acceptable mode flags are 'valid' (0)," |
| " 'same' (1), or 'full' (2).") |
| val = mode |
| return val |
|
|
|
|
| def _bvalfromboundary(boundary): |
| try: |
| val = _boundarydict[boundary] << 2 |
| except KeyError: |
| if val not in [0, 1, 2]: |
| raise ValueError("Acceptable boundary flags are 'fill', 'wrap'" |
| " (or 'circular'), \n and 'symm'" |
| " (or 'symmetric').") |
| val = boundary << 2 |
| return val |
|
|
|
|
| def _inputs_swap_needed(mode, shape1, shape2): |
| """ |
| If in 'valid' mode, returns whether or not the input arrays need to be |
| swapped depending on whether `shape1` is at least as large as `shape2` in |
| every dimension. |
| |
| This is important for some of the correlation and convolution |
| implementations in this module, where the larger array input needs to come |
| before the smaller array input when operating in this mode. |
| |
| Note that if the mode provided is not 'valid', False is immediately |
| returned. |
| """ |
| if mode == 'valid': |
| ok1, ok2 = True, True |
|
|
| for d1, d2 in zip(shape1, shape2): |
| if not d1 >= d2: |
| ok1 = False |
| if not d2 >= d1: |
| ok2 = False |
|
|
| if not (ok1 or ok2): |
| raise ValueError("For 'valid' mode, one must be at least " |
| "as large as the other in every dimension") |
|
|
| return not ok1 |
|
|
| return False |
|
|
|
|
| def correlate(in1, in2, mode='full', method='auto'): |
| r""" |
| Cross-correlate two N-dimensional arrays. |
| |
| Cross-correlate `in1` and `in2`, with the output size determined by the |
| `mode` argument. |
| |
| Parameters |
| ---------- |
| in1 : array_like |
| First input. |
| in2 : array_like |
| Second input. Should have the same number of dimensions as `in1`. |
| mode : str {'full', 'valid', 'same'}, optional |
| A string indicating the size of the output: |
| |
| ``full`` |
| The output is the full discrete linear cross-correlation |
| of the inputs. (Default) |
| ``valid`` |
| The output consists only of those elements that do not |
| rely on the zero-padding. In 'valid' mode, either `in1` or `in2` |
| must be at least as large as the other in every dimension. |
| ``same`` |
| The output is the same size as `in1`, centered |
| with respect to the 'full' output. |
| method : str {'auto', 'direct', 'fft'}, optional |
| A string indicating which method to use to calculate the correlation. |
| |
| ``direct`` |
| The correlation is determined directly from sums, the definition of |
| correlation. |
| ``fft`` |
| The Fast Fourier Transform is used to perform the correlation more |
| quickly (only available for numerical arrays.) |
| ``auto`` |
| Automatically chooses direct or Fourier method based on an estimate |
| of which is faster (default). See `convolve` Notes for more detail. |
| |
| .. versionadded:: 0.19.0 |
| |
| Returns |
| ------- |
| correlate : array |
| An N-dimensional array containing a subset of the discrete linear |
| cross-correlation of `in1` with `in2`. |
| |
| See Also |
| -------- |
| choose_conv_method : contains more documentation on `method`. |
| |
| Notes |
| ----- |
| The correlation z of two d-dimensional arrays x and y is defined as:: |
| |
| z[...,k,...] = sum[..., i_l, ...] x[..., i_l,...] * conj(y[..., i_l - k,...]) |
| |
| This way, if x and y are 1-D arrays and ``z = correlate(x, y, 'full')`` then |
| |
| .. math:: |
| |
| z[k] = (x * y)(k - N + 1) |
| = \sum_{l=0}^{||x||-1}x_l y_{l-k+N-1}^{*} |
| |
| for :math:`k = 0, 1, ..., ||x|| + ||y|| - 2` |
| |
| where :math:`||x||` is the length of ``x``, :math:`N = \max(||x||,||y||)`, |
| and :math:`y_m` is 0 when m is outside the range of y. |
| |
| ``method='fft'`` only works for numerical arrays as it relies on |
| `fftconvolve`. In certain cases (i.e., arrays of objects or when |
| rounding integers can lose precision), ``method='direct'`` is always used. |
| |
| Examples |
| -------- |
| Implement a matched filter using cross-correlation, to recover a signal |
| that has passed through a noisy channel. |
| |
| >>> from scipy import signal |
| >>> sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 128) |
| >>> sig_noise = sig + np.random.randn(len(sig)) |
| >>> corr = signal.correlate(sig_noise, np.ones(128), mode='same') / 128 |
| |
| >>> import matplotlib.pyplot as plt |
| >>> clock = np.arange(64, len(sig), 128) |
| >>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True) |
| >>> ax_orig.plot(sig) |
| >>> ax_orig.plot(clock, sig[clock], 'ro') |
| >>> ax_orig.set_title('Original signal') |
| >>> ax_noise.plot(sig_noise) |
| >>> ax_noise.set_title('Signal with noise') |
| >>> ax_corr.plot(corr) |
| >>> ax_corr.plot(clock, corr[clock], 'ro') |
| >>> ax_corr.axhline(0.5, ls=':') |
| >>> ax_corr.set_title('Cross-correlated with rectangular pulse') |
| >>> ax_orig.margins(0, 0.1) |
| >>> fig.tight_layout() |
| >>> fig.show() |
| |
| """ |
| in1 = asarray(in1) |
| in2 = asarray(in2) |
|
|
| if in1.ndim == in2.ndim == 0: |
| return in1 * in2 |
| elif in1.ndim != in2.ndim: |
| raise ValueError("in1 and in2 should have the same dimensionality") |
|
|
| # Don't use _valfrommode, since correlate should not accept numeric modes |
| try: |
| val = _modedict[mode] |
| except KeyError: |
| raise ValueError("Acceptable mode flags are 'valid'," |
| " 'same', or 'full'.") |
|
|
| # this either calls fftconvolve or this function with method=='direct' |
| if method in ('fft', 'auto'): |
| return convolve(in1, _reverse_and_conj(in2), mode, method) |
|
|
| # fastpath to faster numpy.correlate for 1d inputs when possible |
| if _np_conv_ok(in1, in2, mode): |
| return np.correlate(in1, in2, mode) |
|
|
| # _correlateND is far slower when in2.size > in1.size, so swap them |
| # and then undo the effect afterward if mode == 'full'. Also, it fails |
| # with 'valid' mode if in2 is larger than in1, so swap those, too. |
| # Don't swap inputs for 'same' mode, since shape of in1 matters. |
| swapped_inputs = ((mode == 'full') and (in2.size > in1.size) or |
| _inputs_swap_needed(mode, in1.shape, in2.shape)) |
|
|
| if swapped_inputs: |
| in1, in2 = in2, in1 |
|
|
| if mode == 'valid': |
| ps = [i - j + 1 for i, j in zip(in1.shape, in2.shape)] |
| out = np.empty(ps, in1.dtype) |
|
|
| z = sigtools._correlateND(in1, in2, out, val) |
|
|
| else: |
| ps = [i + j - 1 for i, j in zip(in1.shape, in2.shape)] |
|
|
| # zero pad input |
| in1zpadded = np.zeros(ps, in1.dtype) |
| sc = [slice(0, i) for i in in1.shape] |
| in1zpadded[sc] = in1.copy() |
|
|
| if mode == 'full': |
| out = np.empty(ps, in1.dtype) |
| elif mode == 'same': |
| out = np.empty(in1.shape, in1.dtype) |
|
|
| z = sigtools._correlateND(in1zpadded, in2, out, val) |
|
|
| if swapped_inputs: |
| # Reverse and conjugate to undo the effect of swapping inputs |
| z = _reverse_and_conj(z) |
|
|
| return z |
|
|
|
|
| def _centered(arr, newshape): |
| # Return the center newshape portion of the array. |
| newshape = asarray(newshape) |
| currshape = array(arr.shape) |
| startind = (currshape - newshape) // 2 |
| endind = startind + newshape |
| myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] |
| return arr[tuple(myslice)] |
https://github.com/scipy/scipy/blob/master/scipy/signal/signaltools.py |