# 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 |
Friday, November 3, 2017
Ok, I'll try to resume the tools fro breaking chaos modulation, starting with spicy
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment