1# Author: Travis Oliphant 2# 1999 -- 2002 3 4import operator 5import math 6import timeit 7from scipy.spatial import cKDTree 8from . import sigtools, dlti 9from ._upfirdn import upfirdn, _output_len, _upfirdn_modes 10from scipy import linalg, fft as sp_fft 11from scipy.fft._helper import _init_nd_shape_and_axes 12from scipy._lib._util import prod as _prod 13import numpy as np 14from scipy.special import lambertw 15from .windows import get_window 16from ._arraytools import axis_slice, axis_reverse, odd_ext, even_ext, const_ext 17from .filter_design import cheby1, _validate_sos 18from .fir_filter_design import firwin 19from ._sosfilt import _sosfilt 20import warnings 21 22 23__all__ = ['correlate', 'correlation_lags', 'correlate2d', 24 'convolve', 'convolve2d', 'fftconvolve', 'oaconvolve', 25 'order_filter', 'medfilt', 'medfilt2d', 'wiener', 'lfilter', 26 'lfiltic', 'sosfilt', 'deconvolve', 'hilbert', 'hilbert2', 27 'cmplx_sort', 'unique_roots', 'invres', 'invresz', 'residue', 28 'residuez', 'resample', 'resample_poly', 'detrend', 29 'lfilter_zi', 'sosfilt_zi', 'sosfiltfilt', 'choose_conv_method', 30 'filtfilt', 'decimate', 'vectorstrength'] 31 32 33_modedict = {'valid': 0, 'same': 1, 'full': 2} 34 35_boundarydict = {'fill': 0, 'pad': 0, 'wrap': 2, 'circular': 2, 'symm': 1, 36 'symmetric': 1, 'reflect': 4} 37 38 39def _valfrommode(mode): 40 try: 41 return _modedict[mode] 42 except KeyError as e: 43 raise ValueError("Acceptable mode flags are 'valid'," 44 " 'same', or 'full'.") from e 45 46 47def _bvalfromboundary(boundary): 48 try: 49 return _boundarydict[boundary] << 2 50 except KeyError as e: 51 raise ValueError("Acceptable boundary flags are 'fill', 'circular' " 52 "(or 'wrap'), and 'symmetric' (or 'symm').") from e 53 54 55def _inputs_swap_needed(mode, shape1, shape2, axes=None): 56 """Determine if inputs arrays need to be swapped in `"valid"` mode. 57 58 If in `"valid"` mode, returns whether or not the input arrays need to be 59 swapped depending on whether `shape1` is at least as large as `shape2` in 60 every calculated dimension. 61 62 This is important for some of the correlation and convolution 63 implementations in this module, where the larger array input needs to come 64 before the smaller array input when operating in this mode. 65 66 Note that if the mode provided is not 'valid', False is immediately 67 returned. 68 69 """ 70 if mode != 'valid': 71 return False 72 73 if not shape1: 74 return False 75 76 if axes is None: 77 axes = range(len(shape1)) 78 79 ok1 = all(shape1[i] >= shape2[i] for i in axes) 80 ok2 = all(shape2[i] >= shape1[i] for i in axes) 81 82 if not (ok1 or ok2): 83 raise ValueError("For 'valid' mode, one must be at least " 84 "as large as the other in every dimension") 85 86 return not ok1 87 88 89def correlate(in1, in2, mode='full', method='auto'): 90 r""" 91 Cross-correlate two N-dimensional arrays. 92 93 Cross-correlate `in1` and `in2`, with the output size determined by the 94 `mode` argument. 95 96 Parameters 97 ---------- 98 in1 : array_like 99 First input. 100 in2 : array_like 101 Second input. Should have the same number of dimensions as `in1`. 102 mode : str {'full', 'valid', 'same'}, optional 103 A string indicating the size of the output: 104 105 ``full`` 106 The output is the full discrete linear cross-correlation 107 of the inputs. (Default) 108 ``valid`` 109 The output consists only of those elements that do not 110 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 111 must be at least as large as the other in every dimension. 112 ``same`` 113 The output is the same size as `in1`, centered 114 with respect to the 'full' output. 115 method : str {'auto', 'direct', 'fft'}, optional 116 A string indicating which method to use to calculate the correlation. 117 118 ``direct`` 119 The correlation is determined directly from sums, the definition of 120 correlation. 121 ``fft`` 122 The Fast Fourier Transform is used to perform the correlation more 123 quickly (only available for numerical arrays.) 124 ``auto`` 125 Automatically chooses direct or Fourier method based on an estimate 126 of which is faster (default). See `convolve` Notes for more detail. 127 128 .. versionadded:: 0.19.0 129 130 Returns 131 ------- 132 correlate : array 133 An N-dimensional array containing a subset of the discrete linear 134 cross-correlation of `in1` with `in2`. 135 136 See Also 137 -------- 138 choose_conv_method : contains more documentation on `method`. 139 correlation_lags : calculates the lag / displacement indices array for 1D 140 cross-correlation. 141 142 Notes 143 ----- 144 The correlation z of two d-dimensional arrays x and y is defined as:: 145 146 z[...,k,...] = sum[..., i_l, ...] x[..., i_l,...] * conj(y[..., i_l - k,...]) 147 148 This way, if x and y are 1-D arrays and ``z = correlate(x, y, 'full')`` 149 then 150 151 .. math:: 152 153 z[k] = (x * y)(k - N + 1) 154 = \sum_{l=0}^{||x||-1}x_l y_{l-k+N-1}^{*} 155 156 for :math:`k = 0, 1, ..., ||x|| + ||y|| - 2` 157 158 where :math:`||x||` is the length of ``x``, :math:`N = \max(||x||,||y||)`, 159 and :math:`y_m` is 0 when m is outside the range of y. 160 161 ``method='fft'`` only works for numerical arrays as it relies on 162 `fftconvolve`. In certain cases (i.e., arrays of objects or when 163 rounding integers can lose precision), ``method='direct'`` is always used. 164 165 When using "same" mode with even-length inputs, the outputs of `correlate` 166 and `correlate2d` differ: There is a 1-index offset between them. 167 168 Examples 169 -------- 170 Implement a matched filter using cross-correlation, to recover a signal 171 that has passed through a noisy channel. 172 173 >>> from scipy import signal 174 >>> import matplotlib.pyplot as plt 175 >>> rng = np.random.default_rng() 176 177 >>> sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 128) 178 >>> sig_noise = sig + rng.standard_normal(len(sig)) 179 >>> corr = signal.correlate(sig_noise, np.ones(128), mode='same') / 128 180 181 >>> clock = np.arange(64, len(sig), 128) 182 >>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True) 183 >>> ax_orig.plot(sig) 184 >>> ax_orig.plot(clock, sig[clock], 'ro') 185 >>> ax_orig.set_title('Original signal') 186 >>> ax_noise.plot(sig_noise) 187 >>> ax_noise.set_title('Signal with noise') 188 >>> ax_corr.plot(corr) 189 >>> ax_corr.plot(clock, corr[clock], 'ro') 190 >>> ax_corr.axhline(0.5, ls=':') 191 >>> ax_corr.set_title('Cross-correlated with rectangular pulse') 192 >>> ax_orig.margins(0, 0.1) 193 >>> fig.tight_layout() 194 >>> plt.show() 195 196 Compute the cross-correlation of a noisy signal with the original signal. 197 198 >>> x = np.arange(128) / 128 199 >>> sig = np.sin(2 * np.pi * x) 200 >>> sig_noise = sig + rng.standard_normal(len(sig)) 201 >>> corr = signal.correlate(sig_noise, sig) 202 >>> lags = signal.correlation_lags(len(sig), len(sig_noise)) 203 >>> corr /= np.max(corr) 204 205 >>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, figsize=(4.8, 4.8)) 206 >>> ax_orig.plot(sig) 207 >>> ax_orig.set_title('Original signal') 208 >>> ax_orig.set_xlabel('Sample Number') 209 >>> ax_noise.plot(sig_noise) 210 >>> ax_noise.set_title('Signal with noise') 211 >>> ax_noise.set_xlabel('Sample Number') 212 >>> ax_corr.plot(lags, corr) 213 >>> ax_corr.set_title('Cross-correlated signal') 214 >>> ax_corr.set_xlabel('Lag') 215 >>> ax_orig.margins(0, 0.1) 216 >>> ax_noise.margins(0, 0.1) 217 >>> ax_corr.margins(0, 0.1) 218 >>> fig.tight_layout() 219 >>> plt.show() 220 221 """ 222 in1 = np.asarray(in1) 223 in2 = np.asarray(in2) 224 225 if in1.ndim == in2.ndim == 0: 226 return in1 * in2.conj() 227 elif in1.ndim != in2.ndim: 228 raise ValueError("in1 and in2 should have the same dimensionality") 229 230 # Don't use _valfrommode, since correlate should not accept numeric modes 231 try: 232 val = _modedict[mode] 233 except KeyError as e: 234 raise ValueError("Acceptable mode flags are 'valid'," 235 " 'same', or 'full'.") from e 236 237 # this either calls fftconvolve or this function with method=='direct' 238 if method in ('fft', 'auto'): 239 return convolve(in1, _reverse_and_conj(in2), mode, method) 240 241 elif method == 'direct': 242 # fastpath to faster numpy.correlate for 1d inputs when possible 243 if _np_conv_ok(in1, in2, mode): 244 return np.correlate(in1, in2, mode) 245 246 # _correlateND is far slower when in2.size > in1.size, so swap them 247 # and then undo the effect afterward if mode == 'full'. Also, it fails 248 # with 'valid' mode if in2 is larger than in1, so swap those, too. 249 # Don't swap inputs for 'same' mode, since shape of in1 matters. 250 swapped_inputs = ((mode == 'full') and (in2.size > in1.size) or 251 _inputs_swap_needed(mode, in1.shape, in2.shape)) 252 253 if swapped_inputs: 254 in1, in2 = in2, in1 255 256 if mode == 'valid': 257 ps = [i - j + 1 for i, j in zip(in1.shape, in2.shape)] 258 out = np.empty(ps, in1.dtype) 259 260 z = sigtools._correlateND(in1, in2, out, val) 261 262 else: 263 ps = [i + j - 1 for i, j in zip(in1.shape, in2.shape)] 264 265 # zero pad input 266 in1zpadded = np.zeros(ps, in1.dtype) 267 sc = tuple(slice(0, i) for i in in1.shape) 268 in1zpadded[sc] = in1.copy() 269 270 if mode == 'full': 271 out = np.empty(ps, in1.dtype) 272 elif mode == 'same': 273 out = np.empty(in1.shape, in1.dtype) 274 275 z = sigtools._correlateND(in1zpadded, in2, out, val) 276 277 if swapped_inputs: 278 # Reverse and conjugate to undo the effect of swapping inputs 279 z = _reverse_and_conj(z) 280 281 return z 282 283 else: 284 raise ValueError("Acceptable method flags are 'auto'," 285 " 'direct', or 'fft'.") 286 287 288def correlation_lags(in1_len, in2_len, mode='full'): 289 r""" 290 Calculates the lag / displacement indices array for 1D cross-correlation. 291 292 Parameters 293 ---------- 294 in1_size : int 295 First input size. 296 in2_size : int 297 Second input size. 298 mode : str {'full', 'valid', 'same'}, optional 299 A string indicating the size of the output. 300 See the documentation `correlate` for more information. 301 302 See Also 303 -------- 304 correlate : Compute the N-dimensional cross-correlation. 305 306 Returns 307 ------- 308 lags : array 309 Returns an array containing cross-correlation lag/displacement indices. 310 Indices can be indexed with the np.argmax of the correlation to return 311 the lag/displacement. 312 313 Notes 314 ----- 315 Cross-correlation for continuous functions :math:`f` and :math:`g` is 316 defined as: 317 318 .. math :: 319 320 \left ( f\star g \right )\left ( \tau \right ) 321 \triangleq \int_{t_0}^{t_0 +T} 322 \overline{f\left ( t \right )}g\left ( t+\tau \right )dt 323 324 Where :math:`\tau` is defined as the displacement, also known as the lag. 325 326 Cross correlation for discrete functions :math:`f` and :math:`g` is 327 defined as: 328 329 .. math :: 330 \left ( f\star g \right )\left [ n \right ] 331 \triangleq \sum_{-\infty}^{\infty} 332 \overline{f\left [ m \right ]}g\left [ m+n \right ] 333 334 Where :math:`n` is the lag. 335 336 Examples 337 -------- 338 Cross-correlation of a signal with its time-delayed self. 339 340 >>> from scipy import signal 341 >>> from numpy.random import default_rng 342 >>> rng = default_rng() 343 >>> x = rng.standard_normal(1000) 344 >>> y = np.concatenate([rng.standard_normal(100), x]) 345 >>> correlation = signal.correlate(x, y, mode="full") 346 >>> lags = signal.correlation_lags(x.size, y.size, mode="full") 347 >>> lag = lags[np.argmax(correlation)] 348 """ 349 350 # calculate lag ranges in different modes of operation 351 if mode == "full": 352 # the output is the full discrete linear convolution 353 # of the inputs. (Default) 354 lags = np.arange(-in2_len + 1, in1_len) 355 elif mode == "same": 356 # the output is the same size as `in1`, centered 357 # with respect to the 'full' output. 358 # calculate the full output 359 lags = np.arange(-in2_len + 1, in1_len) 360 # determine the midpoint in the full output 361 mid = lags.size // 2 362 # determine lag_bound to be used with respect 363 # to the midpoint 364 lag_bound = in1_len // 2 365 # calculate lag ranges for even and odd scenarios 366 if in1_len % 2 == 0: 367 lags = lags[(mid-lag_bound):(mid+lag_bound)] 368 else: 369 lags = lags[(mid-lag_bound):(mid+lag_bound)+1] 370 elif mode == "valid": 371 # the output consists only of those elements that do not 372 # rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 373 # must be at least as large as the other in every dimension. 374 375 # the lag_bound will be either negative or positive 376 # this let's us infer how to present the lag range 377 lag_bound = in1_len - in2_len 378 if lag_bound >= 0: 379 lags = np.arange(lag_bound + 1) 380 else: 381 lags = np.arange(lag_bound, 1) 382 return lags 383 384 385def _centered(arr, newshape): 386 # Return the center newshape portion of the array. 387 newshape = np.asarray(newshape) 388 currshape = np.array(arr.shape) 389 startind = (currshape - newshape) // 2 390 endind = startind + newshape 391 myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] 392 return arr[tuple(myslice)] 393 394 395def _init_freq_conv_axes(in1, in2, mode, axes, sorted_axes=False): 396 """Handle the axes argument for frequency-domain convolution. 397 398 Returns the inputs and axes in a standard form, eliminating redundant axes, 399 swapping the inputs if necessary, and checking for various potential 400 errors. 401 402 Parameters 403 ---------- 404 in1 : array 405 First input. 406 in2 : array 407 Second input. 408 mode : str {'full', 'valid', 'same'}, optional 409 A string indicating the size of the output. 410 See the documentation `fftconvolve` for more information. 411 axes : list of ints 412 Axes over which to compute the FFTs. 413 sorted_axes : bool, optional 414 If `True`, sort the axes. 415 Default is `False`, do not sort. 416 417 Returns 418 ------- 419 in1 : array 420 The first input, possible swapped with the second input. 421 in2 : array 422 The second input, possible swapped with the first input. 423 axes : list of ints 424 Axes over which to compute the FFTs. 425 426 """ 427 s1 = in1.shape 428 s2 = in2.shape 429 noaxes = axes is None 430 431 _, axes = _init_nd_shape_and_axes(in1, shape=None, axes=axes) 432 433 if not noaxes and not len(axes): 434 raise ValueError("when provided, axes cannot be empty") 435 436 # Axes of length 1 can rely on broadcasting rules for multipy, 437 # no fft needed. 438 axes = [a for a in axes if s1[a] != 1 and s2[a] != 1] 439 440 if sorted_axes: 441 axes.sort() 442 443 if not all(s1[a] == s2[a] or s1[a] == 1 or s2[a] == 1 444 for a in range(in1.ndim) if a not in axes): 445 raise ValueError("incompatible shapes for in1 and in2:" 446 " {0} and {1}".format(s1, s2)) 447 448 # Check that input sizes are compatible with 'valid' mode. 449 if _inputs_swap_needed(mode, s1, s2, axes=axes): 450 # Convolution is commutative; order doesn't have any effect on output. 451 in1, in2 = in2, in1 452 453 return in1, in2, axes 454 455 456def _freq_domain_conv(in1, in2, axes, shape, calc_fast_len=False): 457 """Convolve two arrays in the frequency domain. 458 459 This function implements only base the FFT-related operations. 460 Specifically, it converts the signals to the frequency domain, multiplies 461 them, then converts them back to the time domain. Calculations of axes, 462 shapes, convolution mode, etc. are implemented in higher level-functions, 463 such as `fftconvolve` and `oaconvolve`. Those functions should be used 464 instead of this one. 465 466 Parameters 467 ---------- 468 in1 : array_like 469 First input. 470 in2 : array_like 471 Second input. Should have the same number of dimensions as `in1`. 472 axes : array_like of ints 473 Axes over which to compute the FFTs. 474 shape : array_like of ints 475 The sizes of the FFTs. 476 calc_fast_len : bool, optional 477 If `True`, set each value of `shape` to the next fast FFT length. 478 Default is `False`, use `axes` as-is. 479 480 Returns 481 ------- 482 out : array 483 An N-dimensional array containing the discrete linear convolution of 484 `in1` with `in2`. 485 486 """ 487 if not len(axes): 488 return in1 * in2 489 490 complex_result = (in1.dtype.kind == 'c' or in2.dtype.kind == 'c') 491 492 if calc_fast_len: 493 # Speed up FFT by padding to optimal size. 494 fshape = [ 495 sp_fft.next_fast_len(shape[a], not complex_result) for a in axes] 496 else: 497 fshape = shape 498 499 if not complex_result: 500 fft, ifft = sp_fft.rfftn, sp_fft.irfftn 501 else: 502 fft, ifft = sp_fft.fftn, sp_fft.ifftn 503 504 sp1 = fft(in1, fshape, axes=axes) 505 sp2 = fft(in2, fshape, axes=axes) 506 507 ret = ifft(sp1 * sp2, fshape, axes=axes) 508 509 if calc_fast_len: 510 fslice = tuple([slice(sz) for sz in shape]) 511 ret = ret[fslice] 512 513 return ret 514 515 516def _apply_conv_mode(ret, s1, s2, mode, axes): 517 """Calculate the convolution result shape based on the `mode` argument. 518 519 Returns the result sliced to the correct size for the given mode. 520 521 Parameters 522 ---------- 523 ret : array 524 The result array, with the appropriate shape for the 'full' mode. 525 s1 : list of int 526 The shape of the first input. 527 s2 : list of int 528 The shape of the second input. 529 mode : str {'full', 'valid', 'same'} 530 A string indicating the size of the output. 531 See the documentation `fftconvolve` for more information. 532 axes : list of ints 533 Axes over which to compute the convolution. 534 535 Returns 536 ------- 537 ret : array 538 A copy of `res`, sliced to the correct size for the given `mode`. 539 540 """ 541 if mode == "full": 542 return ret.copy() 543 elif mode == "same": 544 return _centered(ret, s1).copy() 545 elif mode == "valid": 546 shape_valid = [ret.shape[a] if a not in axes else s1[a] - s2[a] + 1 547 for a in range(ret.ndim)] 548 return _centered(ret, shape_valid).copy() 549 else: 550 raise ValueError("acceptable mode flags are 'valid'," 551 " 'same', or 'full'") 552 553 554def fftconvolve(in1, in2, mode="full", axes=None): 555 """Convolve two N-dimensional arrays using FFT. 556 557 Convolve `in1` and `in2` using the fast Fourier transform method, with 558 the output size determined by the `mode` argument. 559 560 This is generally much faster than `convolve` for large arrays (n > ~500), 561 but can be slower when only a few output values are needed, and can only 562 output float arrays (int or object array inputs will be cast to float). 563 564 As of v0.19, `convolve` automatically chooses this method or the direct 565 method based on an estimation of which is faster. 566 567 Parameters 568 ---------- 569 in1 : array_like 570 First input. 571 in2 : array_like 572 Second input. Should have the same number of dimensions as `in1`. 573 mode : str {'full', 'valid', 'same'}, optional 574 A string indicating the size of the output: 575 576 ``full`` 577 The output is the full discrete linear convolution 578 of the inputs. (Default) 579 ``valid`` 580 The output consists only of those elements that do not 581 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 582 must be at least as large as the other in every dimension. 583 ``same`` 584 The output is the same size as `in1`, centered 585 with respect to the 'full' output. 586 axes : int or array_like of ints or None, optional 587 Axes over which to compute the convolution. 588 The default is over all axes. 589 590 Returns 591 ------- 592 out : array 593 An N-dimensional array containing a subset of the discrete linear 594 convolution of `in1` with `in2`. 595 596 See Also 597 -------- 598 convolve : Uses the direct convolution or FFT convolution algorithm 599 depending on which is faster. 600 oaconvolve : Uses the overlap-add method to do convolution, which is 601 generally faster when the input arrays are large and 602 significantly different in size. 603 604 Examples 605 -------- 606 Autocorrelation of white noise is an impulse. 607 608 >>> from scipy import signal 609 >>> rng = np.random.default_rng() 610 >>> sig = rng.standard_normal(1000) 611 >>> autocorr = signal.fftconvolve(sig, sig[::-1], mode='full') 612 613 >>> import matplotlib.pyplot as plt 614 >>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1) 615 >>> ax_orig.plot(sig) 616 >>> ax_orig.set_title('White noise') 617 >>> ax_mag.plot(np.arange(-len(sig)+1,len(sig)), autocorr) 618 >>> ax_mag.set_title('Autocorrelation') 619 >>> fig.tight_layout() 620 >>> fig.show() 621 622 Gaussian blur implemented using FFT convolution. Notice the dark borders 623 around the image, due to the zero-padding beyond its boundaries. 624 The `convolve2d` function allows for other types of image boundaries, 625 but is far slower. 626 627 >>> from scipy import misc 628 >>> face = misc.face(gray=True) 629 >>> kernel = np.outer(signal.windows.gaussian(70, 8), 630 ... signal.windows.gaussian(70, 8)) 631 >>> blurred = signal.fftconvolve(face, kernel, mode='same') 632 633 >>> fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1, 634 ... figsize=(6, 15)) 635 >>> ax_orig.imshow(face, cmap='gray') 636 >>> ax_orig.set_title('Original') 637 >>> ax_orig.set_axis_off() 638 >>> ax_kernel.imshow(kernel, cmap='gray') 639 >>> ax_kernel.set_title('Gaussian kernel') 640 >>> ax_kernel.set_axis_off() 641 >>> ax_blurred.imshow(blurred, cmap='gray') 642 >>> ax_blurred.set_title('Blurred') 643 >>> ax_blurred.set_axis_off() 644 >>> fig.show() 645 646 """ 647 in1 = np.asarray(in1) 648 in2 = np.asarray(in2) 649 650 if in1.ndim == in2.ndim == 0: # scalar inputs 651 return in1 * in2 652 elif in1.ndim != in2.ndim: 653 raise ValueError("in1 and in2 should have the same dimensionality") 654 elif in1.size == 0 or in2.size == 0: # empty arrays 655 return np.array([]) 656 657 in1, in2, axes = _init_freq_conv_axes(in1, in2, mode, axes, 658 sorted_axes=False) 659 660 s1 = in1.shape 661 s2 = in2.shape 662 663 shape = [max((s1[i], s2[i])) if i not in axes else s1[i] + s2[i] - 1 664 for i in range(in1.ndim)] 665 666 ret = _freq_domain_conv(in1, in2, axes, shape, calc_fast_len=True) 667 668 return _apply_conv_mode(ret, s1, s2, mode, axes) 669 670 671def _calc_oa_lens(s1, s2): 672 """Calculate the optimal FFT lengths for overlapp-add convolution. 673 674 The calculation is done for a single dimension. 675 676 Parameters 677 ---------- 678 s1 : int 679 Size of the dimension for the first array. 680 s2 : int 681 Size of the dimension for the second array. 682 683 Returns 684 ------- 685 block_size : int 686 The size of the FFT blocks. 687 overlap : int 688 The amount of overlap between two blocks. 689 in1_step : int 690 The size of each step for the first array. 691 in2_step : int 692 The size of each step for the first array. 693 694 """ 695 # Set up the arguments for the conventional FFT approach. 696 fallback = (s1+s2-1, None, s1, s2) 697 698 # Use conventional FFT convolve if sizes are same. 699 if s1 == s2 or s1 == 1 or s2 == 1: 700 return fallback 701 702 if s2 > s1: 703 s1, s2 = s2, s1 704 swapped = True 705 else: 706 swapped = False 707 708 # There cannot be a useful block size if s2 is more than half of s1. 709 if s2 >= s1/2: 710 return fallback 711 712 # Derivation of optimal block length 713 # For original formula see: 714 # https://en.wikipedia.org/wiki/Overlap-add_method 715 # 716 # Formula: 717 # K = overlap = s2-1 718 # N = block_size 719 # C = complexity 720 # e = exponential, exp(1) 721 # 722 # C = (N*(log2(N)+1))/(N-K) 723 # C = (N*log2(2N))/(N-K) 724 # C = N/(N-K) * log2(2N) 725 # C1 = N/(N-K) 726 # C2 = log2(2N) = ln(2N)/ln(2) 727 # 728 # dC1/dN = (1*(N-K)-N)/(N-K)^2 = -K/(N-K)^2 729 # dC2/dN = 2/(2*N*ln(2)) = 1/(N*ln(2)) 730 # 731 # dC/dN = dC1/dN*C2 + dC2/dN*C1 732 # dC/dN = -K*ln(2N)/(ln(2)*(N-K)^2) + N/(N*ln(2)*(N-K)) 733 # dC/dN = -K*ln(2N)/(ln(2)*(N-K)^2) + 1/(ln(2)*(N-K)) 734 # dC/dN = -K*ln(2N)/(ln(2)*(N-K)^2) + (N-K)/(ln(2)*(N-K)^2) 735 # dC/dN = (-K*ln(2N) + (N-K)/(ln(2)*(N-K)^2) 736 # dC/dN = (N - K*ln(2N) - K)/(ln(2)*(N-K)^2) 737 # 738 # Solve for minimum, where dC/dN = 0 739 # 0 = (N - K*ln(2N) - K)/(ln(2)*(N-K)^2) 740 # 0 * ln(2)*(N-K)^2 = N - K*ln(2N) - K 741 # 0 = N - K*ln(2N) - K 742 # 0 = N - K*(ln(2N) + 1) 743 # 0 = N - K*ln(2Ne) 744 # N = K*ln(2Ne) 745 # N/K = ln(2Ne) 746 # 747 # e^(N/K) = e^ln(2Ne) 748 # e^(N/K) = 2Ne 749 # 1/e^(N/K) = 1/(2*N*e) 750 # e^(N/-K) = 1/(2*N*e) 751 # e^(N/-K) = K/N*1/(2*K*e) 752 # N/K*e^(N/-K) = 1/(2*e*K) 753 # N/-K*e^(N/-K) = -1/(2*e*K) 754 # 755 # Using Lambert W function 756 # https://en.wikipedia.org/wiki/Lambert_W_function 757 # x = W(y) It is the solution to y = x*e^x 758 # x = N/-K 759 # y = -1/(2*e*K) 760 # 761 # N/-K = W(-1/(2*e*K)) 762 # 763 # N = -K*W(-1/(2*e*K)) 764 overlap = s2-1 765 opt_size = -overlap*lambertw(-1/(2*math.e*overlap), k=-1).real 766 block_size = sp_fft.next_fast_len(math.ceil(opt_size)) 767 768 # Use conventional FFT convolve if there is only going to be one block. 769 if block_size >= s1: 770 return fallback 771 772 if not swapped: 773 in1_step = block_size-s2+1 774 in2_step = s2 775 else: 776 in1_step = s2 777 in2_step = block_size-s2+1 778 779 return block_size, overlap, in1_step, in2_step 780 781 782def oaconvolve(in1, in2, mode="full", axes=None): 783 """Convolve two N-dimensional arrays using the overlap-add method. 784 785 Convolve `in1` and `in2` using the overlap-add method, with 786 the output size determined by the `mode` argument. 787 788 This is generally much faster than `convolve` for large arrays (n > ~500), 789 and generally much faster than `fftconvolve` when one array is much 790 larger than the other, but can be slower when only a few output values are 791 needed or when the arrays are very similar in shape, and can only 792 output float arrays (int or object array inputs will be cast to float). 793 794 Parameters 795 ---------- 796 in1 : array_like 797 First input. 798 in2 : array_like 799 Second input. Should have the same number of dimensions as `in1`. 800 mode : str {'full', 'valid', 'same'}, optional 801 A string indicating the size of the output: 802 803 ``full`` 804 The output is the full discrete linear convolution 805 of the inputs. (Default) 806 ``valid`` 807 The output consists only of those elements that do not 808 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 809 must be at least as large as the other in every dimension. 810 ``same`` 811 The output is the same size as `in1`, centered 812 with respect to the 'full' output. 813 axes : int or array_like of ints or None, optional 814 Axes over which to compute the convolution. 815 The default is over all axes. 816 817 Returns 818 ------- 819 out : array 820 An N-dimensional array containing a subset of the discrete linear 821 convolution of `in1` with `in2`. 822 823 See Also 824 -------- 825 convolve : Uses the direct convolution or FFT convolution algorithm 826 depending on which is faster. 827 fftconvolve : An implementation of convolution using FFT. 828 829 Notes 830 ----- 831 .. versionadded:: 1.4.0 832 833 Examples 834 -------- 835 Convolve a 100,000 sample signal with a 512-sample filter. 836 837 >>> from scipy import signal 838 >>> rng = np.random.default_rng() 839 >>> sig = rng.standard_normal(100000) 840 >>> filt = signal.firwin(512, 0.01) 841 >>> fsig = signal.oaconvolve(sig, filt) 842 843 >>> import matplotlib.pyplot as plt 844 >>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1) 845 >>> ax_orig.plot(sig) 846 >>> ax_orig.set_title('White noise') 847 >>> ax_mag.plot(fsig) 848 >>> ax_mag.set_title('Filtered noise') 849 >>> fig.tight_layout() 850 >>> fig.show() 851 852 References 853 ---------- 854 .. [1] Wikipedia, "Overlap-add_method". 855 https://en.wikipedia.org/wiki/Overlap-add_method 856 .. [2] Richard G. Lyons. Understanding Digital Signal Processing, 857 Third Edition, 2011. Chapter 13.10. 858 ISBN 13: 978-0137-02741-5 859 860 """ 861 in1 = np.asarray(in1) 862 in2 = np.asarray(in2) 863 864 if in1.ndim == in2.ndim == 0: # scalar inputs 865 return in1 * in2 866 elif in1.ndim != in2.ndim: 867 raise ValueError("in1 and in2 should have the same dimensionality") 868 elif in1.size == 0 or in2.size == 0: # empty arrays 869 return np.array([]) 870 elif in1.shape == in2.shape: # Equivalent to fftconvolve 871 return fftconvolve(in1, in2, mode=mode, axes=axes) 872 873 in1, in2, axes = _init_freq_conv_axes(in1, in2, mode, axes, 874 sorted_axes=True) 875 876 s1 = in1.shape 877 s2 = in2.shape 878 879 if not axes: 880 ret = in1 * in2 881 return _apply_conv_mode(ret, s1, s2, mode, axes) 882 883 # Calculate this now since in1 is changed later 884 shape_final = [None if i not in axes else 885 s1[i] + s2[i] - 1 for i in range(in1.ndim)] 886 887 # Calculate the block sizes for the output, steps, first and second inputs. 888 # It is simpler to calculate them all together than doing them in separate 889 # loops due to all the special cases that need to be handled. 890 optimal_sizes = ((-1, -1, s1[i], s2[i]) if i not in axes else 891 _calc_oa_lens(s1[i], s2[i]) for i in range(in1.ndim)) 892 block_size, overlaps, \ 893 in1_step, in2_step = zip(*optimal_sizes) 894 895 # Fall back to fftconvolve if there is only one block in every dimension. 896 if in1_step == s1 and in2_step == s2: 897 return fftconvolve(in1, in2, mode=mode, axes=axes) 898 899 # Figure out the number of steps and padding. 900 # This would get too complicated in a list comprehension. 901 nsteps1 = [] 902 nsteps2 = [] 903 pad_size1 = [] 904 pad_size2 = [] 905 for i in range(in1.ndim): 906 if i not in axes: 907 pad_size1 += [(0, 0)] 908 pad_size2 += [(0, 0)] 909 continue 910 911 if s1[i] > in1_step[i]: 912 curnstep1 = math.ceil((s1[i]+1)/in1_step[i]) 913 if (block_size[i] - overlaps[i])*curnstep1 < shape_final[i]: 914 curnstep1 += 1 915 916 curpad1 = curnstep1*in1_step[i] - s1[i] 917 else: 918 curnstep1 = 1 919 curpad1 = 0 920 921 if s2[i] > in2_step[i]: 922 curnstep2 = math.ceil((s2[i]+1)/in2_step[i]) 923 if (block_size[i] - overlaps[i])*curnstep2 < shape_final[i]: 924 curnstep2 += 1 925 926 curpad2 = curnstep2*in2_step[i] - s2[i] 927 else: 928 curnstep2 = 1 929 curpad2 = 0 930 931 nsteps1 += [curnstep1] 932 nsteps2 += [curnstep2] 933 pad_size1 += [(0, curpad1)] 934 pad_size2 += [(0, curpad2)] 935 936 # Pad the array to a size that can be reshaped to the desired shape 937 # if necessary. 938 if not all(curpad == (0, 0) for curpad in pad_size1): 939 in1 = np.pad(in1, pad_size1, mode='constant', constant_values=0) 940 941 if not all(curpad == (0, 0) for curpad in pad_size2): 942 in2 = np.pad(in2, pad_size2, mode='constant', constant_values=0) 943 944 # Reshape the overlap-add parts to input block sizes. 945 split_axes = [iax+i for i, iax in enumerate(axes)] 946 fft_axes = [iax+1 for iax in split_axes] 947 948 # We need to put each new dimension before the corresponding dimension 949 # being reshaped in order to get the data in the right layout at the end. 950 reshape_size1 = list(in1_step) 951 reshape_size2 = list(in2_step) 952 for i, iax in enumerate(split_axes): 953 reshape_size1.insert(iax, nsteps1[i]) 954 reshape_size2.insert(iax, nsteps2[i]) 955 956 in1 = in1.reshape(*reshape_size1) 957 in2 = in2.reshape(*reshape_size2) 958 959 # Do the convolution. 960 fft_shape = [block_size[i] for i in axes] 961 ret = _freq_domain_conv(in1, in2, fft_axes, fft_shape, calc_fast_len=False) 962 963 # Do the overlap-add. 964 for ax, ax_fft, ax_split in zip(axes, fft_axes, split_axes): 965 overlap = overlaps[ax] 966 if overlap is None: 967 continue 968 969 ret, overpart = np.split(ret, [-overlap], ax_fft) 970 overpart = np.split(overpart, [-1], ax_split)[0] 971 972 ret_overpart = np.split(ret, [overlap], ax_fft)[0] 973 ret_overpart = np.split(ret_overpart, [1], ax_split)[1] 974 ret_overpart += overpart 975 976 # Reshape back to the correct dimensionality. 977 shape_ret = [ret.shape[i] if i not in fft_axes else 978 ret.shape[i]*ret.shape[i-1] 979 for i in range(ret.ndim) if i not in split_axes] 980 ret = ret.reshape(*shape_ret) 981 982 # Slice to the correct size. 983 slice_final = tuple([slice(islice) for islice in shape_final]) 984 ret = ret[slice_final] 985 986 return _apply_conv_mode(ret, s1, s2, mode, axes) 987 988 989def _numeric_arrays(arrays, kinds='buifc'): 990 """ 991 See if a list of arrays are all numeric. 992 993 Parameters 994 ---------- 995 ndarrays : array or list of arrays 996 arrays to check if numeric. 997 numeric_kinds : string-like 998 The dtypes of the arrays to be checked. If the dtype.kind of 999 the ndarrays are not in this string the function returns False and 1000 otherwise returns True. 1001 """ 1002 if type(arrays) == np.ndarray: 1003 return arrays.dtype.kind in kinds 1004 for array_ in arrays: 1005 if array_.dtype.kind not in kinds: 1006 return False 1007 return True 1008 1009 1010def _conv_ops(x_shape, h_shape, mode): 1011 """ 1012 Find the number of operations required for direct/fft methods of 1013 convolution. The direct operations were recorded by making a dummy class to 1014 record the number of operations by overriding ``__mul__`` and ``__add__``. 1015 The FFT operations rely on the (well-known) computational complexity of the 1016 FFT (and the implementation of ``_freq_domain_conv``). 1017 1018 """ 1019 if mode == "full": 1020 out_shape = [n + k - 1 for n, k in zip(x_shape, h_shape)] 1021 elif mode == "valid": 1022 out_shape = [abs(n - k) + 1 for n, k in zip(x_shape, h_shape)] 1023 elif mode == "same": 1024 out_shape = x_shape 1025 else: 1026 raise ValueError("Acceptable mode flags are 'valid'," 1027 " 'same', or 'full', not mode={}".format(mode)) 1028 1029 s1, s2 = x_shape, h_shape 1030 if len(x_shape) == 1: 1031 s1, s2 = s1[0], s2[0] 1032 if mode == "full": 1033 direct_ops = s1 * s2 1034 elif mode == "valid": 1035 direct_ops = (s2 - s1 + 1) * s1 if s2 >= s1 else (s1 - s2 + 1) * s2 1036 elif mode == "same": 1037 direct_ops = (s1 * s2 if s1 < s2 else 1038 s1 * s2 - (s2 // 2) * ((s2 + 1) // 2)) 1039 else: 1040 if mode == "full": 1041 direct_ops = min(_prod(s1), _prod(s2)) * _prod(out_shape) 1042 elif mode == "valid": 1043 direct_ops = min(_prod(s1), _prod(s2)) * _prod(out_shape) 1044 elif mode == "same": 1045 direct_ops = _prod(s1) * _prod(s2) 1046 1047 full_out_shape = [n + k - 1 for n, k in zip(x_shape, h_shape)] 1048 N = _prod(full_out_shape) 1049 fft_ops = 3 * N * np.log(N) # 3 separate FFTs of size full_out_shape 1050 return fft_ops, direct_ops 1051 1052 1053def _fftconv_faster(x, h, mode): 1054 """ 1055 See if using fftconvolve or convolve is faster. 1056 1057 Parameters 1058 ---------- 1059 x : np.ndarray 1060 Signal 1061 h : np.ndarray 1062 Kernel 1063 mode : str 1064 Mode passed to convolve 1065 1066 Returns 1067 ------- 1068 fft_faster : bool 1069 1070 Notes 1071 ----- 1072 See docstring of `choose_conv_method` for details on tuning hardware. 1073 1074 See pull request 11031 for more detail: 1075 https://github.com/scipy/scipy/pull/11031. 1076 1077 """ 1078 fft_ops, direct_ops = _conv_ops(x.shape, h.shape, mode) 1079 offset = -1e-3 if x.ndim == 1 else -1e-4 1080 constants = { 1081 "valid": (1.89095737e-9, 2.1364985e-10, offset), 1082 "full": (1.7649070e-9, 2.1414831e-10, offset), 1083 "same": (3.2646654e-9, 2.8478277e-10, offset) 1084 if h.size <= x.size 1085 else (3.21635404e-9, 1.1773253e-8, -1e-5), 1086 } if x.ndim == 1 else { 1087 "valid": (1.85927e-9, 2.11242e-8, offset), 1088 "full": (1.99817e-9, 1.66174e-8, offset), 1089 "same": (2.04735e-9, 1.55367e-8, offset), 1090 } 1091 O_fft, O_direct, O_offset = constants[mode] 1092 return O_fft * fft_ops < O_direct * direct_ops + O_offset 1093 1094 1095def _reverse_and_conj(x): 1096 """ 1097 Reverse array `x` in all dimensions and perform the complex conjugate 1098 """ 1099 reverse = (slice(None, None, -1),) * x.ndim 1100 return x[reverse].conj() 1101 1102 1103def _np_conv_ok(volume, kernel, mode): 1104 """ 1105 See if numpy supports convolution of `volume` and `kernel` (i.e. both are 1106 1D ndarrays and of the appropriate shape). NumPy's 'same' mode uses the 1107 size of the larger input, while SciPy's uses the size of the first input. 1108 1109 Invalid mode strings will return False and be caught by the calling func. 1110 """ 1111 if volume.ndim == kernel.ndim == 1: 1112 if mode in ('full', 'valid'): 1113 return True 1114 elif mode == 'same': 1115 return volume.size >= kernel.size 1116 else: 1117 return False 1118 1119 1120def _timeit_fast(stmt="pass", setup="pass", repeat=3): 1121 """ 1122 Returns the time the statement/function took, in seconds. 1123 1124 Faster, less precise version of IPython's timeit. `stmt` can be a statement 1125 written as a string or a callable. 1126 1127 Will do only 1 loop (like IPython's timeit) with no repetitions 1128 (unlike IPython) for very slow functions. For fast functions, only does 1129 enough loops to take 5 ms, which seems to produce similar results (on 1130 Windows at least), and avoids doing an extraneous cycle that isn't 1131 measured. 1132 1133 """ 1134 timer = timeit.Timer(stmt, setup) 1135 1136 # determine number of calls per rep so total time for 1 rep >= 5 ms 1137 x = 0 1138 for p in range(0, 10): 1139 number = 10**p 1140 x = timer.timeit(number) # seconds 1141 if x >= 5e-3 / 10: # 5 ms for final test, 1/10th that for this one 1142 break 1143 if x > 1: # second 1144 # If it's macroscopic, don't bother with repetitions 1145 best = x 1146 else: 1147 number *= 10 1148 r = timer.repeat(repeat, number) 1149 best = min(r) 1150 1151 sec = best / number 1152 return sec 1153 1154 1155def choose_conv_method(in1, in2, mode='full', measure=False): 1156 """ 1157 Find the fastest convolution/correlation method. 1158 1159 This primarily exists to be called during the ``method='auto'`` option in 1160 `convolve` and `correlate`. It can also be used to determine the value of 1161 ``method`` for many different convolutions of the same dtype/shape. 1162 In addition, it supports timing the convolution to adapt the value of 1163 ``method`` to a particular set of inputs and/or hardware. 1164 1165 Parameters 1166 ---------- 1167 in1 : array_like 1168 The first argument passed into the convolution function. 1169 in2 : array_like 1170 The second argument passed into the convolution function. 1171 mode : str {'full', 'valid', 'same'}, optional 1172 A string indicating the size of the output: 1173 1174 ``full`` 1175 The output is the full discrete linear convolution 1176 of the inputs. (Default) 1177 ``valid`` 1178 The output consists only of those elements that do not 1179 rely on the zero-padding. 1180 ``same`` 1181 The output is the same size as `in1`, centered 1182 with respect to the 'full' output. 1183 measure : bool, optional 1184 If True, run and time the convolution of `in1` and `in2` with both 1185 methods and return the fastest. If False (default), predict the fastest 1186 method using precomputed values. 1187 1188 Returns 1189 ------- 1190 method : str 1191 A string indicating which convolution method is fastest, either 1192 'direct' or 'fft' 1193 times : dict, optional 1194 A dictionary containing the times (in seconds) needed for each method. 1195 This value is only returned if ``measure=True``. 1196 1197 See Also 1198 -------- 1199 convolve 1200 correlate 1201 1202 Notes 1203 ----- 1204 Generally, this method is 99% accurate for 2D signals and 85% accurate 1205 for 1D signals for randomly chosen input sizes. For precision, use 1206 ``measure=True`` to find the fastest method by timing the convolution. 1207 This can be used to avoid the minimal overhead of finding the fastest 1208 ``method`` later, or to adapt the value of ``method`` to a particular set 1209 of inputs. 1210 1211 Experiments were run on an Amazon EC2 r5a.2xlarge machine to test this 1212 function. These experiments measured the ratio between the time required 1213 when using ``method='auto'`` and the time required for the fastest method 1214 (i.e., ``ratio = time_auto / min(time_fft, time_direct)``). In these 1215 experiments, we found: 1216 1217 * There is a 95% chance of this ratio being less than 1.5 for 1D signals 1218 and a 99% chance of being less than 2.5 for 2D signals. 1219 * The ratio was always less than 2.5/5 for 1D/2D signals respectively. 1220 * This function is most inaccurate for 1D convolutions that take between 1 1221 and 10 milliseconds with ``method='direct'``. A good proxy for this 1222 (at least in our experiments) is ``1e6 <= in1.size * in2.size <= 1e7``. 1223 1224 The 2D results almost certainly generalize to 3D/4D/etc because the 1225 implementation is the same (the 1D implementation is different). 1226 1227 All the numbers above are specific to the EC2 machine. However, we did find 1228 that this function generalizes fairly decently across hardware. The speed 1229 tests were of similar quality (and even slightly better) than the same 1230 tests performed on the machine to tune this function's numbers (a mid-2014 1231 15-inch MacBook Pro with 16GB RAM and a 2.5GHz Intel i7 processor). 1232 1233 There are cases when `fftconvolve` supports the inputs but this function 1234 returns `direct` (e.g., to protect against floating point integer 1235 precision). 1236 1237 .. versionadded:: 0.19 1238 1239 Examples 1240 -------- 1241 Estimate the fastest method for a given input: 1242 1243 >>> from scipy import signal 1244 >>> rng = np.random.default_rng() 1245 >>> img = rng.random((32, 32)) 1246 >>> filter = rng.random((8, 8)) 1247 >>> method = signal.choose_conv_method(img, filter, mode='same') 1248 >>> method 1249 'fft' 1250 1251 This can then be applied to other arrays of the same dtype and shape: 1252 1253 >>> img2 = rng.random((32, 32)) 1254 >>> filter2 = rng.random((8, 8)) 1255 >>> corr2 = signal.correlate(img2, filter2, mode='same', method=method) 1256 >>> conv2 = signal.convolve(img2, filter2, mode='same', method=method) 1257 1258 The output of this function (``method``) works with `correlate` and 1259 `convolve`. 1260 1261 """ 1262 volume = np.asarray(in1) 1263 kernel = np.asarray(in2) 1264 1265 if measure: 1266 times = {} 1267 for method in ['fft', 'direct']: 1268 times[method] = _timeit_fast(lambda: convolve(volume, kernel, 1269 mode=mode, method=method)) 1270 1271 chosen_method = 'fft' if times['fft'] < times['direct'] else 'direct' 1272 return chosen_method, times 1273 1274 # for integer input, 1275 # catch when more precision required than float provides (representing an 1276 # integer as float can lose precision in fftconvolve if larger than 2**52) 1277 if any([_numeric_arrays([x], kinds='ui') for x in [volume, kernel]]): 1278 max_value = int(np.abs(volume).max()) * int(np.abs(kernel).max()) 1279 max_value *= int(min(volume.size, kernel.size)) 1280 if max_value > 2**np.finfo('float').nmant - 1: 1281 return 'direct' 1282 1283 if _numeric_arrays([volume, kernel], kinds='b'): 1284 return 'direct' 1285 1286 if _numeric_arrays([volume, kernel]): 1287 if _fftconv_faster(volume, kernel, mode): 1288 return 'fft' 1289 1290 return 'direct' 1291 1292 1293def convolve(in1, in2, mode='full', method='auto'): 1294 """ 1295 Convolve two N-dimensional arrays. 1296 1297 Convolve `in1` and `in2`, with the output size determined by the 1298 `mode` argument. 1299 1300 Parameters 1301 ---------- 1302 in1 : array_like 1303 First input. 1304 in2 : array_like 1305 Second input. Should have the same number of dimensions as `in1`. 1306 mode : str {'full', 'valid', 'same'}, optional 1307 A string indicating the size of the output: 1308 1309 ``full`` 1310 The output is the full discrete linear convolution 1311 of the inputs. (Default) 1312 ``valid`` 1313 The output consists only of those elements that do not 1314 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 1315 must be at least as large as the other in every dimension. 1316 ``same`` 1317 The output is the same size as `in1`, centered 1318 with respect to the 'full' output. 1319 method : str {'auto', 'direct', 'fft'}, optional 1320 A string indicating which method to use to calculate the convolution. 1321 1322 ``direct`` 1323 The convolution is determined directly from sums, the definition of 1324 convolution. 1325 ``fft`` 1326 The Fourier Transform is used to perform the convolution by calling 1327 `fftconvolve`. 1328 ``auto`` 1329 Automatically chooses direct or Fourier method based on an estimate 1330 of which is faster (default). See Notes for more detail. 1331 1332 .. versionadded:: 0.19.0 1333 1334 Returns 1335 ------- 1336 convolve : array 1337 An N-dimensional array containing a subset of the discrete linear 1338 convolution of `in1` with `in2`. 1339 1340 See Also 1341 -------- 1342 numpy.polymul : performs polynomial multiplication (same operation, but 1343 also accepts poly1d objects) 1344 choose_conv_method : chooses the fastest appropriate convolution method 1345 fftconvolve : Always uses the FFT method. 1346 oaconvolve : Uses the overlap-add method to do convolution, which is 1347 generally faster when the input arrays are large and 1348 significantly different in size. 1349 1350 Notes 1351 ----- 1352 By default, `convolve` and `correlate` use ``method='auto'``, which calls 1353 `choose_conv_method` to choose the fastest method using pre-computed 1354 values (`choose_conv_method` can also measure real-world timing with a 1355 keyword argument). Because `fftconvolve` relies on floating point numbers, 1356 there are certain constraints that may force `method=direct` (more detail 1357 in `choose_conv_method` docstring). 1358 1359 Examples 1360 -------- 1361 Smooth a square pulse using a Hann window: 1362 1363 >>> from scipy import signal 1364 >>> sig = np.repeat([0., 1., 0.], 100) 1365 >>> win = signal.windows.hann(50) 1366 >>> filtered = signal.convolve(sig, win, mode='same') / sum(win) 1367 1368 >>> import matplotlib.pyplot as plt 1369 >>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True) 1370 >>> ax_orig.plot(sig) 1371 >>> ax_orig.set_title('Original pulse') 1372 >>> ax_orig.margins(0, 0.1) 1373 >>> ax_win.plot(win) 1374 >>> ax_win.set_title('Filter impulse response') 1375 >>> ax_win.margins(0, 0.1) 1376 >>> ax_filt.plot(filtered) 1377 >>> ax_filt.set_title('Filtered signal') 1378 >>> ax_filt.margins(0, 0.1) 1379 >>> fig.tight_layout() 1380 >>> fig.show() 1381 1382 """ 1383 volume = np.asarray(in1) 1384 kernel = np.asarray(in2) 1385 1386 if volume.ndim == kernel.ndim == 0: 1387 return volume * kernel 1388 elif volume.ndim != kernel.ndim: 1389 raise ValueError("volume and kernel should have the same " 1390 "dimensionality") 1391 1392 if _inputs_swap_needed(mode, volume.shape, kernel.shape): 1393 # Convolution is commutative; order doesn't have any effect on output 1394 volume, kernel = kernel, volume 1395 1396 if method == 'auto': 1397 method = choose_conv_method(volume, kernel, mode=mode) 1398 1399 if method == 'fft': 1400 out = fftconvolve(volume, kernel, mode=mode) 1401 result_type = np.result_type(volume, kernel) 1402 if result_type.kind in {'u', 'i'}: 1403 out = np.around(out) 1404 return out.astype(result_type) 1405 elif method == 'direct': 1406 # fastpath to faster numpy.convolve for 1d inputs when possible 1407 if _np_conv_ok(volume, kernel, mode): 1408 return np.convolve(volume, kernel, mode) 1409 1410 return correlate(volume, _reverse_and_conj(kernel), mode, 'direct') 1411 else: 1412 raise ValueError("Acceptable method flags are 'auto'," 1413 " 'direct', or 'fft'.") 1414 1415 1416def order_filter(a, domain, rank): 1417 """ 1418 Perform an order filter on an N-D array. 1419 1420 Perform an order filter on the array in. The domain argument acts as a 1421 mask centered over each pixel. The non-zero elements of domain are 1422 used to select elements surrounding each input pixel which are placed 1423 in a list. The list is sorted, and the output for that pixel is the 1424 element corresponding to rank in the sorted list. 1425 1426 Parameters 1427 ---------- 1428 a : ndarray 1429 The N-dimensional input array. 1430 domain : array_like 1431 A mask array with the same number of dimensions as `a`. 1432 Each dimension should have an odd number of elements. 1433 rank : int 1434 A non-negative integer which selects the element from the 1435 sorted list (0 corresponds to the smallest element, 1 is the 1436 next smallest element, etc.). 1437 1438 Returns 1439 ------- 1440 out : ndarray 1441 The results of the order filter in an array with the same 1442 shape as `a`. 1443 1444 Examples 1445 -------- 1446 >>> from scipy import signal 1447 >>> x = np.arange(25).reshape(5, 5) 1448 >>> domain = np.identity(3) 1449 >>> x 1450 array([[ 0, 1, 2, 3, 4], 1451 [ 5, 6, 7, 8, 9], 1452 [10, 11, 12, 13, 14], 1453 [15, 16, 17, 18, 19], 1454 [20, 21, 22, 23, 24]]) 1455 >>> signal.order_filter(x, domain, 0) 1456 array([[ 0., 0., 0., 0., 0.], 1457 [ 0., 0., 1., 2., 0.], 1458 [ 0., 5., 6., 7., 0.], 1459 [ 0., 10., 11., 12., 0.], 1460 [ 0., 0., 0., 0., 0.]]) 1461 >>> signal.order_filter(x, domain, 2) 1462 array([[ 6., 7., 8., 9., 4.], 1463 [ 11., 12., 13., 14., 9.], 1464 [ 16., 17., 18., 19., 14.], 1465 [ 21., 22., 23., 24., 19.], 1466 [ 20., 21., 22., 23., 24.]]) 1467 1468 """ 1469 domain = np.asarray(domain) 1470 size = domain.shape 1471 for k in range(len(size)): 1472 if (size[k] % 2) != 1: 1473 raise ValueError("Each dimension of domain argument " 1474 " should have an odd number of elements.") 1475 return sigtools._order_filterND(a, domain, rank) 1476 1477 1478def medfilt(volume, kernel_size=None): 1479 """ 1480 Perform a median filter on an N-dimensional array. 1481 1482 Apply a median filter to the input array using a local window-size 1483 given by `kernel_size`. The array will automatically be zero-padded. 1484 1485 Parameters 1486 ---------- 1487 volume : array_like 1488 An N-dimensional input array. 1489 kernel_size : array_like, optional 1490 A scalar or an N-length list giving the size of the median filter 1491 window in each dimension. Elements of `kernel_size` should be odd. 1492 If `kernel_size` is a scalar, then this scalar is used as the size in 1493 each dimension. Default size is 3 for each dimension. 1494 1495 Returns 1496 ------- 1497 out : ndarray 1498 An array the same size as input containing the median filtered 1499 result. 1500 1501 Warns 1502 ----- 1503 UserWarning 1504 If array size is smaller than kernel size along any dimension 1505 1506 See Also 1507 -------- 1508 scipy.ndimage.median_filter 1509 scipy.signal.medfilt2d 1510 1511 Notes 1512 ----- 1513 The more general function `scipy.ndimage.median_filter` has a more 1514 efficient implementation of a median filter and therefore runs much faster. 1515 1516 For 2-dimensional images with ``uint8``, ``float32`` or ``float64`` dtypes, 1517 the specialised function `scipy.signal.medfilt2d` may be faster. 1518 1519 """ 1520 volume = np.atleast_1d(volume) 1521 if kernel_size is None: 1522 kernel_size = [3] * volume.ndim 1523 kernel_size = np.asarray(kernel_size) 1524 if kernel_size.shape == (): 1525 kernel_size = np.repeat(kernel_size.item(), volume.ndim) 1526 1527 for k in range(volume.ndim): 1528 if (kernel_size[k] % 2) != 1: 1529 raise ValueError("Each element of kernel_size should be odd.") 1530 if any(k > s for k, s in zip(kernel_size, volume.shape)): 1531 warnings.warn('kernel_size exceeds volume extent: the volume will be ' 1532 'zero-padded.') 1533 1534 domain = np.ones(kernel_size, dtype=volume.dtype) 1535 1536 numels = np.prod(kernel_size, axis=0) 1537 order = numels // 2 1538 return sigtools._order_filterND(volume, domain, order) 1539 1540 1541def wiener(im, mysize=None, noise=None): 1542 """ 1543 Perform a Wiener filter on an N-dimensional array. 1544 1545 Apply a Wiener filter to the N-dimensional array `im`. 1546 1547 Parameters 1548 ---------- 1549 im : ndarray 1550 An N-dimensional array. 1551 mysize : int or array_like, optional 1552 A scalar or an N-length list giving the size of the Wiener filter 1553 window in each dimension. Elements of mysize should be odd. 1554 If mysize is a scalar, then this scalar is used as the size 1555 in each dimension. 1556 noise : float, optional 1557 The noise-power to use. If None, then noise is estimated as the 1558 average of the local variance of the input. 1559 1560 Returns 1561 ------- 1562 out : ndarray 1563 Wiener filtered result with the same shape as `im`. 1564 1565 Examples 1566 -------- 1567 1568 >>> from scipy.misc import face 1569 >>> from scipy.signal.signaltools import wiener 1570 >>> import matplotlib.pyplot as plt 1571 >>> import numpy as np 1572 >>> rng = np.random.default_rng() 1573 >>> img = rng.random((40, 40)) #Create a random image 1574 >>> filtered_img = wiener(img, (5, 5)) #Filter the image 1575 >>> f, (plot1, plot2) = plt.subplots(1, 2) 1576 >>> plot1.imshow(img) 1577 >>> plot2.imshow(filtered_img) 1578 >>> plt.show() 1579 1580 Notes 1581 ----- 1582 This implementation is similar to wiener2 in Matlab/Octave. 1583 For more details see [1]_ 1584 1585 References 1586 ---------- 1587 .. [1] Lim, Jae S., Two-Dimensional Signal and Image Processing, 1588 Englewood Cliffs, NJ, Prentice Hall, 1990, p. 548. 1589 1590 1591 """ 1592 im = np.asarray(im) 1593 if mysize is None: 1594 mysize = [3] * im.ndim 1595 mysize = np.asarray(mysize) 1596 if mysize.shape == (): 1597 mysize = np.repeat(mysize.item(), im.ndim) 1598 1599 # Estimate the local mean 1600 lMean = correlate(im, np.ones(mysize), 'same') / np.prod(mysize, axis=0) 1601 1602 # Estimate the local variance 1603 lVar = (correlate(im ** 2, np.ones(mysize), 'same') / 1604 np.prod(mysize, axis=0) - lMean ** 2) 1605 1606 # Estimate the noise power if needed. 1607 if noise is None: 1608 noise = np.mean(np.ravel(lVar), axis=0) 1609 1610 res = (im - lMean) 1611 res *= (1 - noise / lVar) 1612 res += lMean 1613 out = np.where(lVar < noise, lMean, res) 1614 1615 return out 1616 1617 1618def convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0): 1619 """ 1620 Convolve two 2-dimensional arrays. 1621 1622 Convolve `in1` and `in2` with output size determined by `mode`, and 1623 boundary conditions determined by `boundary` and `fillvalue`. 1624 1625 Parameters 1626 ---------- 1627 in1 : array_like 1628 First input. 1629 in2 : array_like 1630 Second input. Should have the same number of dimensions as `in1`. 1631 mode : str {'full', 'valid', 'same'}, optional 1632 A string indicating the size of the output: 1633 1634 ``full`` 1635 The output is the full discrete linear convolution 1636 of the inputs. (Default) 1637 ``valid`` 1638 The output consists only of those elements that do not 1639 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 1640 must be at least as large as the other in every dimension. 1641 ``same`` 1642 The output is the same size as `in1`, centered 1643 with respect to the 'full' output. 1644 boundary : str {'fill', 'wrap', 'symm'}, optional 1645 A flag indicating how to handle boundaries: 1646 1647 ``fill`` 1648 pad input arrays with fillvalue. (default) 1649 ``wrap`` 1650 circular boundary conditions. 1651 ``symm`` 1652 symmetrical boundary conditions. 1653 1654 fillvalue : scalar, optional 1655 Value to fill pad input arrays with. Default is 0. 1656 1657 Returns 1658 ------- 1659 out : ndarray 1660 A 2-dimensional array containing a subset of the discrete linear 1661 convolution of `in1` with `in2`. 1662 1663 Examples 1664 -------- 1665 Compute the gradient of an image by 2D convolution with a complex Scharr 1666 operator. (Horizontal operator is real, vertical is imaginary.) Use 1667 symmetric boundary condition to avoid creating edges at the image 1668 boundaries. 1669 1670 >>> from scipy import signal 1671 >>> from scipy import misc 1672 >>> ascent = misc.ascent() 1673 >>> scharr = np.array([[ -3-3j, 0-10j, +3 -3j], 1674 ... [-10+0j, 0+ 0j, +10 +0j], 1675 ... [ -3+3j, 0+10j, +3 +3j]]) # Gx + j*Gy 1676 >>> grad = signal.convolve2d(ascent, scharr, boundary='symm', mode='same') 1677 1678 >>> import matplotlib.pyplot as plt 1679 >>> fig, (ax_orig, ax_mag, ax_ang) = plt.subplots(3, 1, figsize=(6, 15)) 1680 >>> ax_orig.imshow(ascent, cmap='gray') 1681 >>> ax_orig.set_title('Original') 1682 >>> ax_orig.set_axis_off() 1683 >>> ax_mag.imshow(np.absolute(grad), cmap='gray') 1684 >>> ax_mag.set_title('Gradient magnitude') 1685 >>> ax_mag.set_axis_off() 1686 >>> ax_ang.imshow(np.angle(grad), cmap='hsv') # hsv is cyclic, like angles 1687 >>> ax_ang.set_title('Gradient orientation') 1688 >>> ax_ang.set_axis_off() 1689 >>> fig.show() 1690 1691 """ 1692 in1 = np.asarray(in1) 1693 in2 = np.asarray(in2) 1694 1695 if not in1.ndim == in2.ndim == 2: 1696 raise ValueError('convolve2d inputs must both be 2-D arrays') 1697 1698 if _inputs_swap_needed(mode, in1.shape, in2.shape): 1699 in1, in2 = in2, in1 1700 1701 val = _valfrommode(mode) 1702 bval = _bvalfromboundary(boundary) 1703 out = sigtools._convolve2d(in1, in2, 1, val, bval, fillvalue) 1704 return out 1705 1706 1707def correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0): 1708 """ 1709 Cross-correlate two 2-dimensional arrays. 1710 1711 Cross correlate `in1` and `in2` with output size determined by `mode`, and 1712 boundary conditions determined by `boundary` and `fillvalue`. 1713 1714 Parameters 1715 ---------- 1716 in1 : array_like 1717 First input. 1718 in2 : array_like 1719 Second input. Should have the same number of dimensions as `in1`. 1720 mode : str {'full', 'valid', 'same'}, optional 1721 A string indicating the size of the output: 1722 1723 ``full`` 1724 The output is the full discrete linear cross-correlation 1725 of the inputs. (Default) 1726 ``valid`` 1727 The output consists only of those elements that do not 1728 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 1729 must be at least as large as the other in every dimension. 1730 ``same`` 1731 The output is the same size as `in1`, centered 1732 with respect to the 'full' output. 1733 boundary : str {'fill', 'wrap', 'symm'}, optional 1734 A flag indicating how to handle boundaries: 1735 1736 ``fill`` 1737 pad input arrays with fillvalue. (default) 1738 ``wrap`` 1739 circular boundary conditions. 1740 ``symm`` 1741 symmetrical boundary conditions. 1742 1743 fillvalue : scalar, optional 1744 Value to fill pad input arrays with. Default is 0. 1745 1746 Returns 1747 ------- 1748 correlate2d : ndarray 1749 A 2-dimensional array containing a subset of the discrete linear 1750 cross-correlation of `in1` with `in2`. 1751 1752 Notes 1753 ----- 1754 When using "same" mode with even-length inputs, the outputs of `correlate` 1755 and `correlate2d` differ: There is a 1-index offset between them. 1756 1757 Examples 1758 -------- 1759 Use 2D cross-correlation to find the location of a template in a noisy 1760 image: 1761 1762 >>> from scipy import signal 1763 >>> from scipy import misc 1764 >>> rng = np.random.default_rng() 1765 >>> face = misc.face(gray=True) - misc.face(gray=True).mean() 1766 >>> template = np.copy(face[300:365, 670:750]) # right eye 1767 >>> template -= template.mean() 1768 >>> face = face + rng.standard_normal(face.shape) * 50 # add noise 1769 >>> corr = signal.correlate2d(face, template, boundary='symm', mode='same') 1770 >>> y, x = np.unravel_index(np.argmax(corr), corr.shape) # find the match 1771 1772 >>> import matplotlib.pyplot as plt 1773 >>> fig, (ax_orig, ax_template, ax_corr) = plt.subplots(3, 1, 1774 ... figsize=(6, 15)) 1775 >>> ax_orig.imshow(face, cmap='gray') 1776 >>> ax_orig.set_title('Original') 1777 >>> ax_orig.set_axis_off() 1778 >>> ax_template.imshow(template, cmap='gray') 1779 >>> ax_template.set_title('Template') 1780 >>> ax_template.set_axis_off() 1781 >>> ax_corr.imshow(corr, cmap='gray') 1782 >>> ax_corr.set_title('Cross-correlation') 1783 >>> ax_corr.set_axis_off() 1784 >>> ax_orig.plot(x, y, 'ro') 1785 >>> fig.show() 1786 1787 """ 1788 in1 = np.asarray(in1) 1789 in2 = np.asarray(in2) 1790 1791 if not in1.ndim == in2.ndim == 2: 1792 raise ValueError('correlate2d inputs must both be 2-D arrays') 1793 1794 swapped_inputs = _inputs_swap_needed(mode, in1.shape, in2.shape) 1795 if swapped_inputs: 1796 in1, in2 = in2, in1 1797 1798 val = _valfrommode(mode) 1799 bval = _bvalfromboundary(boundary) 1800 out = sigtools._convolve2d(in1, in2.conj(), 0, val, bval, fillvalue) 1801 1802 if swapped_inputs: 1803 out = out[::-1, ::-1] 1804 1805 return out 1806 1807 1808def medfilt2d(input, kernel_size=3): 1809 """ 1810 Median filter a 2-dimensional array. 1811 1812 Apply a median filter to the `input` array using a local window-size 1813 given by `kernel_size` (must be odd). The array is zero-padded 1814 automatically. 1815 1816 Parameters 1817 ---------- 1818 input : array_like 1819 A 2-dimensional input array. 1820 kernel_size : array_like, optional 1821 A scalar or a list of length 2, giving the size of the 1822 median filter window in each dimension. Elements of 1823 `kernel_size` should be odd. If `kernel_size` is a scalar, 1824 then this scalar is used as the size in each dimension. 1825 Default is a kernel of size (3, 3). 1826 1827 Returns 1828 ------- 1829 out : ndarray 1830 An array the same size as input containing the median filtered 1831 result. 1832 1833 See also 1834 -------- 1835 scipy.ndimage.median_filter 1836 1837 Notes 1838 ----- 1839 This is faster than `medfilt` when the input dtype is ``uint8``, 1840 ``float32``, or ``float64``; for other types, this falls back to 1841 `medfilt`; you should use `scipy.ndimage.median_filter` instead as it is 1842 much faster. In some situations, `scipy.ndimage.median_filter` may be 1843 faster than this function. 1844 1845 """ 1846 image = np.asarray(input) 1847 1848 # checking dtype.type, rather than just dtype, is necessary for 1849 # excluding np.longdouble with MS Visual C. 1850 if image.dtype.type not in (np.ubyte, np.single, np.double): 1851 return medfilt(image, kernel_size) 1852 1853 if kernel_size is None: 1854 kernel_size = [3] * 2 1855 kernel_size = np.asarray(kernel_size) 1856 if kernel_size.shape == (): 1857 kernel_size = np.repeat(kernel_size.item(), 2) 1858 1859 for size in kernel_size: 1860 if (size % 2) != 1: 1861 raise ValueError("Each element of kernel_size should be odd.") 1862 1863 return sigtools._medfilt2d(image, kernel_size) 1864 1865 1866def lfilter(b, a, x, axis=-1, zi=None): 1867 """ 1868 Filter data along one-dimension with an IIR or FIR filter. 1869 1870 Filter a data sequence, `x`, using a digital filter. This works for many 1871 fundamental data types (including Object type). The filter is a direct 1872 form II transposed implementation of the standard difference equation 1873 (see Notes). 1874 1875 The function `sosfilt` (and filter design using ``output='sos'``) should be 1876 preferred over `lfilter` for most filtering tasks, as second-order sections 1877 have fewer numerical problems. 1878 1879 Parameters 1880 ---------- 1881 b : array_like 1882 The numerator coefficient vector in a 1-D sequence. 1883 a : array_like 1884 The denominator coefficient vector in a 1-D sequence. If ``a[0]`` 1885 is not 1, then both `a` and `b` are normalized by ``a[0]``. 1886 x : array_like 1887 An N-dimensional input array. 1888 axis : int, optional 1889 The axis of the input data array along which to apply the 1890 linear filter. The filter is applied to each subarray along 1891 this axis. Default is -1. 1892 zi : array_like, optional 1893 Initial conditions for the filter delays. It is a vector 1894 (or array of vectors for an N-dimensional input) of length 1895 ``max(len(a), len(b)) - 1``. If `zi` is None or is not given then 1896 initial rest is assumed. See `lfiltic` for more information. 1897 1898 Returns 1899 ------- 1900 y : array 1901 The output of the digital filter. 1902 zf : array, optional 1903 If `zi` is None, this is not returned, otherwise, `zf` holds the 1904 final filter delay values. 1905 1906 See Also 1907 -------- 1908 lfiltic : Construct initial conditions for `lfilter`. 1909 lfilter_zi : Compute initial state (steady state of step response) for 1910 `lfilter`. 1911 filtfilt : A forward-backward filter, to obtain a filter with linear phase. 1912 savgol_filter : A Savitzky-Golay filter. 1913 sosfilt: Filter data using cascaded second-order sections. 1914 sosfiltfilt: A forward-backward filter using second-order sections. 1915 1916 Notes 1917 ----- 1918 The filter function is implemented as a direct II transposed structure. 1919 This means that the filter implements:: 1920 1921 a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M] 1922 - a[1]*y[n-1] - ... - a[N]*y[n-N] 1923 1924 where `M` is the degree of the numerator, `N` is the degree of the 1925 denominator, and `n` is the sample number. It is implemented using 1926 the following difference equations (assuming M = N):: 1927 1928 a[0]*y[n] = b[0] * x[n] + d[0][n-1] 1929 d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n-1] 1930 d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n-1] 1931 ... 1932 d[N-2][n] = b[N-1]*x[n] - a[N-1]*y[n] + d[N-1][n-1] 1933 d[N-1][n] = b[N] * x[n] - a[N] * y[n] 1934 1935 where `d` are the state variables. 1936 1937 The rational transfer function describing this filter in the 1938 z-transform domain is:: 1939 1940 -1 -M 1941 b[0] + b[1]z + ... + b[M] z 1942 Y(z) = -------------------------------- X(z) 1943 -1 -N 1944 a[0] + a[1]z + ... + a[N] z 1945 1946 Examples 1947 -------- 1948 Generate a noisy signal to be filtered: 1949 1950 >>> from scipy import signal 1951 >>> import matplotlib.pyplot as plt 1952 >>> rng = np.random.default_rng() 1953 >>> t = np.linspace(-1, 1, 201) 1954 >>> x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) + 1955 ... 0.1*np.sin(2*np.pi*1.25*t + 1) + 1956 ... 0.18*np.cos(2*np.pi*3.85*t)) 1957 >>> xn = x + rng.standard_normal(len(t)) * 0.08 1958 1959 Create an order 3 lowpass butterworth filter: 1960 1961 >>> b, a = signal.butter(3, 0.05) 1962 1963 Apply the filter to xn. Use lfilter_zi to choose the initial condition of 1964 the filter: 1965 1966 >>> zi = signal.lfilter_zi(b, a) 1967 >>> z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0]) 1968 1969 Apply the filter again, to have a result filtered at an order the same as 1970 filtfilt: 1971 1972 >>> z2, _ = signal.lfilter(b, a, z, zi=zi*z[0]) 1973 1974 Use filtfilt to apply the filter: 1975 1976 >>> y = signal.filtfilt(b, a, xn) 1977 1978 Plot the original signal and the various filtered versions: 1979 1980 >>> plt.figure 1981 >>> plt.plot(t, xn, 'b', alpha=0.75) 1982 >>> plt.plot(t, z, 'r--', t, z2, 'r', t, y, 'k') 1983 >>> plt.legend(('noisy signal', 'lfilter, once', 'lfilter, twice', 1984 ... 'filtfilt'), loc='best') 1985 >>> plt.grid(True) 1986 >>> plt.show() 1987 1988 """ 1989 a = np.atleast_1d(a) 1990 if len(a) == 1: 1991 # This path only supports types fdgFDGO to mirror _linear_filter below. 1992 # Any of b, a, x, or zi can set the dtype, but there is no default 1993 # casting of other types; instead a NotImplementedError is raised. 1994 b = np.asarray(b) 1995 a = np.asarray(a) 1996 if b.ndim != 1 and a.ndim != 1: 1997 raise ValueError('object of too small depth for desired array') 1998 x = _validate_x(x) 1999 inputs = [b, a, x] 2000 if zi is not None: 2001 # _linear_filter does not broadcast zi, but does do expansion of 2002 # singleton dims. 2003 zi = np.asarray(zi) 2004 if zi.ndim != x.ndim: 2005 raise ValueError('object of too small depth for desired array') 2006 expected_shape = list(x.shape) 2007 expected_shape[axis] = b.shape[0] - 1 2008 expected_shape = tuple(expected_shape) 2009 # check the trivial case where zi is the right shape first 2010 if zi.shape != expected_shape: 2011 strides = zi.ndim * [None] 2012 if axis < 0: 2013 axis += zi.ndim 2014 for k in range(zi.ndim): 2015 if k == axis and zi.shape[k] == expected_shape[k]: 2016 strides[k] = zi.strides[k] 2017 elif k != axis and zi.shape[k] == expected_shape[k]: 2018 strides[k] = zi.strides[k] 2019 elif k != axis and zi.shape[k] == 1: 2020 strides[k] = 0 2021 else: 2022 raise ValueError('Unexpected shape for zi: expected ' 2023 '%s, found %s.' % 2024 (expected_shape, zi.shape)) 2025 zi = np.lib.stride_tricks.as_strided(zi, expected_shape, 2026 strides) 2027 inputs.append(zi) 2028 dtype = np.result_type(*inputs) 2029 2030 if dtype.char not in 'fdgFDGO': 2031 raise NotImplementedError("input type '%s' not supported" % dtype) 2032 2033 b = np.array(b, dtype=dtype) 2034 a = np.array(a, dtype=dtype, copy=False) 2035 b /= a[0] 2036 x = np.array(x, dtype=dtype, copy=False) 2037 2038 out_full = np.apply_along_axis(lambda y: np.convolve(b, y), axis, x) 2039 ind = out_full.ndim * [slice(None)] 2040 if zi is not None: 2041 ind[axis] = slice(zi.shape[axis]) 2042 out_full[tuple(ind)] += zi 2043 2044 ind[axis] = slice(out_full.shape[axis] - len(b) + 1) 2045 out = out_full[tuple(ind)] 2046 2047 if zi is None: 2048 return out 2049 else: 2050 ind[axis] = slice(out_full.shape[axis] - len(b) + 1, None) 2051 zf = out_full[tuple(ind)] 2052 return out, zf 2053 else: 2054 if zi is None: 2055 return sigtools._linear_filter(b, a, x, axis) 2056 else: 2057 return sigtools._linear_filter(b, a, x, axis, zi) 2058 2059 2060def lfiltic(b, a, y, x=None): 2061 """ 2062 Construct initial conditions for lfilter given input and output vectors. 2063 2064 Given a linear filter (b, a) and initial conditions on the output `y` 2065 and the input `x`, return the initial conditions on the state vector zi 2066 which is used by `lfilter` to generate the output given the input. 2067 2068 Parameters 2069 ---------- 2070 b : array_like 2071 Linear filter term. 2072 a : array_like 2073 Linear filter term. 2074 y : array_like 2075 Initial conditions. 2076 2077 If ``N = len(a) - 1``, then ``y = {y[-1], y[-2], ..., y[-N]}``. 2078 2079 If `y` is too short, it is padded with zeros. 2080 x : array_like, optional 2081 Initial conditions. 2082 2083 If ``M = len(b) - 1``, then ``x = {x[-1], x[-2], ..., x[-M]}``. 2084 2085 If `x` is not given, its initial conditions are assumed zero. 2086 2087 If `x` is too short, it is padded with zeros. 2088 2089 Returns 2090 ------- 2091 zi : ndarray 2092 The state vector ``zi = {z_0[-1], z_1[-1], ..., z_K-1[-1]}``, 2093 where ``K = max(M, N)``. 2094 2095 See Also 2096 -------- 2097 lfilter, lfilter_zi 2098 2099 """ 2100 N = np.size(a) - 1 2101 M = np.size(b) - 1 2102 K = max(M, N) 2103 y = np.asarray(y) 2104 2105 if x is None: 2106 result_type = np.result_type(np.asarray(b), np.asarray(a), y) 2107 if result_type.kind in 'bui': 2108 result_type = np.float64 2109 x = np.zeros(M, dtype=result_type) 2110 else: 2111 x = np.asarray(x) 2112 2113 result_type = np.result_type(np.asarray(b), np.asarray(a), y, x) 2114 if result_type.kind in 'bui': 2115 result_type = np.float64 2116 x = x.astype(result_type) 2117 2118 L = np.size(x) 2119 if L < M: 2120 x = np.r_[x, np.zeros(M - L)] 2121 2122 y = y.astype(result_type) 2123 zi = np.zeros(K, result_type) 2124 2125 L = np.size(y) 2126 if L < N: 2127 y = np.r_[y, np.zeros(N - L)] 2128 2129 for m in range(M): 2130 zi[m] = np.sum(b[m + 1:] * x[:M - m], axis=0) 2131 2132 for m in range(N): 2133 zi[m] -= np.sum(a[m + 1:] * y[:N - m], axis=0) 2134 2135 return zi 2136 2137 2138def deconvolve(signal, divisor): 2139 """Deconvolves ``divisor`` out of ``signal`` using inverse filtering. 2140 2141 Returns the quotient and remainder such that 2142 ``signal = convolve(divisor, quotient) + remainder`` 2143 2144 Parameters 2145 ---------- 2146 signal : array_like 2147 Signal data, typically a recorded signal 2148 divisor : array_like 2149 Divisor data, typically an impulse response or filter that was 2150 applied to the original signal 2151 2152 Returns 2153 ------- 2154 quotient : ndarray 2155 Quotient, typically the recovered original signal 2156 remainder : ndarray 2157 Remainder 2158 2159 Examples 2160 -------- 2161 Deconvolve a signal that's been filtered: 2162 2163 >>> from scipy import signal 2164 >>> original = [0, 1, 0, 0, 1, 1, 0, 0] 2165 >>> impulse_response = [2, 1] 2166 >>> recorded = signal.convolve(impulse_response, original) 2167 >>> recorded 2168 array([0, 2, 1, 0, 2, 3, 1, 0, 0]) 2169 >>> recovered, remainder = signal.deconvolve(recorded, impulse_response) 2170 >>> recovered 2171 array([ 0., 1., 0., 0., 1., 1., 0., 0.]) 2172 2173 See Also 2174 -------- 2175 numpy.polydiv : performs polynomial division (same operation, but 2176 also accepts poly1d objects) 2177 2178 """ 2179 num = np.atleast_1d(signal) 2180 den = np.atleast_1d(divisor) 2181 N = len(num) 2182 D = len(den) 2183 if D > N: 2184 quot = [] 2185 rem = num 2186 else: 2187 input = np.zeros(N - D + 1, float) 2188 input[0] = 1 2189 quot = lfilter(num, den, input) 2190 rem = num - convolve(den, quot, mode='full') 2191 return quot, rem 2192 2193 2194def hilbert(x, N=None, axis=-1): 2195 """ 2196 Compute the analytic signal, using the Hilbert transform. 2197 2198 The transformation is done along the last axis by default. 2199 2200 Parameters 2201 ---------- 2202 x : array_like 2203 Signal data. Must be real. 2204 N : int, optional 2205 Number of Fourier components. Default: ``x.shape[axis]`` 2206 axis : int, optional 2207 Axis along which to do the transformation. Default: -1. 2208 2209 Returns 2210 ------- 2211 xa : ndarray 2212 Analytic signal of `x`, of each 1-D array along `axis` 2213 2214 Notes 2215 ----- 2216 The analytic signal ``x_a(t)`` of signal ``x(t)`` is: 2217 2218 .. math:: x_a = F^{-1}(F(x) 2U) = x + i y 2219 2220 where `F` is the Fourier transform, `U` the unit step function, 2221 and `y` the Hilbert transform of `x`. [1]_ 2222 2223 In other words, the negative half of the frequency spectrum is zeroed 2224 out, turning the real-valued signal into a complex signal. The Hilbert 2225 transformed signal can be obtained from ``np.imag(hilbert(x))``, and the 2226 original signal from ``np.real(hilbert(x))``. 2227 2228 Examples 2229 -------- 2230 In this example we use the Hilbert transform to determine the amplitude 2231 envelope and instantaneous frequency of an amplitude-modulated signal. 2232 2233 >>> import numpy as np 2234 >>> import matplotlib.pyplot as plt 2235 >>> from scipy.signal import hilbert, chirp 2236 2237 >>> duration = 1.0 2238 >>> fs = 400.0 2239 >>> samples = int(fs*duration) 2240 >>> t = np.arange(samples) / fs 2241 2242 We create a chirp of which the frequency increases from 20 Hz to 100 Hz and 2243 apply an amplitude modulation. 2244 2245 >>> signal = chirp(t, 20.0, t[-1], 100.0) 2246 >>> signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) ) 2247 2248 The amplitude envelope is given by magnitude of the analytic signal. The 2249 instantaneous frequency can be obtained by differentiating the 2250 instantaneous phase in respect to time. The instantaneous phase corresponds 2251 to the phase angle of the analytic signal. 2252 2253 >>> analytic_signal = hilbert(signal) 2254 >>> amplitude_envelope = np.abs(analytic_signal) 2255 >>> instantaneous_phase = np.unwrap(np.angle(analytic_signal)) 2256 >>> instantaneous_frequency = (np.diff(instantaneous_phase) / 2257 ... (2.0*np.pi) * fs) 2258 2259 >>> fig, (ax0, ax1) = plt.subplots(nrows=2) 2260 >>> ax0.plot(t, signal, label='signal') 2261 >>> ax0.plot(t, amplitude_envelope, label='envelope') 2262 >>> ax0.set_xlabel("time in seconds") 2263 >>> ax0.legend() 2264 >>> ax1.plot(t[1:], instantaneous_frequency) 2265 >>> ax1.set_xlabel("time in seconds") 2266 >>> ax1.set_ylim(0.0, 120.0) 2267 >>> fig.tight_layout() 2268 2269 References 2270 ---------- 2271 .. [1] Wikipedia, "Analytic signal". 2272 https://en.wikipedia.org/wiki/Analytic_signal 2273 .. [2] Leon Cohen, "Time-Frequency Analysis", 1995. Chapter 2. 2274 .. [3] Alan V. Oppenheim, Ronald W. Schafer. Discrete-Time Signal 2275 Processing, Third Edition, 2009. Chapter 12. 2276 ISBN 13: 978-1292-02572-8 2277 2278 """ 2279 x = np.asarray(x) 2280 if np.iscomplexobj(x): 2281 raise ValueError("x must be real.") 2282 if N is None: 2283 N = x.shape[axis] 2284 if N <= 0: 2285 raise ValueError("N must be positive.") 2286 2287 Xf = sp_fft.fft(x, N, axis=axis) 2288 h = np.zeros(N) 2289 if N % 2 == 0: 2290 h[0] = h[N // 2] = 1 2291 h[1:N // 2] = 2 2292 else: 2293 h[0] = 1 2294 h[1:(N + 1) // 2] = 2 2295 2296 if x.ndim > 1: 2297 ind = [np.newaxis] * x.ndim 2298 ind[axis] = slice(None) 2299 h = h[tuple(ind)] 2300 x = sp_fft.ifft(Xf * h, axis=axis) 2301 return x 2302 2303 2304def hilbert2(x, N=None): 2305 """ 2306 Compute the '2-D' analytic signal of `x` 2307 2308 Parameters 2309 ---------- 2310 x : array_like 2311 2-D signal data. 2312 N : int or tuple of two ints, optional 2313 Number of Fourier components. Default is ``x.shape`` 2314 2315 Returns 2316 ------- 2317 xa : ndarray 2318 Analytic signal of `x` taken along axes (0,1). 2319 2320 References 2321 ---------- 2322 .. [1] Wikipedia, "Analytic signal", 2323 https://en.wikipedia.org/wiki/Analytic_signal 2324 2325 """ 2326 x = np.atleast_2d(x) 2327 if x.ndim > 2: 2328 raise ValueError("x must be 2-D.") 2329 if np.iscomplexobj(x): 2330 raise ValueError("x must be real.") 2331 if N is None: 2332 N = x.shape 2333 elif isinstance(N, int): 2334 if N <= 0: 2335 raise ValueError("N must be positive.") 2336 N = (N, N) 2337 elif len(N) != 2 or np.any(np.asarray(N) <= 0): 2338 raise ValueError("When given as a tuple, N must hold exactly " 2339 "two positive integers") 2340 2341 Xf = sp_fft.fft2(x, N, axes=(0, 1)) 2342 h1 = np.zeros(N[0], 'd') 2343 h2 = np.zeros(N[1], 'd') 2344 for p in range(2): 2345 h = eval("h%d" % (p + 1)) 2346 N1 = N[p] 2347 if N1 % 2 == 0: 2348 h[0] = h[N1 // 2] = 1 2349 h[1:N1 // 2] = 2 2350 else: 2351 h[0] = 1 2352 h[1:(N1 + 1) // 2] = 2 2353 exec("h%d = h" % (p + 1), globals(), locals()) 2354 2355 h = h1[:, np.newaxis] * h2[np.newaxis, :] 2356 k = x.ndim 2357 while k > 2: 2358 h = h[:, np.newaxis] 2359 k -= 1 2360 x = sp_fft.ifft2(Xf * h, axes=(0, 1)) 2361 return x 2362 2363 2364def cmplx_sort(p): 2365 """Sort roots based on magnitude. 2366 2367 Parameters 2368 ---------- 2369 p : array_like 2370 The roots to sort, as a 1-D array. 2371 2372 Returns 2373 ------- 2374 p_sorted : ndarray 2375 Sorted roots. 2376 indx : ndarray 2377 Array of indices needed to sort the input `p`. 2378 2379 Examples 2380 -------- 2381 >>> from scipy import signal 2382 >>> vals = [1, 4, 1+1.j, 3] 2383 >>> p_sorted, indx = signal.cmplx_sort(vals) 2384 >>> p_sorted 2385 array([1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j]) 2386 >>> indx 2387 array([0, 2, 3, 1]) 2388 """ 2389 p = np.asarray(p) 2390 indx = np.argsort(abs(p)) 2391 return np.take(p, indx, 0), indx 2392 2393 2394def unique_roots(p, tol=1e-3, rtype='min'): 2395 """Determine unique roots and their multiplicities from a list of roots. 2396 2397 Parameters 2398 ---------- 2399 p : array_like 2400 The list of roots. 2401 tol : float, optional 2402 The tolerance for two roots to be considered equal in terms of 2403 the distance between them. Default is 1e-3. Refer to Notes about 2404 the details on roots grouping. 2405 rtype : {'max', 'maximum', 'min', 'minimum', 'avg', 'mean'}, optional 2406 How to determine the returned root if multiple roots are within 2407 `tol` of each other. 2408 2409 - 'max', 'maximum': pick the maximum of those roots 2410 - 'min', 'minimum': pick the minimum of those roots 2411 - 'avg', 'mean': take the average of those roots 2412 2413 When finding minimum or maximum among complex roots they are compared 2414 first by the real part and then by the imaginary part. 2415 2416 Returns 2417 ------- 2418 unique : ndarray 2419 The list of unique roots. 2420 multiplicity : ndarray 2421 The multiplicity of each root. 2422 2423 Notes 2424 ----- 2425 If we have 3 roots ``a``, ``b`` and ``c``, such that ``a`` is close to 2426 ``b`` and ``b`` is close to ``c`` (distance is less than `tol`), then it 2427 doesn't necessarily mean that ``a`` is close to ``c``. It means that roots 2428 grouping is not unique. In this function we use "greedy" grouping going 2429 through the roots in the order they are given in the input `p`. 2430 2431 This utility function is not specific to roots but can be used for any 2432 sequence of values for which uniqueness and multiplicity has to be 2433 determined. For a more general routine, see `numpy.unique`. 2434 2435 Examples 2436 -------- 2437 >>> from scipy import signal 2438 >>> vals = [0, 1.3, 1.31, 2.8, 1.25, 2.2, 10.3] 2439 >>> uniq, mult = signal.unique_roots(vals, tol=2e-2, rtype='avg') 2440 2441 Check which roots have multiplicity larger than 1: 2442 2443 >>> uniq[mult > 1] 2444 array([ 1.305]) 2445 """ 2446 if rtype in ['max', 'maximum']: 2447 reduce = np.max 2448 elif rtype in ['min', 'minimum']: 2449 reduce = np.min 2450 elif rtype in ['avg', 'mean']: 2451 reduce = np.mean 2452 else: 2453 raise ValueError("`rtype` must be one of " 2454 "{'max', 'maximum', 'min', 'minimum', 'avg', 'mean'}") 2455 2456 p = np.asarray(p) 2457 2458 points = np.empty((len(p), 2)) 2459 points[:, 0] = np.real(p) 2460 points[:, 1] = np.imag(p) 2461 tree = cKDTree(points) 2462 2463 p_unique = [] 2464 p_multiplicity = [] 2465 used = np.zeros(len(p), dtype=bool) 2466 for i in range(len(p)): 2467 if used[i]: 2468 continue 2469 2470 group = tree.query_ball_point(points[i], tol) 2471 group = [x for x in group if not used[x]] 2472 2473 p_unique.append(reduce(p[group])) 2474 p_multiplicity.append(len(group)) 2475 2476 used[group] = True 2477 2478 return np.asarray(p_unique), np.asarray(p_multiplicity) 2479 2480 2481def invres(r, p, k, tol=1e-3, rtype='avg'): 2482 """Compute b(s) and a(s) from partial fraction expansion. 2483 2484 If `M` is the degree of numerator `b` and `N` the degree of denominator 2485 `a`:: 2486 2487 b(s) b[0] s**(M) + b[1] s**(M-1) + ... + b[M] 2488 H(s) = ------ = ------------------------------------------ 2489 a(s) a[0] s**(N) + a[1] s**(N-1) + ... + a[N] 2490 2491 then the partial-fraction expansion H(s) is defined as:: 2492 2493 r[0] r[1] r[-1] 2494 = -------- + -------- + ... + --------- + k(s) 2495 (s-p[0]) (s-p[1]) (s-p[-1]) 2496 2497 If there are any repeated roots (closer together than `tol`), then H(s) 2498 has terms like:: 2499 2500 r[i] r[i+1] r[i+n-1] 2501 -------- + ----------- + ... + ----------- 2502 (s-p[i]) (s-p[i])**2 (s-p[i])**n 2503 2504 This function is used for polynomials in positive powers of s or z, 2505 such as analog filters or digital filters in controls engineering. For 2506 negative powers of z (typical for digital filters in DSP), use `invresz`. 2507 2508 Parameters 2509 ---------- 2510 r : array_like 2511 Residues corresponding to the poles. For repeated poles, the residues 2512 must be ordered to correspond to ascending by power fractions. 2513 p : array_like 2514 Poles. Equal poles must be adjacent. 2515 k : array_like 2516 Coefficients of the direct polynomial term. 2517 tol : float, optional 2518 The tolerance for two roots to be considered equal in terms of 2519 the distance between them. Default is 1e-3. See `unique_roots` 2520 for further details. 2521 rtype : {'avg', 'min', 'max'}, optional 2522 Method for computing a root to represent a group of identical roots. 2523 Default is 'avg'. See `unique_roots` for further details. 2524 2525 Returns 2526 ------- 2527 b : ndarray 2528 Numerator polynomial coefficients. 2529 a : ndarray 2530 Denominator polynomial coefficients. 2531 2532 See Also 2533 -------- 2534 residue, invresz, unique_roots 2535 2536 """ 2537 r = np.atleast_1d(r) 2538 p = np.atleast_1d(p) 2539 k = np.trim_zeros(np.atleast_1d(k), 'f') 2540 2541 unique_poles, multiplicity = _group_poles(p, tol, rtype) 2542 factors, denominator = _compute_factors(unique_poles, multiplicity, 2543 include_powers=True) 2544 2545 if len(k) == 0: 2546 numerator = 0 2547 else: 2548 numerator = np.polymul(k, denominator) 2549 2550 for residue, factor in zip(r, factors): 2551 numerator = np.polyadd(numerator, residue * factor) 2552 2553 return numerator, denominator 2554 2555 2556def _compute_factors(roots, multiplicity, include_powers=False): 2557 """Compute the total polynomial divided by factors for each root.""" 2558 current = np.array([1]) 2559 suffixes = [current] 2560 for pole, mult in zip(roots[-1:0:-1], multiplicity[-1:0:-1]): 2561 monomial = np.array([1, -pole]) 2562 for _ in range(mult): 2563 current = np.polymul(current, monomial) 2564 suffixes.append(current) 2565 suffixes = suffixes[::-1] 2566 2567 factors = [] 2568 current = np.array([1]) 2569 for pole, mult, suffix in zip(roots, multiplicity, suffixes): 2570 monomial = np.array([1, -pole]) 2571 block = [] 2572 for i in range(mult): 2573 if i == 0 or include_powers: 2574 block.append(np.polymul(current, suffix)) 2575 current = np.polymul(current, monomial) 2576 factors.extend(reversed(block)) 2577 2578 return factors, current 2579 2580 2581def _compute_residues(poles, multiplicity, numerator): 2582 denominator_factors, _ = _compute_factors(poles, multiplicity) 2583 numerator = numerator.astype(poles.dtype) 2584 2585 residues = [] 2586 for pole, mult, factor in zip(poles, multiplicity, 2587 denominator_factors): 2588 if mult == 1: 2589 residues.append(np.polyval(numerator, pole) / 2590 np.polyval(factor, pole)) 2591 else: 2592 numer = numerator.copy() 2593 monomial = np.array([1, -pole]) 2594 factor, d = np.polydiv(factor, monomial) 2595 2596 block = [] 2597 for _ in range(mult): 2598 numer, n = np.polydiv(numer, monomial) 2599 r = n[0] / d[0] 2600 numer = np.polysub(numer, r * factor) 2601 block.append(r) 2602 2603 residues.extend(reversed(block)) 2604 2605 return np.asarray(residues) 2606 2607 2608def residue(b, a, tol=1e-3, rtype='avg'): 2609 """Compute partial-fraction expansion of b(s) / a(s). 2610 2611 If `M` is the degree of numerator `b` and `N` the degree of denominator 2612 `a`:: 2613 2614 b(s) b[0] s**(M) + b[1] s**(M-1) + ... + b[M] 2615 H(s) = ------ = ------------------------------------------ 2616 a(s) a[0] s**(N) + a[1] s**(N-1) + ... + a[N] 2617 2618 then the partial-fraction expansion H(s) is defined as:: 2619 2620 r[0] r[1] r[-1] 2621 = -------- + -------- + ... + --------- + k(s) 2622 (s-p[0]) (s-p[1]) (s-p[-1]) 2623 2624 If there are any repeated roots (closer together than `tol`), then H(s) 2625 has terms like:: 2626 2627 r[i] r[i+1] r[i+n-1] 2628 -------- + ----------- + ... + ----------- 2629 (s-p[i]) (s-p[i])**2 (s-p[i])**n 2630 2631 This function is used for polynomials in positive powers of s or z, 2632 such as analog filters or digital filters in controls engineering. For 2633 negative powers of z (typical for digital filters in DSP), use `residuez`. 2634 2635 See Notes for details about the algorithm. 2636 2637 Parameters 2638 ---------- 2639 b : array_like 2640 Numerator polynomial coefficients. 2641 a : array_like 2642 Denominator polynomial coefficients. 2643 tol : float, optional 2644 The tolerance for two roots to be considered equal in terms of 2645 the distance between them. Default is 1e-3. See `unique_roots` 2646 for further details. 2647 rtype : {'avg', 'min', 'max'}, optional 2648 Method for computing a root to represent a group of identical roots. 2649 Default is 'avg'. See `unique_roots` for further details. 2650 2651 Returns 2652 ------- 2653 r : ndarray 2654 Residues corresponding to the poles. For repeated poles, the residues 2655 are ordered to correspond to ascending by power fractions. 2656 p : ndarray 2657 Poles ordered by magnitude in ascending order. 2658 k : ndarray 2659 Coefficients of the direct polynomial term. 2660 2661 See Also 2662 -------- 2663 invres, residuez, numpy.poly, unique_roots 2664 2665 Notes 2666 ----- 2667 The "deflation through subtraction" algorithm is used for 2668 computations --- method 6 in [1]_. 2669 2670 The form of partial fraction expansion depends on poles multiplicity in 2671 the exact mathematical sense. However there is no way to exactly 2672 determine multiplicity of roots of a polynomial in numerical computing. 2673 Thus you should think of the result of `residue` with given `tol` as 2674 partial fraction expansion computed for the denominator composed of the 2675 computed poles with empirically determined multiplicity. The choice of 2676 `tol` can drastically change the result if there are close poles. 2677 2678 References 2679 ---------- 2680 .. [1] J. F. Mahoney, B. D. Sivazlian, "Partial fractions expansion: a 2681 review of computational methodology and efficiency", Journal of 2682 Computational and Applied Mathematics, Vol. 9, 1983. 2683 """ 2684 b = np.asarray(b) 2685 a = np.asarray(a) 2686 if (np.issubdtype(b.dtype, np.complexfloating) 2687 or np.issubdtype(a.dtype, np.complexfloating)): 2688 b = b.astype(complex) 2689 a = a.astype(complex) 2690 else: 2691 b = b.astype(float) 2692 a = a.astype(float) 2693 2694 b = np.trim_zeros(np.atleast_1d(b), 'f') 2695 a = np.trim_zeros(np.atleast_1d(a), 'f') 2696 2697 if a.size == 0: 2698 raise ValueError("Denominator `a` is zero.") 2699 2700 poles = np.roots(a) 2701 if b.size == 0: 2702 return np.zeros(poles.shape), cmplx_sort(poles)[0], np.array([]) 2703 2704 if len(b) < len(a): 2705 k = np.empty(0) 2706 else: 2707 k, b = np.polydiv(b, a) 2708 2709 unique_poles, multiplicity = unique_roots(poles, tol=tol, rtype=rtype) 2710 unique_poles, order = cmplx_sort(unique_poles) 2711 multiplicity = multiplicity[order] 2712 2713 residues = _compute_residues(unique_poles, multiplicity, b) 2714 2715 index = 0 2716 for pole, mult in zip(unique_poles, multiplicity): 2717 poles[index:index + mult] = pole 2718 index += mult 2719 2720 return residues / a[0], poles, k 2721 2722 2723def residuez(b, a, tol=1e-3, rtype='avg'): 2724 """Compute partial-fraction expansion of b(z) / a(z). 2725 2726 If `M` is the degree of numerator `b` and `N` the degree of denominator 2727 `a`:: 2728 2729 b(z) b[0] + b[1] z**(-1) + ... + b[M] z**(-M) 2730 H(z) = ------ = ------------------------------------------ 2731 a(z) a[0] + a[1] z**(-1) + ... + a[N] z**(-N) 2732 2733 then the partial-fraction expansion H(z) is defined as:: 2734 2735 r[0] r[-1] 2736 = --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ... 2737 (1-p[0]z**(-1)) (1-p[-1]z**(-1)) 2738 2739 If there are any repeated roots (closer than `tol`), then the partial 2740 fraction expansion has terms like:: 2741 2742 r[i] r[i+1] r[i+n-1] 2743 -------------- + ------------------ + ... + ------------------ 2744 (1-p[i]z**(-1)) (1-p[i]z**(-1))**2 (1-p[i]z**(-1))**n 2745 2746 This function is used for polynomials in negative powers of z, 2747 such as digital filters in DSP. For positive powers, use `residue`. 2748 2749 See Notes of `residue` for details about the algorithm. 2750 2751 Parameters 2752 ---------- 2753 b : array_like 2754 Numerator polynomial coefficients. 2755 a : array_like 2756 Denominator polynomial coefficients. 2757 tol : float, optional 2758 The tolerance for two roots to be considered equal in terms of 2759 the distance between them. Default is 1e-3. See `unique_roots` 2760 for further details. 2761 rtype : {'avg', 'min', 'max'}, optional 2762 Method for computing a root to represent a group of identical roots. 2763 Default is 'avg'. See `unique_roots` for further details. 2764 2765 Returns 2766 ------- 2767 r : ndarray 2768 Residues corresponding to the poles. For repeated poles, the residues 2769 are ordered to correspond to ascending by power fractions. 2770 p : ndarray 2771 Poles ordered by magnitude in ascending order. 2772 k : ndarray 2773 Coefficients of the direct polynomial term. 2774 2775 See Also 2776 -------- 2777 invresz, residue, unique_roots 2778 """ 2779 b = np.asarray(b) 2780 a = np.asarray(a) 2781 if (np.issubdtype(b.dtype, np.complexfloating) 2782 or np.issubdtype(a.dtype, np.complexfloating)): 2783 b = b.astype(complex) 2784 a = a.astype(complex) 2785 else: 2786 b = b.astype(float) 2787 a = a.astype(float) 2788 2789 b = np.trim_zeros(np.atleast_1d(b), 'b') 2790 a = np.trim_zeros(np.atleast_1d(a), 'b') 2791 2792 if a.size == 0: 2793 raise ValueError("Denominator `a` is zero.") 2794 elif a[0] == 0: 2795 raise ValueError("First coefficient of determinant `a` must be " 2796 "non-zero.") 2797 2798 poles = np.roots(a) 2799 if b.size == 0: 2800 return np.zeros(poles.shape), cmplx_sort(poles)[0], np.array([]) 2801 2802 b_rev = b[::-1] 2803 a_rev = a[::-1] 2804 2805 if len(b_rev) < len(a_rev): 2806 k_rev = np.empty(0) 2807 else: 2808 k_rev, b_rev = np.polydiv(b_rev, a_rev) 2809 2810 unique_poles, multiplicity = unique_roots(poles, tol=tol, rtype=rtype) 2811 unique_poles, order = cmplx_sort(unique_poles) 2812 multiplicity = multiplicity[order] 2813 2814 residues = _compute_residues(1 / unique_poles, multiplicity, b_rev) 2815 2816 index = 0 2817 powers = np.empty(len(residues), dtype=int) 2818 for pole, mult in zip(unique_poles, multiplicity): 2819 poles[index:index + mult] = pole 2820 powers[index:index + mult] = 1 + np.arange(mult) 2821 index += mult 2822 2823 residues *= (-poles) ** powers / a_rev[0] 2824 2825 return residues, poles, k_rev[::-1] 2826 2827 2828def _group_poles(poles, tol, rtype): 2829 if rtype in ['max', 'maximum']: 2830 reduce = np.max 2831 elif rtype in ['min', 'minimum']: 2832 reduce = np.min 2833 elif rtype in ['avg', 'mean']: 2834 reduce = np.mean 2835 else: 2836 raise ValueError("`rtype` must be one of " 2837 "{'max', 'maximum', 'min', 'minimum', 'avg', 'mean'}") 2838 2839 unique = [] 2840 multiplicity = [] 2841 2842 pole = poles[0] 2843 block = [pole] 2844 for i in range(1, len(poles)): 2845 if abs(poles[i] - pole) <= tol: 2846 block.append(pole) 2847 else: 2848 unique.append(reduce(block)) 2849 multiplicity.append(len(block)) 2850 pole = poles[i] 2851 block = [pole] 2852 2853 unique.append(reduce(block)) 2854 multiplicity.append(len(block)) 2855 2856 return np.asarray(unique), np.asarray(multiplicity) 2857 2858 2859def invresz(r, p, k, tol=1e-3, rtype='avg'): 2860 """Compute b(z) and a(z) from partial fraction expansion. 2861 2862 If `M` is the degree of numerator `b` and `N` the degree of denominator 2863 `a`:: 2864 2865 b(z) b[0] + b[1] z**(-1) + ... + b[M] z**(-M) 2866 H(z) = ------ = ------------------------------------------ 2867 a(z) a[0] + a[1] z**(-1) + ... + a[N] z**(-N) 2868 2869 then the partial-fraction expansion H(z) is defined as:: 2870 2871 r[0] r[-1] 2872 = --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ... 2873 (1-p[0]z**(-1)) (1-p[-1]z**(-1)) 2874 2875 If there are any repeated roots (closer than `tol`), then the partial 2876 fraction expansion has terms like:: 2877 2878 r[i] r[i+1] r[i+n-1] 2879 -------------- + ------------------ + ... + ------------------ 2880 (1-p[i]z**(-1)) (1-p[i]z**(-1))**2 (1-p[i]z**(-1))**n 2881 2882 This function is used for polynomials in negative powers of z, 2883 such as digital filters in DSP. For positive powers, use `invres`. 2884 2885 Parameters 2886 ---------- 2887 r : array_like 2888 Residues corresponding to the poles. For repeated poles, the residues 2889 must be ordered to correspond to ascending by power fractions. 2890 p : array_like 2891 Poles. Equal poles must be adjacent. 2892 k : array_like 2893 Coefficients of the direct polynomial term. 2894 tol : float, optional 2895 The tolerance for two roots to be considered equal in terms of 2896 the distance between them. Default is 1e-3. See `unique_roots` 2897 for further details. 2898 rtype : {'avg', 'min', 'max'}, optional 2899 Method for computing a root to represent a group of identical roots. 2900 Default is 'avg'. See `unique_roots` for further details. 2901 2902 Returns 2903 ------- 2904 b : ndarray 2905 Numerator polynomial coefficients. 2906 a : ndarray 2907 Denominator polynomial coefficients. 2908 2909 See Also 2910 -------- 2911 residuez, unique_roots, invres 2912 2913 """ 2914 r = np.atleast_1d(r) 2915 p = np.atleast_1d(p) 2916 k = np.trim_zeros(np.atleast_1d(k), 'b') 2917 2918 unique_poles, multiplicity = _group_poles(p, tol, rtype) 2919 factors, denominator = _compute_factors(unique_poles, multiplicity, 2920 include_powers=True) 2921 2922 if len(k) == 0: 2923 numerator = 0 2924 else: 2925 numerator = np.polymul(k[::-1], denominator[::-1]) 2926 2927 for residue, factor in zip(r, factors): 2928 numerator = np.polyadd(numerator, residue * factor[::-1]) 2929 2930 return numerator[::-1], denominator 2931 2932 2933def resample(x, num, t=None, axis=0, window=None, domain='time'): 2934 """ 2935 Resample `x` to `num` samples using Fourier method along the given axis. 2936 2937 The resampled signal starts at the same value as `x` but is sampled 2938 with a spacing of ``len(x) / num * (spacing of x)``. Because a 2939 Fourier method is used, the signal is assumed to be periodic. 2940 2941 Parameters 2942 ---------- 2943 x : array_like 2944 The data to be resampled. 2945 num : int 2946 The number of samples in the resampled signal. 2947 t : array_like, optional 2948 If `t` is given, it is assumed to be the equally spaced sample 2949 positions associated with the signal data in `x`. 2950 axis : int, optional 2951 The axis of `x` that is resampled. Default is 0. 2952 window : array_like, callable, string, float, or tuple, optional 2953 Specifies the window applied to the signal in the Fourier 2954 domain. See below for details. 2955 domain : string, optional 2956 A string indicating the domain of the input `x`: 2957 ``time`` Consider the input `x` as time-domain (Default), 2958 ``freq`` Consider the input `x` as frequency-domain. 2959 2960 Returns 2961 ------- 2962 resampled_x or (resampled_x, resampled_t) 2963 Either the resampled array, or, if `t` was given, a tuple 2964 containing the resampled array and the corresponding resampled 2965 positions. 2966 2967 See Also 2968 -------- 2969 decimate : Downsample the signal after applying an FIR or IIR filter. 2970 resample_poly : Resample using polyphase filtering and an FIR filter. 2971 2972 Notes 2973 ----- 2974 The argument `window` controls a Fourier-domain window that tapers 2975 the Fourier spectrum before zero-padding to alleviate ringing in 2976 the resampled values for sampled signals you didn't intend to be 2977 interpreted as band-limited. 2978 2979 If `window` is a function, then it is called with a vector of inputs 2980 indicating the frequency bins (i.e. fftfreq(x.shape[axis]) ). 2981 2982 If `window` is an array of the same length as `x.shape[axis]` it is 2983 assumed to be the window to be applied directly in the Fourier 2984 domain (with dc and low-frequency first). 2985 2986 For any other type of `window`, the function `scipy.signal.get_window` 2987 is called to generate the window. 2988 2989 The first sample of the returned vector is the same as the first 2990 sample of the input vector. The spacing between samples is changed 2991 from ``dx`` to ``dx * len(x) / num``. 2992 2993 If `t` is not None, then it is used solely to calculate the resampled 2994 positions `resampled_t` 2995 2996 As noted, `resample` uses FFT transformations, which can be very 2997 slow if the number of input or output samples is large and prime; 2998 see `scipy.fft.fft`. 2999 3000 Examples 3001 -------- 3002 Note that the end of the resampled data rises to meet the first 3003 sample of the next cycle: 3004 3005 >>> from scipy import signal 3006 3007 >>> x = np.linspace(0, 10, 20, endpoint=False) 3008 >>> y = np.cos(-x**2/6.0) 3009 >>> f = signal.resample(y, 100) 3010 >>> xnew = np.linspace(0, 10, 100, endpoint=False) 3011 3012 >>> import matplotlib.pyplot as plt 3013 >>> plt.plot(x, y, 'go-', xnew, f, '.-', 10, y[0], 'ro') 3014 >>> plt.legend(['data', 'resampled'], loc='best') 3015 >>> plt.show() 3016 """ 3017 3018 if domain not in ('time', 'freq'): 3019 raise ValueError("Acceptable domain flags are 'time' or" 3020 " 'freq', not domain={}".format(domain)) 3021 3022 x = np.asarray(x) 3023 Nx = x.shape[axis] 3024 3025 # Check if we can use faster real FFT 3026 real_input = np.isrealobj(x) 3027 3028 if domain == 'time': 3029 # Forward transform 3030 if real_input: 3031 X = sp_fft.rfft(x, axis=axis) 3032 else: # Full complex FFT 3033 X = sp_fft.fft(x, axis=axis) 3034 else: # domain == 'freq' 3035 X = x 3036 3037 # Apply window to spectrum 3038 if window is not None: 3039 if callable(window): 3040 W = window(sp_fft.fftfreq(Nx)) 3041 elif isinstance(window, np.ndarray): 3042 if window.shape != (Nx,): 3043 raise ValueError('window must have the same length as data') 3044 W = window 3045 else: 3046 W = sp_fft.ifftshift(get_window(window, Nx)) 3047 3048 newshape_W = [1] * x.ndim 3049 newshape_W[axis] = X.shape[axis] 3050 if real_input: 3051 # Fold the window back on itself to mimic complex behavior 3052 W_real = W.copy() 3053 W_real[1:] += W_real[-1:0:-1] 3054 W_real[1:] *= 0.5 3055 X *= W_real[:newshape_W[axis]].reshape(newshape_W) 3056 else: 3057 X *= W.reshape(newshape_W) 3058 3059 # Copy each half of the original spectrum to the output spectrum, either 3060 # truncating high frequences (downsampling) or zero-padding them 3061 # (upsampling) 3062 3063 # Placeholder array for output spectrum 3064 newshape = list(x.shape) 3065 if real_input: 3066 newshape[axis] = num // 2 + 1 3067 else: 3068 newshape[axis] = num 3069 Y = np.zeros(newshape, X.dtype) 3070 3071 # Copy positive frequency components (and Nyquist, if present) 3072 N = min(num, Nx) 3073 nyq = N // 2 + 1 # Slice index that includes Nyquist if present 3074 sl = [slice(None)] * x.ndim 3075 sl[axis] = slice(0, nyq) 3076 Y[tuple(sl)] = X[tuple(sl)] 3077 if not real_input: 3078 # Copy negative frequency components 3079 if N > 2: # (slice expression doesn't collapse to empty array) 3080 sl[axis] = slice(nyq - N, None) 3081 Y[tuple(sl)] = X[tuple(sl)] 3082 3083 # Split/join Nyquist component(s) if present 3084 # So far we have set Y[+N/2]=X[+N/2] 3085 if N % 2 == 0: 3086 if num < Nx: # downsampling 3087 if real_input: 3088 sl[axis] = slice(N//2, N//2 + 1) 3089 Y[tuple(sl)] *= 2. 3090 else: 3091 # select the component of Y at frequency +N/2, 3092 # add the component of X at -N/2 3093 sl[axis] = slice(-N//2, -N//2 + 1) 3094 Y[tuple(sl)] += X[tuple(sl)] 3095 elif Nx < num: # upsampling 3096 # select the component at frequency +N/2 and halve it 3097 sl[axis] = slice(N//2, N//2 + 1) 3098 Y[tuple(sl)] *= 0.5 3099 if not real_input: 3100 temp = Y[tuple(sl)] 3101 # set the component at -N/2 equal to the component at +N/2 3102 sl[axis] = slice(num-N//2, num-N//2 + 1) 3103 Y[tuple(sl)] = temp 3104 3105 # Inverse transform 3106 if real_input: 3107 y = sp_fft.irfft(Y, num, axis=axis) 3108 else: 3109 y = sp_fft.ifft(Y, axis=axis, overwrite_x=True) 3110 3111 y *= (float(num) / float(Nx)) 3112 3113 if t is None: 3114 return y 3115 else: 3116 new_t = np.arange(0, num) * (t[1] - t[0]) * Nx / float(num) + t[0] 3117 return y, new_t 3118 3119 3120def resample_poly(x, up, down, axis=0, window=('kaiser', 5.0), 3121 padtype='constant', cval=None): 3122 """ 3123 Resample `x` along the given axis using polyphase filtering. 3124 3125 The signal `x` is upsampled by the factor `up`, a zero-phase low-pass 3126 FIR filter is applied, and then it is downsampled by the factor `down`. 3127 The resulting sample rate is ``up / down`` times the original sample 3128 rate. By default, values beyond the boundary of the signal are assumed 3129 to be zero during the filtering step. 3130 3131 Parameters 3132 ---------- 3133 x : array_like 3134 The data to be resampled. 3135 up : int 3136 The upsampling factor. 3137 down : int 3138 The downsampling factor. 3139 axis : int, optional 3140 The axis of `x` that is resampled. Default is 0. 3141 window : string, tuple, or array_like, optional 3142 Desired window to use to design the low-pass filter, or the FIR filter 3143 coefficients to employ. See below for details. 3144 padtype : string, optional 3145 `constant`, `line`, `mean`, `median`, `maximum`, `minimum` or any of 3146 the other signal extension modes supported by `scipy.signal.upfirdn`. 3147 Changes assumptions on values beyond the boundary. If `constant`, 3148 assumed to be `cval` (default zero). If `line` assumed to continue a 3149 linear trend defined by the first and last points. `mean`, `median`, 3150 `maximum` and `minimum` work as in `np.pad` and assume that the values 3151 beyond the boundary are the mean, median, maximum or minimum 3152 respectively of the array along the axis. 3153 3154 .. versionadded:: 1.4.0 3155 cval : float, optional 3156 Value to use if `padtype='constant'`. Default is zero. 3157 3158 .. versionadded:: 1.4.0 3159 3160 Returns 3161 ------- 3162 resampled_x : array 3163 The resampled array. 3164 3165 See Also 3166 -------- 3167 decimate : Downsample the signal after applying an FIR or IIR filter. 3168 resample : Resample up or down using the FFT method. 3169 3170 Notes 3171 ----- 3172 This polyphase method will likely be faster than the Fourier method 3173 in `scipy.signal.resample` when the number of samples is large and 3174 prime, or when the number of samples is large and `up` and `down` 3175 share a large greatest common denominator. The length of the FIR 3176 filter used will depend on ``max(up, down) // gcd(up, down)``, and 3177 the number of operations during polyphase filtering will depend on 3178 the filter length and `down` (see `scipy.signal.upfirdn` for details). 3179 3180 The argument `window` specifies the FIR low-pass filter design. 3181 3182 If `window` is an array_like it is assumed to be the FIR filter 3183 coefficients. Note that the FIR filter is applied after the upsampling 3184 step, so it should be designed to operate on a signal at a sampling 3185 frequency higher than the original by a factor of `up//gcd(up, down)`. 3186 This function's output will be centered with respect to this array, so it 3187 is best to pass a symmetric filter with an odd number of samples if, as 3188 is usually the case, a zero-phase filter is desired. 3189 3190 For any other type of `window`, the functions `scipy.signal.get_window` 3191 and `scipy.signal.firwin` are called to generate the appropriate filter 3192 coefficients. 3193 3194 The first sample of the returned vector is the same as the first 3195 sample of the input vector. The spacing between samples is changed 3196 from ``dx`` to ``dx * down / float(up)``. 3197 3198 Examples 3199 -------- 3200 By default, the end of the resampled data rises to meet the first 3201 sample of the next cycle for the FFT method, and gets closer to zero 3202 for the polyphase method: 3203 3204 >>> from scipy import signal 3205 3206 >>> x = np.linspace(0, 10, 20, endpoint=False) 3207 >>> y = np.cos(-x**2/6.0) 3208 >>> f_fft = signal.resample(y, 100) 3209 >>> f_poly = signal.resample_poly(y, 100, 20) 3210 >>> xnew = np.linspace(0, 10, 100, endpoint=False) 3211 3212 >>> import matplotlib.pyplot as plt 3213 >>> plt.plot(xnew, f_fft, 'b.-', xnew, f_poly, 'r.-') 3214 >>> plt.plot(x, y, 'ko-') 3215 >>> plt.plot(10, y[0], 'bo', 10, 0., 'ro') # boundaries 3216 >>> plt.legend(['resample', 'resamp_poly', 'data'], loc='best') 3217 >>> plt.show() 3218 3219 This default behaviour can be changed by using the padtype option: 3220 3221 >>> import numpy as np 3222 >>> from scipy import signal 3223 3224 >>> N = 5 3225 >>> x = np.linspace(0, 1, N, endpoint=False) 3226 >>> y = 2 + x**2 - 1.7*np.sin(x) + .2*np.cos(11*x) 3227 >>> y2 = 1 + x**3 + 0.1*np.sin(x) + .1*np.cos(11*x) 3228 >>> Y = np.stack([y, y2], axis=-1) 3229 >>> up = 4 3230 >>> xr = np.linspace(0, 1, N*up, endpoint=False) 3231 3232 >>> y2 = signal.resample_poly(Y, up, 1, padtype='constant') 3233 >>> y3 = signal.resample_poly(Y, up, 1, padtype='mean') 3234 >>> y4 = signal.resample_poly(Y, up, 1, padtype='line') 3235 3236 >>> import matplotlib.pyplot as plt 3237 >>> for i in [0,1]: 3238 ... plt.figure() 3239 ... plt.plot(xr, y4[:,i], 'g.', label='line') 3240 ... plt.plot(xr, y3[:,i], 'y.', label='mean') 3241 ... plt.plot(xr, y2[:,i], 'r.', label='constant') 3242 ... plt.plot(x, Y[:,i], 'k-') 3243 ... plt.legend() 3244 >>> plt.show() 3245 3246 """ 3247 x = np.asarray(x) 3248 if up != int(up): 3249 raise ValueError("up must be an integer") 3250 if down != int(down): 3251 raise ValueError("down must be an integer") 3252 up = int(up) 3253 down = int(down) 3254 if up < 1 or down < 1: 3255 raise ValueError('up and down must be >= 1') 3256 if cval is not None and padtype != 'constant': 3257 raise ValueError('cval has no effect when padtype is ', padtype) 3258 3259 # Determine our up and down factors 3260 # Use a rational approximation to save computation time on really long 3261 # signals 3262 g_ = math.gcd(up, down) 3263 up //= g_ 3264 down //= g_ 3265 if up == down == 1: 3266 return x.copy() 3267 n_in = x.shape[axis] 3268 n_out = n_in * up 3269 n_out = n_out // down + bool(n_out % down) 3270 3271 if isinstance(window, (list, np.ndarray)): 3272 window = np.array(window) # use array to force a copy (we modify it) 3273 if window.ndim > 1: 3274 raise ValueError('window must be 1-D') 3275 half_len = (window.size - 1) // 2 3276 h = window 3277 else: 3278 # Design a linear-phase low-pass FIR filter 3279 max_rate = max(up, down) 3280 f_c = 1. / max_rate # cutoff of FIR filter (rel. to Nyquist) 3281 half_len = 10 * max_rate # reasonable cutoff for sinc-like function 3282 h = firwin(2 * half_len + 1, f_c, window=window) 3283 h *= up 3284 3285 # Zero-pad our filter to put the output samples at the center 3286 n_pre_pad = (down - half_len % down) 3287 n_post_pad = 0 3288 n_pre_remove = (half_len + n_pre_pad) // down 3289 # We should rarely need to do this given our filter lengths... 3290 while _output_len(len(h) + n_pre_pad + n_post_pad, n_in, 3291 up, down) < n_out + n_pre_remove: 3292 n_post_pad += 1 3293 h = np.concatenate((np.zeros(n_pre_pad, dtype=h.dtype), h, 3294 np.zeros(n_post_pad, dtype=h.dtype))) 3295 n_pre_remove_end = n_pre_remove + n_out 3296 3297 # Remove background depending on the padtype option 3298 funcs = {'mean': np.mean, 'median': np.median, 3299 'minimum': np.amin, 'maximum': np.amax} 3300 upfirdn_kwargs = {'mode': 'constant', 'cval': 0} 3301 if padtype in funcs: 3302 background_values = funcs[padtype](x, axis=axis, keepdims=True) 3303 elif padtype in _upfirdn_modes: 3304 upfirdn_kwargs = {'mode': padtype} 3305 if padtype == 'constant': 3306 if cval is None: 3307 cval = 0 3308 upfirdn_kwargs['cval'] = cval 3309 else: 3310 raise ValueError( 3311 'padtype must be one of: maximum, mean, median, minimum, ' + 3312 ', '.join(_upfirdn_modes)) 3313 3314 if padtype in funcs: 3315 x = x - background_values 3316 3317 # filter then remove excess 3318 y = upfirdn(h, x, up, down, axis=axis, **upfirdn_kwargs) 3319 keep = [slice(None), ]*x.ndim 3320 keep[axis] = slice(n_pre_remove, n_pre_remove_end) 3321 y_keep = y[tuple(keep)] 3322 3323 # Add background back 3324 if padtype in funcs: 3325 y_keep += background_values 3326 3327 return y_keep 3328 3329 3330def vectorstrength(events, period): 3331 ''' 3332 Determine the vector strength of the events corresponding to the given 3333 period. 3334 3335 The vector strength is a measure of phase synchrony, how well the 3336 timing of the events is synchronized to a single period of a periodic 3337 signal. 3338 3339 If multiple periods are used, calculate the vector strength of each. 3340 This is called the "resonating vector strength". 3341 3342 Parameters 3343 ---------- 3344 events : 1D array_like 3345 An array of time points containing the timing of the events. 3346 period : float or array_like 3347 The period of the signal that the events should synchronize to. 3348 The period is in the same units as `events`. It can also be an array 3349 of periods, in which case the outputs are arrays of the same length. 3350 3351 Returns 3352 ------- 3353 strength : float or 1D array 3354 The strength of the synchronization. 1.0 is perfect synchronization 3355 and 0.0 is no synchronization. If `period` is an array, this is also 3356 an array with each element containing the vector strength at the 3357 corresponding period. 3358 phase : float or array 3359 The phase that the events are most strongly synchronized to in radians. 3360 If `period` is an array, this is also an array with each element 3361 containing the phase for the corresponding period. 3362 3363 References 3364 ---------- 3365 van Hemmen, JL, Longtin, A, and Vollmayr, AN. Testing resonating vector 3366 strength: Auditory system, electric fish, and noise. 3367 Chaos 21, 047508 (2011); 3368 :doi:`10.1063/1.3670512`. 3369 van Hemmen, JL. Vector strength after Goldberg, Brown, and von Mises: 3370 biological and mathematical perspectives. Biol Cybern. 3371 2013 Aug;107(4):385-96. :doi:`10.1007/s00422-013-0561-7`. 3372 van Hemmen, JL and Vollmayr, AN. Resonating vector strength: what happens 3373 when we vary the "probing" frequency while keeping the spike times 3374 fixed. Biol Cybern. 2013 Aug;107(4):491-94. 3375 :doi:`10.1007/s00422-013-0560-8`. 3376 ''' 3377 events = np.asarray(events) 3378 period = np.asarray(period) 3379 if events.ndim > 1: 3380 raise ValueError('events cannot have dimensions more than 1') 3381 if period.ndim > 1: 3382 raise ValueError('period cannot have dimensions more than 1') 3383 3384 # we need to know later if period was originally a scalar 3385 scalarperiod = not period.ndim 3386 3387 events = np.atleast_2d(events) 3388 period = np.atleast_2d(period) 3389 if (period <= 0).any(): 3390 raise ValueError('periods must be positive') 3391 3392 # this converts the times to vectors 3393 vectors = np.exp(np.dot(2j*np.pi/period.T, events)) 3394 3395 # the vector strength is just the magnitude of the mean of the vectors 3396 # the vector phase is the angle of the mean of the vectors 3397 vectormean = np.mean(vectors, axis=1) 3398 strength = abs(vectormean) 3399 phase = np.angle(vectormean) 3400 3401 # if the original period was a scalar, return scalars 3402 if scalarperiod: 3403 strength = strength[0] 3404 phase = phase[0] 3405 return strength, phase 3406 3407 3408def detrend(data, axis=-1, type='linear', bp=0, overwrite_data=False): 3409 """ 3410 Remove linear trend along axis from data. 3411 3412 Parameters 3413 ---------- 3414 data : array_like 3415 The input data. 3416 axis : int, optional 3417 The axis along which to detrend the data. By default this is the 3418 last axis (-1). 3419 type : {'linear', 'constant'}, optional 3420 The type of detrending. If ``type == 'linear'`` (default), 3421 the result of a linear least-squares fit to `data` is subtracted 3422 from `data`. 3423 If ``type == 'constant'``, only the mean of `data` is subtracted. 3424 bp : array_like of ints, optional 3425 A sequence of break points. If given, an individual linear fit is 3426 performed for each part of `data` between two break points. 3427 Break points are specified as indices into `data`. This parameter 3428 only has an effect when ``type == 'linear'``. 3429 overwrite_data : bool, optional 3430 If True, perform in place detrending and avoid a copy. Default is False 3431 3432 Returns 3433 ------- 3434 ret : ndarray 3435 The detrended input data. 3436 3437 Examples 3438 -------- 3439 >>> from scipy import signal 3440 >>> from numpy.random import default_rng 3441 >>> rng = default_rng() 3442 >>> npoints = 1000 3443 >>> noise = rng.standard_normal(npoints) 3444 >>> x = 3 + 2*np.linspace(0, 1, npoints) + noise 3445 >>> (signal.detrend(x) - noise).max() 3446 0.06 # random 3447 3448 """ 3449 if type not in ['linear', 'l', 'constant', 'c']: 3450 raise ValueError("Trend type must be 'linear' or 'constant'.") 3451 data = np.asarray(data) 3452 dtype = data.dtype.char 3453 if dtype not in 'dfDF': 3454 dtype = 'd' 3455 if type in ['constant', 'c']: 3456 ret = data - np.mean(data, axis, keepdims=True) 3457 return ret 3458 else: 3459 dshape = data.shape 3460 N = dshape[axis] 3461 bp = np.sort(np.unique(np.r_[0, bp, N])) 3462 if np.any(bp > N): 3463 raise ValueError("Breakpoints must be less than length " 3464 "of data along given axis.") 3465 Nreg = len(bp) - 1 3466 # Restructure data so that axis is along first dimension and 3467 # all other dimensions are collapsed into second dimension 3468 rnk = len(dshape) 3469 if axis < 0: 3470 axis = axis + rnk 3471 newdims = np.r_[axis, 0:axis, axis + 1:rnk] 3472 newdata = np.reshape(np.transpose(data, tuple(newdims)), 3473 (N, _prod(dshape) // N)) 3474 if not overwrite_data: 3475 newdata = newdata.copy() # make sure we have a copy 3476 if newdata.dtype.char not in 'dfDF': 3477 newdata = newdata.astype(dtype) 3478 # Find leastsq fit and remove it for each piece 3479 for m in range(Nreg): 3480 Npts = bp[m + 1] - bp[m] 3481 A = np.ones((Npts, 2), dtype) 3482 A[:, 0] = np.cast[dtype](np.arange(1, Npts + 1) * 1.0 / Npts) 3483 sl = slice(bp[m], bp[m + 1]) 3484 coef, resids, rank, s = linalg.lstsq(A, newdata[sl]) 3485 newdata[sl] = newdata[sl] - np.dot(A, coef) 3486 # Put data back in original shape. 3487 tdshape = np.take(dshape, newdims, 0) 3488 ret = np.reshape(newdata, tuple(tdshape)) 3489 vals = list(range(1, rnk)) 3490 olddims = vals[:axis] + [0] + vals[axis:] 3491 ret = np.transpose(ret, tuple(olddims)) 3492 return ret 3493 3494 3495def lfilter_zi(b, a): 3496 """ 3497 Construct initial conditions for lfilter for step response steady-state. 3498 3499 Compute an initial state `zi` for the `lfilter` function that corresponds 3500 to the steady state of the step response. 3501 3502 A typical use of this function is to set the initial state so that the 3503 output of the filter starts at the same value as the first element of 3504 the signal to be filtered. 3505 3506 Parameters 3507 ---------- 3508 b, a : array_like (1-D) 3509 The IIR filter coefficients. See `lfilter` for more 3510 information. 3511 3512 Returns 3513 ------- 3514 zi : 1-D ndarray 3515 The initial state for the filter. 3516 3517 See Also 3518 -------- 3519 lfilter, lfiltic, filtfilt 3520 3521 Notes 3522 ----- 3523 A linear filter with order m has a state space representation (A, B, C, D), 3524 for which the output y of the filter can be expressed as:: 3525 3526 z(n+1) = A*z(n) + B*x(n) 3527 y(n) = C*z(n) + D*x(n) 3528 3529 where z(n) is a vector of length m, A has shape (m, m), B has shape 3530 (m, 1), C has shape (1, m) and D has shape (1, 1) (assuming x(n) is 3531 a scalar). lfilter_zi solves:: 3532 3533 zi = A*zi + B 3534 3535 In other words, it finds the initial condition for which the response 3536 to an input of all ones is a constant. 3537 3538 Given the filter coefficients `a` and `b`, the state space matrices 3539 for the transposed direct form II implementation of the linear filter, 3540 which is the implementation used by scipy.signal.lfilter, are:: 3541 3542 A = scipy.linalg.companion(a).T 3543 B = b[1:] - a[1:]*b[0] 3544 3545 assuming `a[0]` is 1.0; if `a[0]` is not 1, `a` and `b` are first 3546 divided by a[0]. 3547 3548 Examples 3549 -------- 3550 The following code creates a lowpass Butterworth filter. Then it 3551 applies that filter to an array whose values are all 1.0; the 3552 output is also all 1.0, as expected for a lowpass filter. If the 3553 `zi` argument of `lfilter` had not been given, the output would have 3554 shown the transient signal. 3555 3556 >>> from numpy import array, ones 3557 >>> from scipy.signal import lfilter, lfilter_zi, butter 3558 >>> b, a = butter(5, 0.25) 3559 >>> zi = lfilter_zi(b, a) 3560 >>> y, zo = lfilter(b, a, ones(10), zi=zi) 3561 >>> y 3562 array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]) 3563 3564 Another example: 3565 3566 >>> x = array([0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0]) 3567 >>> y, zf = lfilter(b, a, x, zi=zi*x[0]) 3568 >>> y 3569 array([ 0.5 , 0.5 , 0.5 , 0.49836039, 0.48610528, 3570 0.44399389, 0.35505241]) 3571 3572 Note that the `zi` argument to `lfilter` was computed using 3573 `lfilter_zi` and scaled by `x[0]`. Then the output `y` has no 3574 transient until the input drops from 0.5 to 0.0. 3575 3576 """ 3577 3578 # FIXME: Can this function be replaced with an appropriate 3579 # use of lfiltic? For example, when b,a = butter(N,Wn), 3580 # lfiltic(b, a, y=numpy.ones_like(a), x=numpy.ones_like(b)). 3581 # 3582 3583 # We could use scipy.signal.normalize, but it uses warnings in 3584 # cases where a ValueError is more appropriate, and it allows 3585 # b to be 2D. 3586 b = np.atleast_1d(b) 3587 if b.ndim != 1: 3588 raise ValueError("Numerator b must be 1-D.") 3589 a = np.atleast_1d(a) 3590 if a.ndim != 1: 3591 raise ValueError("Denominator a must be 1-D.") 3592 3593 while len(a) > 1 and a[0] == 0.0: 3594 a = a[1:] 3595 if a.size < 1: 3596 raise ValueError("There must be at least one nonzero `a` coefficient.") 3597 3598 if a[0] != 1.0: 3599 # Normalize the coefficients so a[0] == 1. 3600 b = b / a[0] 3601 a = a / a[0] 3602 3603 n = max(len(a), len(b)) 3604 3605 # Pad a or b with zeros so they are the same length. 3606 if len(a) < n: 3607 a = np.r_[a, np.zeros(n - len(a))] 3608 elif len(b) < n: 3609 b = np.r_[b, np.zeros(n - len(b))] 3610 3611 IminusA = np.eye(n - 1, dtype=np.result_type(a, b)) - linalg.companion(a).T 3612 B = b[1:] - a[1:] * b[0] 3613 # Solve zi = A*zi + B 3614 zi = np.linalg.solve(IminusA, B) 3615 3616 # For future reference: we could also use the following 3617 # explicit formulas to solve the linear system: 3618 # 3619 # zi = np.zeros(n - 1) 3620 # zi[0] = B.sum() / IminusA[:,0].sum() 3621 # asum = 1.0 3622 # csum = 0.0 3623 # for k in range(1,n-1): 3624 # asum += a[k] 3625 # csum += b[k] - a[k]*b[0] 3626 # zi[k] = asum*zi[0] - csum 3627 3628 return zi 3629 3630 3631def sosfilt_zi(sos): 3632 """ 3633 Construct initial conditions for sosfilt for step response steady-state. 3634 3635 Compute an initial state `zi` for the `sosfilt` function that corresponds 3636 to the steady state of the step response. 3637 3638 A typical use of this function is to set the initial state so that the 3639 output of the filter starts at the same value as the first element of 3640 the signal to be filtered. 3641 3642 Parameters 3643 ---------- 3644 sos : array_like 3645 Array of second-order filter coefficients, must have shape 3646 ``(n_sections, 6)``. See `sosfilt` for the SOS filter format 3647 specification. 3648 3649 Returns 3650 ------- 3651 zi : ndarray 3652 Initial conditions suitable for use with ``sosfilt``, shape 3653 ``(n_sections, 2)``. 3654 3655 See Also 3656 -------- 3657 sosfilt, zpk2sos 3658 3659 Notes 3660 ----- 3661 .. versionadded:: 0.16.0 3662 3663 Examples 3664 -------- 3665 Filter a rectangular pulse that begins at time 0, with and without 3666 the use of the `zi` argument of `scipy.signal.sosfilt`. 3667 3668 >>> from scipy import signal 3669 >>> import matplotlib.pyplot as plt 3670 3671 >>> sos = signal.butter(9, 0.125, output='sos') 3672 >>> zi = signal.sosfilt_zi(sos) 3673 >>> x = (np.arange(250) < 100).astype(int) 3674 >>> f1 = signal.sosfilt(sos, x) 3675 >>> f2, zo = signal.sosfilt(sos, x, zi=zi) 3676 3677 >>> plt.plot(x, 'k--', label='x') 3678 >>> plt.plot(f1, 'b', alpha=0.5, linewidth=2, label='filtered') 3679 >>> plt.plot(f2, 'g', alpha=0.25, linewidth=4, label='filtered with zi') 3680 >>> plt.legend(loc='best') 3681 >>> plt.show() 3682 3683 """ 3684 sos = np.asarray(sos) 3685 if sos.ndim != 2 or sos.shape[1] != 6: 3686 raise ValueError('sos must be shape (n_sections, 6)') 3687 3688 if sos.dtype.kind in 'bui': 3689 sos = sos.astype(np.float64) 3690 3691 n_sections = sos.shape[0] 3692 zi = np.empty((n_sections, 2), dtype=sos.dtype) 3693 scale = 1.0 3694 for section in range(n_sections): 3695 b = sos[section, :3] 3696 a = sos[section, 3:] 3697 zi[section] = scale * lfilter_zi(b, a) 3698 # If H(z) = B(z)/A(z) is this section's transfer function, then 3699 # b.sum()/a.sum() is H(1), the gain at omega=0. That's the steady 3700 # state value of this section's step response. 3701 scale *= b.sum() / a.sum() 3702 3703 return zi 3704 3705 3706def _filtfilt_gust(b, a, x, axis=-1, irlen=None): 3707 """Forward-backward IIR filter that uses Gustafsson's method. 3708 3709 Apply the IIR filter defined by `(b,a)` to `x` twice, first forward 3710 then backward, using Gustafsson's initial conditions [1]_. 3711 3712 Let ``y_fb`` be the result of filtering first forward and then backward, 3713 and let ``y_bf`` be the result of filtering first backward then forward. 3714 Gustafsson's method is to compute initial conditions for the forward 3715 pass and the backward pass such that ``y_fb == y_bf``. 3716 3717 Parameters 3718 ---------- 3719 b : scalar or 1-D ndarray 3720 Numerator coefficients of the filter. 3721 a : scalar or 1-D ndarray 3722 Denominator coefficients of the filter. 3723 x : ndarray 3724 Data to be filtered. 3725 axis : int, optional 3726 Axis of `x` to be filtered. Default is -1. 3727 irlen : int or None, optional 3728 The length of the nonnegligible part of the impulse response. 3729 If `irlen` is None, or if the length of the signal is less than 3730 ``2 * irlen``, then no part of the impulse response is ignored. 3731 3732 Returns 3733 ------- 3734 y : ndarray 3735 The filtered data. 3736 x0 : ndarray 3737 Initial condition for the forward filter. 3738 x1 : ndarray 3739 Initial condition for the backward filter. 3740 3741 Notes 3742 ----- 3743 Typically the return values `x0` and `x1` are not needed by the 3744 caller. The intended use of these return values is in unit tests. 3745 3746 References 3747 ---------- 3748 .. [1] F. Gustaffson. Determining the initial states in forward-backward 3749 filtering. Transactions on Signal Processing, 46(4):988-992, 1996. 3750 3751 """ 3752 # In the comments, "Gustafsson's paper" and [1] refer to the 3753 # paper referenced in the docstring. 3754 3755 b = np.atleast_1d(b) 3756 a = np.atleast_1d(a) 3757 3758 order = max(len(b), len(a)) - 1 3759 if order == 0: 3760 # The filter is just scalar multiplication, with no state. 3761 scale = (b[0] / a[0])**2 3762 y = scale * x 3763 return y, np.array([]), np.array([]) 3764 3765 if axis != -1 or axis != x.ndim - 1: 3766 # Move the axis containing the data to the end. 3767 x = np.swapaxes(x, axis, x.ndim - 1) 3768 3769 # n is the number of samples in the data to be filtered. 3770 n = x.shape[-1] 3771 3772 if irlen is None or n <= 2*irlen: 3773 m = n 3774 else: 3775 m = irlen 3776 3777 # Create Obs, the observability matrix (called O in the paper). 3778 # This matrix can be interpreted as the operator that propagates 3779 # an arbitrary initial state to the output, assuming the input is 3780 # zero. 3781 # In Gustafsson's paper, the forward and backward filters are not 3782 # necessarily the same, so he has both O_f and O_b. We use the same 3783 # filter in both directions, so we only need O. The same comment 3784 # applies to S below. 3785 Obs = np.zeros((m, order)) 3786 zi = np.zeros(order) 3787 zi[0] = 1 3788 Obs[:, 0] = lfilter(b, a, np.zeros(m), zi=zi)[0] 3789 for k in range(1, order): 3790 Obs[k:, k] = Obs[:-k, 0] 3791 3792 # Obsr is O^R (Gustafsson's notation for row-reversed O) 3793 Obsr = Obs[::-1] 3794 3795 # Create S. S is the matrix that applies the filter to the reversed 3796 # propagated initial conditions. That is, 3797 # out = S.dot(zi) 3798 # is the same as 3799 # tmp, _ = lfilter(b, a, zeros(), zi=zi) # Propagate ICs. 3800 # out = lfilter(b, a, tmp[::-1]) # Reverse and filter. 3801 3802 # Equations (5) & (6) of [1] 3803 S = lfilter(b, a, Obs[::-1], axis=0) 3804 3805 # Sr is S^R (row-reversed S) 3806 Sr = S[::-1] 3807 3808 # M is [(S^R - O), (O^R - S)] 3809 if m == n: 3810 M = np.hstack((Sr - Obs, Obsr - S)) 3811 else: 3812 # Matrix described in section IV of [1]. 3813 M = np.zeros((2*m, 2*order)) 3814 M[:m, :order] = Sr - Obs 3815 M[m:, order:] = Obsr - S 3816 3817 # Naive forward-backward and backward-forward filters. 3818 # These have large transients because the filters use zero initial 3819 # conditions. 3820 y_f = lfilter(b, a, x) 3821 y_fb = lfilter(b, a, y_f[..., ::-1])[..., ::-1] 3822 3823 y_b = lfilter(b, a, x[..., ::-1])[..., ::-1] 3824 y_bf = lfilter(b, a, y_b) 3825 3826 delta_y_bf_fb = y_bf - y_fb 3827 if m == n: 3828 delta = delta_y_bf_fb 3829 else: 3830 start_m = delta_y_bf_fb[..., :m] 3831 end_m = delta_y_bf_fb[..., -m:] 3832 delta = np.concatenate((start_m, end_m), axis=-1) 3833 3834 # ic_opt holds the "optimal" initial conditions. 3835 # The following code computes the result shown in the formula 3836 # of the paper between equations (6) and (7). 3837 if delta.ndim == 1: 3838 ic_opt = linalg.lstsq(M, delta)[0] 3839 else: 3840 # Reshape delta so it can be used as an array of multiple 3841 # right-hand-sides in linalg.lstsq. 3842 delta2d = delta.reshape(-1, delta.shape[-1]).T 3843 ic_opt0 = linalg.lstsq(M, delta2d)[0].T 3844 ic_opt = ic_opt0.reshape(delta.shape[:-1] + (M.shape[-1],)) 3845 3846 # Now compute the filtered signal using equation (7) of [1]. 3847 # First, form [S^R, O^R] and call it W. 3848 if m == n: 3849 W = np.hstack((Sr, Obsr)) 3850 else: 3851 W = np.zeros((2*m, 2*order)) 3852 W[:m, :order] = Sr 3853 W[m:, order:] = Obsr 3854 3855 # Equation (7) of [1] says 3856 # Y_fb^opt = Y_fb^0 + W * [x_0^opt; x_{N-1}^opt] 3857 # `wic` is (almost) the product on the right. 3858 # W has shape (m, 2*order), and ic_opt has shape (..., 2*order), 3859 # so we can't use W.dot(ic_opt). Instead, we dot ic_opt with W.T, 3860 # so wic has shape (..., m). 3861 wic = ic_opt.dot(W.T) 3862 3863 # `wic` is "almost" the product of W and the optimal ICs in equation 3864 # (7)--if we're using a truncated impulse response (m < n), `wic` 3865 # contains only the adjustments required for the ends of the signal. 3866 # Here we form y_opt, taking this into account if necessary. 3867 y_opt = y_fb 3868 if m == n: 3869 y_opt += wic 3870 else: 3871 y_opt[..., :m] += wic[..., :m] 3872 y_opt[..., -m:] += wic[..., -m:] 3873 3874 x0 = ic_opt[..., :order] 3875 x1 = ic_opt[..., -order:] 3876 if axis != -1 or axis != x.ndim - 1: 3877 # Restore the data axis to its original position. 3878 x0 = np.swapaxes(x0, axis, x.ndim - 1) 3879 x1 = np.swapaxes(x1, axis, x.ndim - 1) 3880 y_opt = np.swapaxes(y_opt, axis, x.ndim - 1) 3881 3882 return y_opt, x0, x1 3883 3884 3885def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad', 3886 irlen=None): 3887 """ 3888 Apply a digital filter forward and backward to a signal. 3889 3890 This function applies a linear digital filter twice, once forward and 3891 once backwards. The combined filter has zero phase and a filter order 3892 twice that of the original. 3893 3894 The function provides options for handling the edges of the signal. 3895 3896 The function `sosfiltfilt` (and filter design using ``output='sos'``) 3897 should be preferred over `filtfilt` for most filtering tasks, as 3898 second-order sections have fewer numerical problems. 3899 3900 Parameters 3901 ---------- 3902 b : (N,) array_like 3903 The numerator coefficient vector of the filter. 3904 a : (N,) array_like 3905 The denominator coefficient vector of the filter. If ``a[0]`` 3906 is not 1, then both `a` and `b` are normalized by ``a[0]``. 3907 x : array_like 3908 The array of data to be filtered. 3909 axis : int, optional 3910 The axis of `x` to which the filter is applied. 3911 Default is -1. 3912 padtype : str or None, optional 3913 Must be 'odd', 'even', 'constant', or None. This determines the 3914 type of extension to use for the padded signal to which the filter 3915 is applied. If `padtype` is None, no padding is used. The default 3916 is 'odd'. 3917 padlen : int or None, optional 3918 The number of elements by which to extend `x` at both ends of 3919 `axis` before applying the filter. This value must be less than 3920 ``x.shape[axis] - 1``. ``padlen=0`` implies no padding. 3921 The default value is ``3 * max(len(a), len(b))``. 3922 method : str, optional 3923 Determines the method for handling the edges of the signal, either 3924 "pad" or "gust". When `method` is "pad", the signal is padded; the 3925 type of padding is determined by `padtype` and `padlen`, and `irlen` 3926 is ignored. When `method` is "gust", Gustafsson's method is used, 3927 and `padtype` and `padlen` are ignored. 3928 irlen : int or None, optional 3929 When `method` is "gust", `irlen` specifies the length of the 3930 impulse response of the filter. If `irlen` is None, no part 3931 of the impulse response is ignored. For a long signal, specifying 3932 `irlen` can significantly improve the performance of the filter. 3933 3934 Returns 3935 ------- 3936 y : ndarray 3937 The filtered output with the same shape as `x`. 3938 3939 See Also 3940 -------- 3941 sosfiltfilt, lfilter_zi, lfilter, lfiltic, savgol_filter, sosfilt 3942 3943 Notes 3944 ----- 3945 When `method` is "pad", the function pads the data along the given axis 3946 in one of three ways: odd, even or constant. The odd and even extensions 3947 have the corresponding symmetry about the end point of the data. The 3948 constant extension extends the data with the values at the end points. On 3949 both the forward and backward passes, the initial condition of the 3950 filter is found by using `lfilter_zi` and scaling it by the end point of 3951 the extended data. 3952 3953 When `method` is "gust", Gustafsson's method [1]_ is used. Initial 3954 conditions are chosen for the forward and backward passes so that the 3955 forward-backward filter gives the same result as the backward-forward 3956 filter. 3957 3958 The option to use Gustaffson's method was added in scipy version 0.16.0. 3959 3960 References 3961 ---------- 3962 .. [1] F. Gustaffson, "Determining the initial states in forward-backward 3963 filtering", Transactions on Signal Processing, Vol. 46, pp. 988-992, 3964 1996. 3965 3966 Examples 3967 -------- 3968 The examples will use several functions from `scipy.signal`. 3969 3970 >>> from scipy import signal 3971 >>> import matplotlib.pyplot as plt 3972 3973 First we create a one second signal that is the sum of two pure sine 3974 waves, with frequencies 5 Hz and 250 Hz, sampled at 2000 Hz. 3975 3976 >>> t = np.linspace(0, 1.0, 2001) 3977 >>> xlow = np.sin(2 * np.pi * 5 * t) 3978 >>> xhigh = np.sin(2 * np.pi * 250 * t) 3979 >>> x = xlow + xhigh 3980 3981 Now create a lowpass Butterworth filter with a cutoff of 0.125 times 3982 the Nyquist frequency, or 125 Hz, and apply it to ``x`` with `filtfilt`. 3983 The result should be approximately ``xlow``, with no phase shift. 3984 3985 >>> b, a = signal.butter(8, 0.125) 3986 >>> y = signal.filtfilt(b, a, x, padlen=150) 3987 >>> np.abs(y - xlow).max() 3988 9.1086182074789912e-06 3989 3990 We get a fairly clean result for this artificial example because 3991 the odd extension is exact, and with the moderately long padding, 3992 the filter's transients have dissipated by the time the actual data 3993 is reached. In general, transient effects at the edges are 3994 unavoidable. 3995 3996 The following example demonstrates the option ``method="gust"``. 3997 3998 First, create a filter. 3999 4000 >>> b, a = signal.ellip(4, 0.01, 120, 0.125) # Filter to be applied. 4001 4002 `sig` is a random input signal to be filtered. 4003 4004 >>> rng = np.random.default_rng() 4005 >>> n = 60 4006 >>> sig = rng.standard_normal(n)**3 + 3*rng.standard_normal(n).cumsum() 4007 4008 Apply `filtfilt` to `sig`, once using the Gustafsson method, and 4009 once using padding, and plot the results for comparison. 4010 4011 >>> fgust = signal.filtfilt(b, a, sig, method="gust") 4012 >>> fpad = signal.filtfilt(b, a, sig, padlen=50) 4013 >>> plt.plot(sig, 'k-', label='input') 4014 >>> plt.plot(fgust, 'b-', linewidth=4, label='gust') 4015 >>> plt.plot(fpad, 'c-', linewidth=1.5, label='pad') 4016 >>> plt.legend(loc='best') 4017 >>> plt.show() 4018 4019 The `irlen` argument can be used to improve the performance 4020 of Gustafsson's method. 4021 4022 Estimate the impulse response length of the filter. 4023 4024 >>> z, p, k = signal.tf2zpk(b, a) 4025 >>> eps = 1e-9 4026 >>> r = np.max(np.abs(p)) 4027 >>> approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r))) 4028 >>> approx_impulse_len 4029 137 4030 4031 Apply the filter to a longer signal, with and without the `irlen` 4032 argument. The difference between `y1` and `y2` is small. For long 4033 signals, using `irlen` gives a significant performance improvement. 4034 4035 >>> x = rng.standard_normal(5000) 4036 >>> y1 = signal.filtfilt(b, a, x, method='gust') 4037 >>> y2 = signal.filtfilt(b, a, x, method='gust', irlen=approx_impulse_len) 4038 >>> print(np.max(np.abs(y1 - y2))) 4039 1.80056858312e-10 4040 4041 """ 4042 b = np.atleast_1d(b) 4043 a = np.atleast_1d(a) 4044 x = np.asarray(x) 4045 4046 if method not in ["pad", "gust"]: 4047 raise ValueError("method must be 'pad' or 'gust'.") 4048 4049 if method == "gust": 4050 y, z1, z2 = _filtfilt_gust(b, a, x, axis=axis, irlen=irlen) 4051 return y 4052 4053 # method == "pad" 4054 edge, ext = _validate_pad(padtype, padlen, x, axis, 4055 ntaps=max(len(a), len(b))) 4056 4057 # Get the steady state of the filter's step response. 4058 zi = lfilter_zi(b, a) 4059 4060 # Reshape zi and create x0 so that zi*x0 broadcasts 4061 # to the correct value for the 'zi' keyword argument 4062 # to lfilter. 4063 zi_shape = [1] * x.ndim 4064 zi_shape[axis] = zi.size 4065 zi = np.reshape(zi, zi_shape) 4066 x0 = axis_slice(ext, stop=1, axis=axis) 4067 4068 # Forward filter. 4069 (y, zf) = lfilter(b, a, ext, axis=axis, zi=zi * x0) 4070 4071 # Backward filter. 4072 # Create y0 so zi*y0 broadcasts appropriately. 4073 y0 = axis_slice(y, start=-1, axis=axis) 4074 (y, zf) = lfilter(b, a, axis_reverse(y, axis=axis), axis=axis, zi=zi * y0) 4075 4076 # Reverse y. 4077 y = axis_reverse(y, axis=axis) 4078 4079 if edge > 0: 4080 # Slice the actual signal from the extended signal. 4081 y = axis_slice(y, start=edge, stop=-edge, axis=axis) 4082 4083 return y 4084 4085 4086def _validate_pad(padtype, padlen, x, axis, ntaps): 4087 """Helper to validate padding for filtfilt""" 4088 if padtype not in ['even', 'odd', 'constant', None]: 4089 raise ValueError(("Unknown value '%s' given to padtype. padtype " 4090 "must be 'even', 'odd', 'constant', or None.") % 4091 padtype) 4092 4093 if padtype is None: 4094 padlen = 0 4095 4096 if padlen is None: 4097 # Original padding; preserved for backwards compatibility. 4098 edge = ntaps * 3 4099 else: 4100 edge = padlen 4101 4102 # x's 'axis' dimension must be bigger than edge. 4103 if x.shape[axis] <= edge: 4104 raise ValueError("The length of the input vector x must be greater " 4105 "than padlen, which is %d." % edge) 4106 4107 if padtype is not None and edge > 0: 4108 # Make an extension of length `edge` at each 4109 # end of the input array. 4110 if padtype == 'even': 4111 ext = even_ext(x, edge, axis=axis) 4112 elif padtype == 'odd': 4113 ext = odd_ext(x, edge, axis=axis) 4114 else: 4115 ext = const_ext(x, edge, axis=axis) 4116 else: 4117 ext = x 4118 return edge, ext 4119 4120 4121def _validate_x(x): 4122 x = np.asarray(x) 4123 if x.ndim == 0: 4124 raise ValueError('x must be at least 1-D') 4125 return x 4126 4127 4128def sosfilt(sos, x, axis=-1, zi=None): 4129 """ 4130 Filter data along one dimension using cascaded second-order sections. 4131 4132 Filter a data sequence, `x`, using a digital IIR filter defined by 4133 `sos`. 4134 4135 Parameters 4136 ---------- 4137 sos : array_like 4138 Array of second-order filter coefficients, must have shape 4139 ``(n_sections, 6)``. Each row corresponds to a second-order 4140 section, with the first three columns providing the numerator 4141 coefficients and the last three providing the denominator 4142 coefficients. 4143 x : array_like 4144 An N-dimensional input array. 4145 axis : int, optional 4146 The axis of the input data array along which to apply the 4147 linear filter. The filter is applied to each subarray along 4148 this axis. Default is -1. 4149 zi : array_like, optional 4150 Initial conditions for the cascaded filter delays. It is a (at 4151 least 2D) vector of shape ``(n_sections, ..., 2, ...)``, where 4152 ``..., 2, ...`` denotes the shape of `x`, but with ``x.shape[axis]`` 4153 replaced by 2. If `zi` is None or is not given then initial rest 4154 (i.e. all zeros) is assumed. 4155 Note that these initial conditions are *not* the same as the initial 4156 conditions given by `lfiltic` or `lfilter_zi`. 4157 4158 Returns 4159 ------- 4160 y : ndarray 4161 The output of the digital filter. 4162 zf : ndarray, optional 4163 If `zi` is None, this is not returned, otherwise, `zf` holds the 4164 final filter delay values. 4165 4166 See Also 4167 -------- 4168 zpk2sos, sos2zpk, sosfilt_zi, sosfiltfilt, sosfreqz 4169 4170 Notes 4171 ----- 4172 The filter function is implemented as a series of second-order filters 4173 with direct-form II transposed structure. It is designed to minimize 4174 numerical precision errors for high-order filters. 4175 4176 .. versionadded:: 0.16.0 4177 4178 Examples 4179 -------- 4180 Plot a 13th-order filter's impulse response using both `lfilter` and 4181 `sosfilt`, showing the instability that results from trying to do a 4182 13th-order filter in a single stage (the numerical error pushes some poles 4183 outside of the unit circle): 4184 4185 >>> import matplotlib.pyplot as plt 4186 >>> from scipy import signal 4187 >>> b, a = signal.ellip(13, 0.009, 80, 0.05, output='ba') 4188 >>> sos = signal.ellip(13, 0.009, 80, 0.05, output='sos') 4189 >>> x = signal.unit_impulse(700) 4190 >>> y_tf = signal.lfilter(b, a, x) 4191 >>> y_sos = signal.sosfilt(sos, x) 4192 >>> plt.plot(y_tf, 'r', label='TF') 4193 >>> plt.plot(y_sos, 'k', label='SOS') 4194 >>> plt.legend(loc='best') 4195 >>> plt.show() 4196 4197 """ 4198 x = _validate_x(x) 4199 sos, n_sections = _validate_sos(sos) 4200 x_zi_shape = list(x.shape) 4201 x_zi_shape[axis] = 2 4202 x_zi_shape = tuple([n_sections] + x_zi_shape) 4203 inputs = [sos, x] 4204 if zi is not None: 4205 inputs.append(np.asarray(zi)) 4206 dtype = np.result_type(*inputs) 4207 if dtype.char not in 'fdgFDGO': 4208 raise NotImplementedError("input type '%s' not supported" % dtype) 4209 if zi is not None: 4210 zi = np.array(zi, dtype) # make a copy so that we can operate in place 4211 if zi.shape != x_zi_shape: 4212 raise ValueError('Invalid zi shape. With axis=%r, an input with ' 4213 'shape %r, and an sos array with %d sections, zi ' 4214 'must have shape %r, got %r.' % 4215 (axis, x.shape, n_sections, x_zi_shape, zi.shape)) 4216 return_zi = True 4217 else: 4218 zi = np.zeros(x_zi_shape, dtype=dtype) 4219 return_zi = False 4220 axis = axis % x.ndim # make positive 4221 x = np.moveaxis(x, axis, -1) 4222 zi = np.moveaxis(zi, [0, axis + 1], [-2, -1]) 4223 x_shape, zi_shape = x.shape, zi.shape 4224 x = np.reshape(x, (-1, x.shape[-1])) 4225 x = np.array(x, dtype, order='C') # make a copy, can modify in place 4226 zi = np.ascontiguousarray(np.reshape(zi, (-1, n_sections, 2))) 4227 sos = sos.astype(dtype, copy=False) 4228 _sosfilt(sos, x, zi) 4229 x.shape = x_shape 4230 x = np.moveaxis(x, -1, axis) 4231 if return_zi: 4232 zi.shape = zi_shape 4233 zi = np.moveaxis(zi, [-2, -1], [0, axis + 1]) 4234 out = (x, zi) 4235 else: 4236 out = x 4237 return out 4238 4239 4240def sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None): 4241 """ 4242 A forward-backward digital filter using cascaded second-order sections. 4243 4244 See `filtfilt` for more complete information about this method. 4245 4246 Parameters 4247 ---------- 4248 sos : array_like 4249 Array of second-order filter coefficients, must have shape 4250 ``(n_sections, 6)``. Each row corresponds to a second-order 4251 section, with the first three columns providing the numerator 4252 coefficients and the last three providing the denominator 4253 coefficients. 4254 x : array_like 4255 The array of data to be filtered. 4256 axis : int, optional 4257 The axis of `x` to which the filter is applied. 4258 Default is -1. 4259 padtype : str or None, optional 4260 Must be 'odd', 'even', 'constant', or None. This determines the 4261 type of extension to use for the padded signal to which the filter 4262 is applied. If `padtype` is None, no padding is used. The default 4263 is 'odd'. 4264 padlen : int or None, optional 4265 The number of elements by which to extend `x` at both ends of 4266 `axis` before applying the filter. This value must be less than 4267 ``x.shape[axis] - 1``. ``padlen=0`` implies no padding. 4268 The default value is:: 4269 4270 3 * (2 * len(sos) + 1 - min((sos[:, 2] == 0).sum(), 4271 (sos[:, 5] == 0).sum())) 4272 4273 The extra subtraction at the end attempts to compensate for poles 4274 and zeros at the origin (e.g. for odd-order filters) to yield 4275 equivalent estimates of `padlen` to those of `filtfilt` for 4276 second-order section filters built with `scipy.signal` functions. 4277 4278 Returns 4279 ------- 4280 y : ndarray 4281 The filtered output with the same shape as `x`. 4282 4283 See Also 4284 -------- 4285 filtfilt, sosfilt, sosfilt_zi, sosfreqz 4286 4287 Notes 4288 ----- 4289 .. versionadded:: 0.18.0 4290 4291 Examples 4292 -------- 4293 >>> from scipy.signal import sosfiltfilt, butter 4294 >>> import matplotlib.pyplot as plt 4295 >>> rng = np.random.default_rng() 4296 4297 Create an interesting signal to filter. 4298 4299 >>> n = 201 4300 >>> t = np.linspace(0, 1, n) 4301 >>> x = 1 + (t < 0.5) - 0.25*t**2 + 0.05*rng.standard_normal(n) 4302 4303 Create a lowpass Butterworth filter, and use it to filter `x`. 4304 4305 >>> sos = butter(4, 0.125, output='sos') 4306 >>> y = sosfiltfilt(sos, x) 4307 4308 For comparison, apply an 8th order filter using `sosfilt`. The filter 4309 is initialized using the mean of the first four values of `x`. 4310 4311 >>> from scipy.signal import sosfilt, sosfilt_zi 4312 >>> sos8 = butter(8, 0.125, output='sos') 4313 >>> zi = x[:4].mean() * sosfilt_zi(sos8) 4314 >>> y2, zo = sosfilt(sos8, x, zi=zi) 4315 4316 Plot the results. Note that the phase of `y` matches the input, while 4317 `y2` has a significant phase delay. 4318 4319 >>> plt.plot(t, x, alpha=0.5, label='x(t)') 4320 >>> plt.plot(t, y, label='y(t)') 4321 >>> plt.plot(t, y2, label='y2(t)') 4322 >>> plt.legend(framealpha=1, shadow=True) 4323 >>> plt.grid(alpha=0.25) 4324 >>> plt.xlabel('t') 4325 >>> plt.show() 4326 4327 """ 4328 sos, n_sections = _validate_sos(sos) 4329 x = _validate_x(x) 4330 4331 # `method` is "pad"... 4332 ntaps = 2 * n_sections + 1 4333 ntaps -= min((sos[:, 2] == 0).sum(), (sos[:, 5] == 0).sum()) 4334 edge, ext = _validate_pad(padtype, padlen, x, axis, 4335 ntaps=ntaps) 4336 4337 # These steps follow the same form as filtfilt with modifications 4338 zi = sosfilt_zi(sos) # shape (n_sections, 2) --> (n_sections, ..., 2, ...) 4339 zi_shape = [1] * x.ndim 4340 zi_shape[axis] = 2 4341 zi.shape = [n_sections] + zi_shape 4342 x_0 = axis_slice(ext, stop=1, axis=axis) 4343 (y, zf) = sosfilt(sos, ext, axis=axis, zi=zi * x_0) 4344 y_0 = axis_slice(y, start=-1, axis=axis) 4345 (y, zf) = sosfilt(sos, axis_reverse(y, axis=axis), axis=axis, zi=zi * y_0) 4346 y = axis_reverse(y, axis=axis) 4347 if edge > 0: 4348 y = axis_slice(y, start=edge, stop=-edge, axis=axis) 4349 return y 4350 4351 4352def decimate(x, q, n=None, ftype='iir', axis=-1, zero_phase=True): 4353 """ 4354 Downsample the signal after applying an anti-aliasing filter. 4355 4356 By default, an order 8 Chebyshev type I filter is used. A 30 point FIR 4357 filter with Hamming window is used if `ftype` is 'fir'. 4358 4359 Parameters 4360 ---------- 4361 x : array_like 4362 The signal to be downsampled, as an N-dimensional array. 4363 q : int 4364 The downsampling factor. When using IIR downsampling, it is recommended 4365 to call `decimate` multiple times for downsampling factors higher than 4366 13. 4367 n : int, optional 4368 The order of the filter (1 less than the length for 'fir'). Defaults to 4369 8 for 'iir' and 20 times the downsampling factor for 'fir'. 4370 ftype : str {'iir', 'fir'} or ``dlti`` instance, optional 4371 If 'iir' or 'fir', specifies the type of lowpass filter. If an instance 4372 of an `dlti` object, uses that object to filter before downsampling. 4373 axis : int, optional 4374 The axis along which to decimate. 4375 zero_phase : bool, optional 4376 Prevent phase shift by filtering with `filtfilt` instead of `lfilter` 4377 when using an IIR filter, and shifting the outputs back by the filter's 4378 group delay when using an FIR filter. The default value of ``True`` is 4379 recommended, since a phase shift is generally not desired. 4380 4381 .. versionadded:: 0.18.0 4382 4383 Returns 4384 ------- 4385 y : ndarray 4386 The down-sampled signal. 4387 4388 See Also 4389 -------- 4390 resample : Resample up or down using the FFT method. 4391 resample_poly : Resample using polyphase filtering and an FIR filter. 4392 4393 Notes 4394 ----- 4395 The ``zero_phase`` keyword was added in 0.18.0. 4396 The possibility to use instances of ``dlti`` as ``ftype`` was added in 4397 0.18.0. 4398 """ 4399 4400 x = np.asarray(x) 4401 q = operator.index(q) 4402 4403 if n is not None: 4404 n = operator.index(n) 4405 4406 if ftype == 'fir': 4407 if n is None: 4408 half_len = 10 * q # reasonable cutoff for our sinc-like function 4409 n = 2 * half_len 4410 b, a = firwin(n+1, 1. / q, window='hamming'), 1. 4411 elif ftype == 'iir': 4412 if n is None: 4413 n = 8 4414 system = dlti(*cheby1(n, 0.05, 0.8 / q)) 4415 b, a = system.num, system.den 4416 elif isinstance(ftype, dlti): 4417 system = ftype._as_tf() # Avoids copying if already in TF form 4418 b, a = system.num, system.den 4419 else: 4420 raise ValueError('invalid ftype') 4421 4422 result_type = x.dtype 4423 if result_type.kind in 'bui': 4424 result_type = np.float64 4425 b = np.asarray(b, dtype=result_type) 4426 a = np.asarray(a, dtype=result_type) 4427 4428 sl = [slice(None)] * x.ndim 4429 a = np.asarray(a) 4430 4431 if a.size == 1: # FIR case 4432 b = b / a 4433 if zero_phase: 4434 y = resample_poly(x, 1, q, axis=axis, window=b) 4435 else: 4436 # upfirdn is generally faster than lfilter by a factor equal to the 4437 # downsampling factor, since it only calculates the needed outputs 4438 n_out = x.shape[axis] // q + bool(x.shape[axis] % q) 4439 y = upfirdn(b, x, up=1, down=q, axis=axis) 4440 sl[axis] = slice(None, n_out, None) 4441 4442 else: # IIR case 4443 if zero_phase: 4444 y = filtfilt(b, a, x, axis=axis) 4445 else: 4446 y = lfilter(b, a, x, axis=axis) 4447 sl[axis] = slice(None, None, q) 4448 4449 return y[tuple(sl)] 4450