Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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', '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: 

43 raise ValueError("Acceptable mode flags are 'valid'," 

44 " 'same', or 'full'.") 

45 

46 

47def _bvalfromboundary(boundary): 

48 try: 

49 return _boundarydict[boundary] << 2 

50 except KeyError: 

51 raise ValueError("Acceptable boundary flags are 'fill', 'circular' " 

52 "(or 'wrap'), and 'symmetric' (or 'symm').") 

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 

140 Notes 

141 ----- 

142 The correlation z of two d-dimensional arrays x and y is defined as:: 

143 

144 z[...,k,...] = sum[..., i_l, ...] x[..., i_l,...] * conj(y[..., i_l - k,...]) 

145 

146 This way, if x and y are 1-D arrays and ``z = correlate(x, y, 'full')`` 

147 then 

148 

149 .. math:: 

150 

151 z[k] = (x * y)(k - N + 1) 

152 = \sum_{l=0}^{||x||-1}x_l y_{l-k+N-1}^{*} 

153 

154 for :math:`k = 0, 1, ..., ||x|| + ||y|| - 2` 

155 

156 where :math:`||x||` is the length of ``x``, :math:`N = \max(||x||,||y||)`, 

157 and :math:`y_m` is 0 when m is outside the range of y. 

158 

159 ``method='fft'`` only works for numerical arrays as it relies on 

160 `fftconvolve`. In certain cases (i.e., arrays of objects or when 

161 rounding integers can lose precision), ``method='direct'`` is always used. 

162 

163 When using "same" mode with even-length inputs, the outputs of `correlate` 

164 and `correlate2d` differ: There is a 1-index offset between them. 

165 

166 Examples 

167 -------- 

168 Implement a matched filter using cross-correlation, to recover a signal 

169 that has passed through a noisy channel. 

170 

171 >>> from scipy import signal 

172 >>> sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 128) 

173 >>> sig_noise = sig + np.random.randn(len(sig)) 

174 >>> corr = signal.correlate(sig_noise, np.ones(128), mode='same') / 128 

175 

176 >>> import matplotlib.pyplot as plt 

177 >>> clock = np.arange(64, len(sig), 128) 

178 >>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True) 

179 >>> ax_orig.plot(sig) 

180 >>> ax_orig.plot(clock, sig[clock], 'ro') 

181 >>> ax_orig.set_title('Original signal') 

182 >>> ax_noise.plot(sig_noise) 

183 >>> ax_noise.set_title('Signal with noise') 

184 >>> ax_corr.plot(corr) 

185 >>> ax_corr.plot(clock, corr[clock], 'ro') 

186 >>> ax_corr.axhline(0.5, ls=':') 

187 >>> ax_corr.set_title('Cross-correlated with rectangular pulse') 

188 >>> ax_orig.margins(0, 0.1) 

189 >>> fig.tight_layout() 

190 >>> fig.show() 

191 

192 """ 

193 in1 = np.asarray(in1) 

194 in2 = np.asarray(in2) 

195 

196 if in1.ndim == in2.ndim == 0: 

197 return in1 * in2.conj() 

198 elif in1.ndim != in2.ndim: 

199 raise ValueError("in1 and in2 should have the same dimensionality") 

200 

201 # Don't use _valfrommode, since correlate should not accept numeric modes 

202 try: 

203 val = _modedict[mode] 

204 except KeyError: 

205 raise ValueError("Acceptable mode flags are 'valid'," 

206 " 'same', or 'full'.") 

207 

208 # this either calls fftconvolve or this function with method=='direct' 

209 if method in ('fft', 'auto'): 

210 return convolve(in1, _reverse_and_conj(in2), mode, method) 

211 

212 elif method == 'direct': 

213 # fastpath to faster numpy.correlate for 1d inputs when possible 

214 if _np_conv_ok(in1, in2, mode): 

215 return np.correlate(in1, in2, mode) 

216 

217 # _correlateND is far slower when in2.size > in1.size, so swap them 

218 # and then undo the effect afterward if mode == 'full'. Also, it fails 

219 # with 'valid' mode if in2 is larger than in1, so swap those, too. 

220 # Don't swap inputs for 'same' mode, since shape of in1 matters. 

221 swapped_inputs = ((mode == 'full') and (in2.size > in1.size) or 

222 _inputs_swap_needed(mode, in1.shape, in2.shape)) 

223 

224 if swapped_inputs: 

225 in1, in2 = in2, in1 

226 

227 if mode == 'valid': 

228 ps = [i - j + 1 for i, j in zip(in1.shape, in2.shape)] 

229 out = np.empty(ps, in1.dtype) 

230 

231 z = sigtools._correlateND(in1, in2, out, val) 

232 

233 else: 

234 ps = [i + j - 1 for i, j in zip(in1.shape, in2.shape)] 

235 

236 # zero pad input 

237 in1zpadded = np.zeros(ps, in1.dtype) 

238 sc = tuple(slice(0, i) for i in in1.shape) 

239 in1zpadded[sc] = in1.copy() 

240 

241 if mode == 'full': 

242 out = np.empty(ps, in1.dtype) 

243 elif mode == 'same': 

244 out = np.empty(in1.shape, in1.dtype) 

245 

246 z = sigtools._correlateND(in1zpadded, in2, out, val) 

247 

248 if swapped_inputs: 

249 # Reverse and conjugate to undo the effect of swapping inputs 

250 z = _reverse_and_conj(z) 

251 

252 return z 

253 

254 else: 

255 raise ValueError("Acceptable method flags are 'auto'," 

256 " 'direct', or 'fft'.") 

257 

258 

259def _centered(arr, newshape): 

260 # Return the center newshape portion of the array. 

261 newshape = np.asarray(newshape) 

262 currshape = np.array(arr.shape) 

263 startind = (currshape - newshape) // 2 

264 endind = startind + newshape 

265 myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] 

266 return arr[tuple(myslice)] 

267 

268 

269def _init_freq_conv_axes(in1, in2, mode, axes, sorted_axes=False): 

270 """Handle the axes argument for frequency-domain convolution. 

271 

272 Returns the inputs and axes in a standard form, eliminating redundant axes, 

273 swapping the inputs if necessary, and checking for various potential 

274 errors. 

275 

276 Parameters 

277 ---------- 

278 in1 : array 

279 First input. 

280 in2 : array 

281 Second input. 

282 mode : str {'full', 'valid', 'same'}, optional 

283 A string indicating the size of the output. 

284 See the documentation `fftconvolve` for more information. 

285 axes : list of ints 

286 Axes over which to compute the FFTs. 

287 sorted_axes : bool, optional 

288 If `True`, sort the axes. 

289 Default is `False`, do not sort. 

290 

291 Returns 

292 ------- 

293 in1 : array 

294 The first input, possible swapped with the second input. 

295 in2 : array 

296 The second input, possible swapped with the first input. 

297 axes : list of ints 

298 Axes over which to compute the FFTs. 

299 

300 """ 

301 s1 = in1.shape 

302 s2 = in2.shape 

303 noaxes = axes is None 

304 

305 _, axes = _init_nd_shape_and_axes(in1, shape=None, axes=axes) 

306 

307 if not noaxes and not len(axes): 

308 raise ValueError("when provided, axes cannot be empty") 

309 

310 # Axes of length 1 can rely on broadcasting rules for multipy, 

311 # no fft needed. 

312 axes = [a for a in axes if s1[a] != 1 and s2[a] != 1] 

313 

314 if sorted_axes: 

315 axes.sort() 

316 

317 if not all(s1[a] == s2[a] or s1[a] == 1 or s2[a] == 1 

318 for a in range(in1.ndim) if a not in axes): 

319 raise ValueError("incompatible shapes for in1 and in2:" 

320 " {0} and {1}".format(s1, s2)) 

321 

322 # Check that input sizes are compatible with 'valid' mode. 

323 if _inputs_swap_needed(mode, s1, s2, axes=axes): 

324 # Convolution is commutative; order doesn't have any effect on output. 

325 in1, in2 = in2, in1 

326 

327 return in1, in2, axes 

328 

329 

330def _freq_domain_conv(in1, in2, axes, shape, calc_fast_len=False): 

331 """Convolve two arrays in the frequency domain. 

332 

333 This function implements only base the FFT-related operations. 

334 Specifically, it converts the signals to the frequency domain, multiplies 

335 them, then converts them back to the time domain. Calculations of axes, 

336 shapes, convolution mode, etc. are implemented in higher level-functions, 

337 such as `fftconvolve` and `oaconvolve`. Those functions should be used 

338 instead of this one. 

339 

340 Parameters 

341 ---------- 

342 in1 : array_like 

343 First input. 

344 in2 : array_like 

345 Second input. Should have the same number of dimensions as `in1`. 

346 axes : array_like of ints 

347 Axes over which to compute the FFTs. 

348 shape : array_like of ints 

349 The sizes of the FFTs. 

350 calc_fast_len : bool, optional 

351 If `True`, set each value of `shape` to the next fast FFT length. 

352 Default is `False`, use `axes` as-is. 

353 

354 Returns 

355 ------- 

356 out : array 

357 An N-dimensional array containing the discrete linear convolution of 

358 `in1` with `in2`. 

359 

360 """ 

361 if not len(axes): 

362 return in1 * in2 

363 

364 complex_result = (in1.dtype.kind == 'c' or in2.dtype.kind == 'c') 

365 

366 if calc_fast_len: 

367 # Speed up FFT by padding to optimal size. 

368 fshape = [ 

369 sp_fft.next_fast_len(shape[a], not complex_result) for a in axes] 

370 else: 

371 fshape = shape 

372 

373 if not complex_result: 

374 fft, ifft = sp_fft.rfftn, sp_fft.irfftn 

375 else: 

376 fft, ifft = sp_fft.fftn, sp_fft.ifftn 

377 

378 sp1 = fft(in1, fshape, axes=axes) 

379 sp2 = fft(in2, fshape, axes=axes) 

380 

381 ret = ifft(sp1 * sp2, fshape, axes=axes) 

382 

383 if calc_fast_len: 

384 fslice = tuple([slice(sz) for sz in shape]) 

385 ret = ret[fslice] 

386 

387 return ret 

388 

389 

390def _apply_conv_mode(ret, s1, s2, mode, axes): 

391 """Calculate the convolution result shape based on the `mode` argument. 

392 

393 Returns the result sliced to the correct size for the given mode. 

394 

395 Parameters 

396 ---------- 

397 ret : array 

398 The result array, with the appropriate shape for the 'full' mode. 

399 s1 : list of int 

400 The shape of the first input. 

401 s2 : list of int 

402 The shape of the second input. 

403 mode : str {'full', 'valid', 'same'} 

404 A string indicating the size of the output. 

405 See the documentation `fftconvolve` for more information. 

406 axes : list of ints 

407 Axes over which to compute the convolution. 

408 

409 Returns 

410 ------- 

411 ret : array 

412 A copy of `res`, sliced to the correct size for the given `mode`. 

413 

414 """ 

415 if mode == "full": 

416 return ret.copy() 

417 elif mode == "same": 

418 return _centered(ret, s1).copy() 

419 elif mode == "valid": 

420 shape_valid = [ret.shape[a] if a not in axes else s1[a] - s2[a] + 1 

421 for a in range(ret.ndim)] 

422 return _centered(ret, shape_valid).copy() 

423 else: 

424 raise ValueError("acceptable mode flags are 'valid'," 

425 " 'same', or 'full'") 

426 

427 

428def fftconvolve(in1, in2, mode="full", axes=None): 

429 """Convolve two N-dimensional arrays using FFT. 

430 

431 Convolve `in1` and `in2` using the fast Fourier transform method, with 

432 the output size determined by the `mode` argument. 

433 

434 This is generally much faster than `convolve` for large arrays (n > ~500), 

435 but can be slower when only a few output values are needed, and can only 

436 output float arrays (int or object array inputs will be cast to float). 

437 

438 As of v0.19, `convolve` automatically chooses this method or the direct 

439 method based on an estimation of which is faster. 

440 

441 Parameters 

442 ---------- 

443 in1 : array_like 

444 First input. 

445 in2 : array_like 

446 Second input. Should have the same number of dimensions as `in1`. 

447 mode : str {'full', 'valid', 'same'}, optional 

448 A string indicating the size of the output: 

449 

450 ``full`` 

451 The output is the full discrete linear convolution 

452 of the inputs. (Default) 

453 ``valid`` 

454 The output consists only of those elements that do not 

455 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 

456 must be at least as large as the other in every dimension. 

457 ``same`` 

458 The output is the same size as `in1`, centered 

459 with respect to the 'full' output. 

460 axes : int or array_like of ints or None, optional 

461 Axes over which to compute the convolution. 

462 The default is over all axes. 

463 

464 Returns 

465 ------- 

466 out : array 

467 An N-dimensional array containing a subset of the discrete linear 

468 convolution of `in1` with `in2`. 

469 

470 See Also 

471 -------- 

472 convolve : Uses the direct convolution or FFT convolution algorithm 

473 depending on which is faster. 

474 oaconvolve : Uses the overlap-add method to do convolution, which is 

475 generally faster when the input arrays are large and 

476 significantly different in size. 

477 

478 Examples 

479 -------- 

480 Autocorrelation of white noise is an impulse. 

481 

482 >>> from scipy import signal 

483 >>> sig = np.random.randn(1000) 

484 >>> autocorr = signal.fftconvolve(sig, sig[::-1], mode='full') 

485 

486 >>> import matplotlib.pyplot as plt 

487 >>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1) 

488 >>> ax_orig.plot(sig) 

489 >>> ax_orig.set_title('White noise') 

490 >>> ax_mag.plot(np.arange(-len(sig)+1,len(sig)), autocorr) 

491 >>> ax_mag.set_title('Autocorrelation') 

492 >>> fig.tight_layout() 

493 >>> fig.show() 

494 

495 Gaussian blur implemented using FFT convolution. Notice the dark borders 

496 around the image, due to the zero-padding beyond its boundaries. 

497 The `convolve2d` function allows for other types of image boundaries, 

498 but is far slower. 

499 

500 >>> from scipy import misc 

501 >>> face = misc.face(gray=True) 

502 >>> kernel = np.outer(signal.gaussian(70, 8), signal.gaussian(70, 8)) 

503 >>> blurred = signal.fftconvolve(face, kernel, mode='same') 

504 

505 >>> fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1, 

506 ... figsize=(6, 15)) 

507 >>> ax_orig.imshow(face, cmap='gray') 

508 >>> ax_orig.set_title('Original') 

509 >>> ax_orig.set_axis_off() 

510 >>> ax_kernel.imshow(kernel, cmap='gray') 

511 >>> ax_kernel.set_title('Gaussian kernel') 

512 >>> ax_kernel.set_axis_off() 

513 >>> ax_blurred.imshow(blurred, cmap='gray') 

514 >>> ax_blurred.set_title('Blurred') 

515 >>> ax_blurred.set_axis_off() 

516 >>> fig.show() 

517 

518 """ 

519 in1 = np.asarray(in1) 

520 in2 = np.asarray(in2) 

521 

522 if in1.ndim == in2.ndim == 0: # scalar inputs 

523 return in1 * in2 

524 elif in1.ndim != in2.ndim: 

525 raise ValueError("in1 and in2 should have the same dimensionality") 

526 elif in1.size == 0 or in2.size == 0: # empty arrays 

527 return np.array([]) 

528 

529 in1, in2, axes = _init_freq_conv_axes(in1, in2, mode, axes, 

530 sorted_axes=False) 

531 

532 s1 = in1.shape 

533 s2 = in2.shape 

534 

535 shape = [max((s1[i], s2[i])) if i not in axes else s1[i] + s2[i] - 1 

536 for i in range(in1.ndim)] 

537 

538 ret = _freq_domain_conv(in1, in2, axes, shape, calc_fast_len=True) 

539 

540 return _apply_conv_mode(ret, s1, s2, mode, axes) 

541 

542 

543def _calc_oa_lens(s1, s2): 

544 """Calculate the optimal FFT lengths for overlapp-add convolution. 

545 

546 The calculation is done for a single dimension. 

547 

548 Parameters 

549 ---------- 

550 s1 : int 

551 Size of the dimension for the first array. 

552 s2 : int 

553 Size of the dimension for the second array. 

554 

555 Returns 

556 ------- 

557 block_size : int 

558 The size of the FFT blocks. 

559 overlap : int 

560 The amount of overlap between two blocks. 

561 in1_step : int 

562 The size of each step for the first array. 

563 in2_step : int 

564 The size of each step for the first array. 

565 

566 """ 

567 # Set up the arguments for the conventional FFT approach. 

568 fallback = (s1+s2-1, None, s1, s2) 

569 

570 # Use conventional FFT convolve if sizes are same. 

571 if s1 == s2 or s1 == 1 or s2 == 1: 

572 return fallback 

573 

574 if s2 > s1: 

575 s1, s2 = s2, s1 

576 swapped = True 

577 else: 

578 swapped = False 

579 

580 # There cannot be a useful block size if s2 is more than half of s1. 

581 if s2 >= s1/2: 

582 return fallback 

583 

584 # Derivation of optimal block length 

585 # For original formula see: 

586 # https://en.wikipedia.org/wiki/Overlap-add_method 

587 # 

588 # Formula: 

589 # K = overlap = s2-1 

590 # N = block_size 

591 # C = complexity 

592 # e = exponential, exp(1) 

593 # 

594 # C = (N*(log2(N)+1))/(N-K) 

595 # C = (N*log2(2N))/(N-K) 

596 # C = N/(N-K) * log2(2N) 

597 # C1 = N/(N-K) 

598 # C2 = log2(2N) = ln(2N)/ln(2) 

599 # 

600 # dC1/dN = (1*(N-K)-N)/(N-K)^2 = -K/(N-K)^2 

601 # dC2/dN = 2/(2*N*ln(2)) = 1/(N*ln(2)) 

602 # 

603 # dC/dN = dC1/dN*C2 + dC2/dN*C1 

604 # dC/dN = -K*ln(2N)/(ln(2)*(N-K)^2) + N/(N*ln(2)*(N-K)) 

605 # dC/dN = -K*ln(2N)/(ln(2)*(N-K)^2) + 1/(ln(2)*(N-K)) 

606 # dC/dN = -K*ln(2N)/(ln(2)*(N-K)^2) + (N-K)/(ln(2)*(N-K)^2) 

607 # dC/dN = (-K*ln(2N) + (N-K)/(ln(2)*(N-K)^2) 

608 # dC/dN = (N - K*ln(2N) - K)/(ln(2)*(N-K)^2) 

609 # 

610 # Solve for minimum, where dC/dN = 0 

611 # 0 = (N - K*ln(2N) - K)/(ln(2)*(N-K)^2) 

612 # 0 * ln(2)*(N-K)^2 = N - K*ln(2N) - K 

613 # 0 = N - K*ln(2N) - K 

614 # 0 = N - K*(ln(2N) + 1) 

615 # 0 = N - K*ln(2Ne) 

616 # N = K*ln(2Ne) 

617 # N/K = ln(2Ne) 

618 # 

619 # e^(N/K) = e^ln(2Ne) 

620 # e^(N/K) = 2Ne 

621 # 1/e^(N/K) = 1/(2*N*e) 

622 # e^(N/-K) = 1/(2*N*e) 

623 # e^(N/-K) = K/N*1/(2*K*e) 

624 # N/K*e^(N/-K) = 1/(2*e*K) 

625 # N/-K*e^(N/-K) = -1/(2*e*K) 

626 # 

627 # Using Lambert W function 

628 # https://en.wikipedia.org/wiki/Lambert_W_function 

629 # x = W(y) It is the solution to y = x*e^x 

630 # x = N/-K 

631 # y = -1/(2*e*K) 

632 # 

633 # N/-K = W(-1/(2*e*K)) 

634 # 

635 # N = -K*W(-1/(2*e*K)) 

636 overlap = s2-1 

637 opt_size = -overlap*lambertw(-1/(2*math.e*overlap), k=-1).real 

638 block_size = sp_fft.next_fast_len(math.ceil(opt_size)) 

639 

640 # Use conventional FFT convolve if there is only going to be one block. 

641 if block_size >= s1: 

642 return fallback 

643 

644 if not swapped: 

645 in1_step = block_size-s2+1 

646 in2_step = s2 

647 else: 

648 in1_step = s2 

649 in2_step = block_size-s2+1 

650 

651 return block_size, overlap, in1_step, in2_step 

652 

653 

654def oaconvolve(in1, in2, mode="full", axes=None): 

655 """Convolve two N-dimensional arrays using the overlap-add method. 

656 

657 Convolve `in1` and `in2` using the overlap-add method, with 

658 the output size determined by the `mode` argument. 

659 

660 This is generally much faster than `convolve` for large arrays (n > ~500), 

661 and generally much faster than `fftconvolve` when one array is much 

662 larger than the other, but can be slower when only a few output values are 

663 needed or when the arrays are very similar in shape, and can only 

664 output float arrays (int or object array inputs will be cast to float). 

665 

666 Parameters 

667 ---------- 

668 in1 : array_like 

669 First input. 

670 in2 : array_like 

671 Second input. Should have the same number of dimensions as `in1`. 

672 mode : str {'full', 'valid', 'same'}, optional 

673 A string indicating the size of the output: 

674 

675 ``full`` 

676 The output is the full discrete linear convolution 

677 of the inputs. (Default) 

678 ``valid`` 

679 The output consists only of those elements that do not 

680 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 

681 must be at least as large as the other in every dimension. 

682 ``same`` 

683 The output is the same size as `in1`, centered 

684 with respect to the 'full' output. 

685 axes : int or array_like of ints or None, optional 

686 Axes over which to compute the convolution. 

687 The default is over all axes. 

688 

689 Returns 

690 ------- 

691 out : array 

692 An N-dimensional array containing a subset of the discrete linear 

693 convolution of `in1` with `in2`. 

694 

695 See Also 

696 -------- 

697 convolve : Uses the direct convolution or FFT convolution algorithm 

698 depending on which is faster. 

699 fftconvolve : An implementation of convolution using FFT. 

700 

701 Notes 

702 ----- 

703 .. versionadded:: 1.4.0 

704 

705 Examples 

706 -------- 

707 Convolve a 100,000 sample signal with a 512-sample filter. 

708 

709 >>> from scipy import signal 

710 >>> sig = np.random.randn(100000) 

711 >>> filt = signal.firwin(512, 0.01) 

712 >>> fsig = signal.oaconvolve(sig, filt) 

713 

714 >>> import matplotlib.pyplot as plt 

715 >>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1) 

716 >>> ax_orig.plot(sig) 

717 >>> ax_orig.set_title('White noise') 

718 >>> ax_mag.plot(fsig) 

719 >>> ax_mag.set_title('Filtered noise') 

720 >>> fig.tight_layout() 

721 >>> fig.show() 

722 

723 References 

724 ---------- 

725 .. [1] Wikipedia, "Overlap-add_method". 

726 https://en.wikipedia.org/wiki/Overlap-add_method 

727 .. [2] Richard G. Lyons. Understanding Digital Signal Processing, 

728 Third Edition, 2011. Chapter 13.10. 

729 ISBN 13: 978-0137-02741-5 

730 

731 """ 

732 in1 = np.asarray(in1) 

733 in2 = np.asarray(in2) 

734 

735 if in1.ndim == in2.ndim == 0: # scalar inputs 

736 return in1 * in2 

737 elif in1.ndim != in2.ndim: 

738 raise ValueError("in1 and in2 should have the same dimensionality") 

739 elif in1.size == 0 or in2.size == 0: # empty arrays 

740 return np.array([]) 

741 elif in1.shape == in2.shape: # Equivalent to fftconvolve 

742 return fftconvolve(in1, in2, mode=mode, axes=axes) 

743 

744 in1, in2, axes = _init_freq_conv_axes(in1, in2, mode, axes, 

745 sorted_axes=True) 

746 

747 s1 = in1.shape 

748 s2 = in2.shape 

749 

750 if not axes: 

751 ret = in1 * in2 

752 return _apply_conv_mode(ret, s1, s2, mode, axes) 

753 

754 # Calculate this now since in1 is changed later 

755 shape_final = [None if i not in axes else 

756 s1[i] + s2[i] - 1 for i in range(in1.ndim)] 

757 

758 # Calculate the block sizes for the output, steps, first and second inputs. 

759 # It is simpler to calculate them all together than doing them in separate 

760 # loops due to all the special cases that need to be handled. 

761 optimal_sizes = ((-1, -1, s1[i], s2[i]) if i not in axes else 

762 _calc_oa_lens(s1[i], s2[i]) for i in range(in1.ndim)) 

763 block_size, overlaps, \ 

764 in1_step, in2_step = zip(*optimal_sizes) 

765 

766 # Fall back to fftconvolve if there is only one block in every dimension. 

767 if in1_step == s1 and in2_step == s2: 

768 return fftconvolve(in1, in2, mode=mode, axes=axes) 

769 

770 # Figure out the number of steps and padding. 

771 # This would get too complicated in a list comprehension. 

772 nsteps1 = [] 

773 nsteps2 = [] 

774 pad_size1 = [] 

775 pad_size2 = [] 

776 for i in range(in1.ndim): 

777 if i not in axes: 

778 pad_size1 += [(0, 0)] 

779 pad_size2 += [(0, 0)] 

780 continue 

781 

782 if s1[i] > in1_step[i]: 

783 curnstep1 = math.ceil((s1[i]+1)/in1_step[i]) 

784 if (block_size[i] - overlaps[i])*curnstep1 < shape_final[i]: 

785 curnstep1 += 1 

786 

787 curpad1 = curnstep1*in1_step[i] - s1[i] 

788 else: 

789 curnstep1 = 1 

790 curpad1 = 0 

791 

792 if s2[i] > in2_step[i]: 

793 curnstep2 = math.ceil((s2[i]+1)/in2_step[i]) 

794 if (block_size[i] - overlaps[i])*curnstep2 < shape_final[i]: 

795 curnstep2 += 1 

796 

797 curpad2 = curnstep2*in2_step[i] - s2[i] 

798 else: 

799 curnstep2 = 1 

800 curpad2 = 0 

801 

802 nsteps1 += [curnstep1] 

803 nsteps2 += [curnstep2] 

804 pad_size1 += [(0, curpad1)] 

805 pad_size2 += [(0, curpad2)] 

806 

807 # Pad the array to a size that can be reshaped to the desired shape 

808 # if necessary. 

809 if not all(curpad == (0, 0) for curpad in pad_size1): 

810 in1 = np.pad(in1, pad_size1, mode='constant', constant_values=0) 

811 

812 if not all(curpad == (0, 0) for curpad in pad_size2): 

813 in2 = np.pad(in2, pad_size2, mode='constant', constant_values=0) 

814 

815 # Reshape the overlap-add parts to input block sizes. 

816 split_axes = [iax+i for i, iax in enumerate(axes)] 

817 fft_axes = [iax+1 for iax in split_axes] 

818 

819 # We need to put each new dimension before the corresponding dimension 

820 # being reshaped in order to get the data in the right layout at the end. 

821 reshape_size1 = list(in1_step) 

822 reshape_size2 = list(in2_step) 

823 for i, iax in enumerate(split_axes): 

824 reshape_size1.insert(iax, nsteps1[i]) 

825 reshape_size2.insert(iax, nsteps2[i]) 

826 

827 in1 = in1.reshape(*reshape_size1) 

828 in2 = in2.reshape(*reshape_size2) 

829 

830 # Do the convolution. 

831 fft_shape = [block_size[i] for i in axes] 

832 ret = _freq_domain_conv(in1, in2, fft_axes, fft_shape, calc_fast_len=False) 

833 

834 # Do the overlap-add. 

835 for ax, ax_fft, ax_split in zip(axes, fft_axes, split_axes): 

836 overlap = overlaps[ax] 

837 if overlap is None: 

838 continue 

839 

840 ret, overpart = np.split(ret, [-overlap], ax_fft) 

841 overpart = np.split(overpart, [-1], ax_split)[0] 

842 

843 ret_overpart = np.split(ret, [overlap], ax_fft)[0] 

844 ret_overpart = np.split(ret_overpart, [1], ax_split)[1] 

845 ret_overpart += overpart 

846 

847 # Reshape back to the correct dimensionality. 

848 shape_ret = [ret.shape[i] if i not in fft_axes else 

849 ret.shape[i]*ret.shape[i-1] 

850 for i in range(ret.ndim) if i not in split_axes] 

851 ret = ret.reshape(*shape_ret) 

852 

853 # Slice to the correct size. 

854 slice_final = tuple([slice(islice) for islice in shape_final]) 

855 ret = ret[slice_final] 

856 

857 return _apply_conv_mode(ret, s1, s2, mode, axes) 

858 

859 

860def _numeric_arrays(arrays, kinds='buifc'): 

861 """ 

862 See if a list of arrays are all numeric. 

863 

864 Parameters 

865 ---------- 

866 ndarrays : array or list of arrays 

867 arrays to check if numeric. 

868 numeric_kinds : string-like 

869 The dtypes of the arrays to be checked. If the dtype.kind of 

870 the ndarrays are not in this string the function returns False and 

871 otherwise returns True. 

872 """ 

873 if type(arrays) == np.ndarray: 

874 return arrays.dtype.kind in kinds 

875 for array_ in arrays: 

876 if array_.dtype.kind not in kinds: 

877 return False 

878 return True 

879 

880 

881def _conv_ops(x_shape, h_shape, mode): 

882 """ 

883 Find the number of operations required for direct/fft methods of 

884 convolution. The direct operations were recorded by making a dummy class to 

885 record the number of operations by overriding ``__mul__`` and ``__add__``. 

886 The FFT operations rely on the (well-known) computational complexity of the 

887 FFT (and the implementation of ``_freq_domain_conv``). 

888 

889 """ 

890 if mode == "full": 

891 out_shape = [n + k - 1 for n, k in zip(x_shape, h_shape)] 

892 elif mode == "valid": 

893 out_shape = [abs(n - k) + 1 for n, k in zip(x_shape, h_shape)] 

894 elif mode == "same": 

895 out_shape = x_shape 

896 else: 

897 raise ValueError("Acceptable mode flags are 'valid'," 

898 " 'same', or 'full', not mode={}".format(mode)) 

899 

900 s1, s2 = x_shape, h_shape 

901 if len(x_shape) == 1: 

902 s1, s2 = s1[0], s2[0] 

903 if mode == "full": 

904 direct_ops = s1 * s2 

905 elif mode == "valid": 

906 direct_ops = (s2 - s1 + 1) * s1 if s2 >= s1 else (s1 - s2 + 1) * s2 

907 elif mode == "same": 

908 direct_ops = (s1 * s2 if s1 < s2 else 

909 s1 * s2 - (s2 // 2) * ((s2 + 1) // 2)) 

910 else: 

911 if mode == "full": 

912 direct_ops = min(_prod(s1), _prod(s2)) * _prod(out_shape) 

913 elif mode == "valid": 

914 direct_ops = min(_prod(s1), _prod(s2)) * _prod(out_shape) 

915 elif mode == "same": 

916 direct_ops = _prod(s1) * _prod(s2) 

917 

918 full_out_shape = [n + k - 1 for n, k in zip(x_shape, h_shape)] 

919 N = _prod(full_out_shape) 

920 fft_ops = 3 * N * np.log(N) # 3 separate FFTs of size full_out_shape 

921 return fft_ops, direct_ops 

922 

923 

924def _fftconv_faster(x, h, mode): 

925 """ 

926 See if using fftconvolve or convolve is faster. 

927 

928 Parameters 

929 ---------- 

930 x : np.ndarray 

931 Signal 

932 h : np.ndarray 

933 Kernel 

934 mode : str 

935 Mode passed to convolve 

936 

937 Returns 

938 ------- 

939 fft_faster : bool 

940 

941 Notes 

942 ----- 

943 See docstring of `choose_conv_method` for details on tuning hardware. 

944 

945 See pull request 11031 for more detail: 

946 https://github.com/scipy/scipy/pull/11031. 

947 

948 """ 

949 fft_ops, direct_ops = _conv_ops(x.shape, h.shape, mode) 

950 offset = -1e-3 if x.ndim == 1 else -1e-4 

951 constants = { 

952 "valid": (1.89095737e-9, 2.1364985e-10, offset), 

953 "full": (1.7649070e-9, 2.1414831e-10, offset), 

954 "same": (3.2646654e-9, 2.8478277e-10, offset) 

955 if h.size <= x.size 

956 else (3.21635404e-9, 1.1773253e-8, -1e-5), 

957 } if x.ndim == 1 else { 

958 "valid": (1.85927e-9, 2.11242e-8, offset), 

959 "full": (1.99817e-9, 1.66174e-8, offset), 

960 "same": (2.04735e-9, 1.55367e-8, offset), 

961 } 

962 O_fft, O_direct, O_offset = constants[mode] 

963 return O_fft * fft_ops < O_direct * direct_ops + O_offset 

964 

965 

966def _reverse_and_conj(x): 

967 """ 

968 Reverse array `x` in all dimensions and perform the complex conjugate 

969 """ 

970 reverse = (slice(None, None, -1),) * x.ndim 

971 return x[reverse].conj() 

972 

973 

974def _np_conv_ok(volume, kernel, mode): 

975 """ 

976 See if numpy supports convolution of `volume` and `kernel` (i.e. both are 

977 1D ndarrays and of the appropriate shape). NumPy's 'same' mode uses the 

978 size of the larger input, while SciPy's uses the size of the first input. 

979 

980 Invalid mode strings will return False and be caught by the calling func. 

981 """ 

982 if volume.ndim == kernel.ndim == 1: 

983 if mode in ('full', 'valid'): 

984 return True 

985 elif mode == 'same': 

986 return volume.size >= kernel.size 

987 else: 

988 return False 

989 

990 

991def _timeit_fast(stmt="pass", setup="pass", repeat=3): 

992 """ 

993 Returns the time the statement/function took, in seconds. 

994 

995 Faster, less precise version of IPython's timeit. `stmt` can be a statement 

996 written as a string or a callable. 

997 

998 Will do only 1 loop (like IPython's timeit) with no repetitions 

999 (unlike IPython) for very slow functions. For fast functions, only does 

1000 enough loops to take 5 ms, which seems to produce similar results (on 

1001 Windows at least), and avoids doing an extraneous cycle that isn't 

1002 measured. 

1003 

1004 """ 

1005 timer = timeit.Timer(stmt, setup) 

1006 

1007 # determine number of calls per rep so total time for 1 rep >= 5 ms 

1008 x = 0 

1009 for p in range(0, 10): 

1010 number = 10**p 

1011 x = timer.timeit(number) # seconds 

1012 if x >= 5e-3 / 10: # 5 ms for final test, 1/10th that for this one 

1013 break 

1014 if x > 1: # second 

1015 # If it's macroscopic, don't bother with repetitions 

1016 best = x 

1017 else: 

1018 number *= 10 

1019 r = timer.repeat(repeat, number) 

1020 best = min(r) 

1021 

1022 sec = best / number 

1023 return sec 

1024 

1025 

1026def choose_conv_method(in1, in2, mode='full', measure=False): 

1027 """ 

1028 Find the fastest convolution/correlation method. 

1029 

1030 This primarily exists to be called during the ``method='auto'`` option in 

1031 `convolve` and `correlate`. It can also be used to determine the value of 

1032 ``method`` for many different convolutions of the same dtype/shape. 

1033 In addition, it supports timing the convolution to adapt the value of 

1034 ``method`` to a particular set of inputs and/or hardware. 

1035 

1036 Parameters 

1037 ---------- 

1038 in1 : array_like 

1039 The first argument passed into the convolution function. 

1040 in2 : array_like 

1041 The second argument passed into the convolution function. 

1042 mode : str {'full', 'valid', 'same'}, optional 

1043 A string indicating the size of the output: 

1044 

1045 ``full`` 

1046 The output is the full discrete linear convolution 

1047 of the inputs. (Default) 

1048 ``valid`` 

1049 The output consists only of those elements that do not 

1050 rely on the zero-padding. 

1051 ``same`` 

1052 The output is the same size as `in1`, centered 

1053 with respect to the 'full' output. 

1054 measure : bool, optional 

1055 If True, run and time the convolution of `in1` and `in2` with both 

1056 methods and return the fastest. If False (default), predict the fastest 

1057 method using precomputed values. 

1058 

1059 Returns 

1060 ------- 

1061 method : str 

1062 A string indicating which convolution method is fastest, either 

1063 'direct' or 'fft' 

1064 times : dict, optional 

1065 A dictionary containing the times (in seconds) needed for each method. 

1066 This value is only returned if ``measure=True``. 

1067 

1068 See Also 

1069 -------- 

1070 convolve 

1071 correlate 

1072 

1073 Notes 

1074 ----- 

1075 Generally, this method is 99% accurate for 2D signals and 85% accurate 

1076 for 1D signals for randomly chosen input sizes. For precision, use 

1077 ``measure=True`` to find the fastest method by timing the convolution. 

1078 This can be used to avoid the minimal overhead of finding the fastest 

1079 ``method`` later, or to adapt the value of ``method`` to a particular set 

1080 of inputs. 

1081 

1082 Experiments were run on an Amazon EC2 r5a.2xlarge machine to test this 

1083 function. These experiments measured the ratio between the time required 

1084 when using ``method='auto'`` and the time required for the fastest method 

1085 (i.e., ``ratio = time_auto / min(time_fft, time_direct)``). In these 

1086 experiments, we found: 

1087 

1088 * There is a 95% chance of this ratio being less than 1.5 for 1D signals 

1089 and a 99% chance of being less than 2.5 for 2D signals. 

1090 * The ratio was always less than 2.5/5 for 1D/2D signals respectively. 

1091 * This function is most inaccurate for 1D convolutions that take between 1 

1092 and 10 milliseconds with ``method='direct'``. A good proxy for this 

1093 (at least in our experiments) is ``1e6 <= in1.size * in2.size <= 1e7``. 

1094 

1095 The 2D results almost certainly generalize to 3D/4D/etc because the 

1096 implementation is the same (the 1D implementation is different). 

1097 

1098 All the numbers above are specific to the EC2 machine. However, we did find 

1099 that this function generalizes fairly decently across hardware. The speed 

1100 tests were of similar quality (and even slightly better) than the same 

1101 tests performed on the machine to tune this function's numbers (a mid-2014 

1102 15-inch MacBook Pro with 16GB RAM and a 2.5GHz Intel i7 processor). 

1103 

1104 There are cases when `fftconvolve` supports the inputs but this function 

1105 returns `direct` (e.g., to protect against floating point integer 

1106 precision). 

1107 

1108 .. versionadded:: 0.19 

1109 

1110 Examples 

1111 -------- 

1112 Estimate the fastest method for a given input: 

1113 

1114 >>> from scipy import signal 

1115 >>> img = np.random.rand(32, 32) 

1116 >>> filter = np.random.rand(8, 8) 

1117 >>> method = signal.choose_conv_method(img, filter, mode='same') 

1118 >>> method 

1119 'fft' 

1120 

1121 This can then be applied to other arrays of the same dtype and shape: 

1122 

1123 >>> img2 = np.random.rand(32, 32) 

1124 >>> filter2 = np.random.rand(8, 8) 

1125 >>> corr2 = signal.correlate(img2, filter2, mode='same', method=method) 

1126 >>> conv2 = signal.convolve(img2, filter2, mode='same', method=method) 

1127 

1128 The output of this function (``method``) works with `correlate` and 

1129 `convolve`. 

1130 

1131 """ 

1132 volume = np.asarray(in1) 

1133 kernel = np.asarray(in2) 

1134 

1135 if measure: 

1136 times = {} 

1137 for method in ['fft', 'direct']: 

1138 times[method] = _timeit_fast(lambda: convolve(volume, kernel, 

1139 mode=mode, method=method)) 

1140 

1141 chosen_method = 'fft' if times['fft'] < times['direct'] else 'direct' 

1142 return chosen_method, times 

1143 

1144 # for integer input, 

1145 # catch when more precision required than float provides (representing an 

1146 # integer as float can lose precision in fftconvolve if larger than 2**52) 

1147 if any([_numeric_arrays([x], kinds='ui') for x in [volume, kernel]]): 

1148 max_value = int(np.abs(volume).max()) * int(np.abs(kernel).max()) 

1149 max_value *= int(min(volume.size, kernel.size)) 

1150 if max_value > 2**np.finfo('float').nmant - 1: 

1151 return 'direct' 

1152 

1153 if _numeric_arrays([volume, kernel], kinds='b'): 

1154 return 'direct' 

1155 

1156 if _numeric_arrays([volume, kernel]): 

1157 if _fftconv_faster(volume, kernel, mode): 

1158 return 'fft' 

1159 

1160 return 'direct' 

1161 

1162 

1163def convolve(in1, in2, mode='full', method='auto'): 

1164 """ 

1165 Convolve two N-dimensional arrays. 

1166 

1167 Convolve `in1` and `in2`, with the output size determined by the 

1168 `mode` argument. 

1169 

1170 Parameters 

1171 ---------- 

1172 in1 : array_like 

1173 First input. 

1174 in2 : array_like 

1175 Second input. Should have the same number of dimensions as `in1`. 

1176 mode : str {'full', 'valid', 'same'}, optional 

1177 A string indicating the size of the output: 

1178 

1179 ``full`` 

1180 The output is the full discrete linear convolution 

1181 of the inputs. (Default) 

1182 ``valid`` 

1183 The output consists only of those elements that do not 

1184 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 

1185 must be at least as large as the other in every dimension. 

1186 ``same`` 

1187 The output is the same size as `in1`, centered 

1188 with respect to the 'full' output. 

1189 method : str {'auto', 'direct', 'fft'}, optional 

1190 A string indicating which method to use to calculate the convolution. 

1191 

1192 ``direct`` 

1193 The convolution is determined directly from sums, the definition of 

1194 convolution. 

1195 ``fft`` 

1196 The Fourier Transform is used to perform the convolution by calling 

1197 `fftconvolve`. 

1198 ``auto`` 

1199 Automatically chooses direct or Fourier method based on an estimate 

1200 of which is faster (default). See Notes for more detail. 

1201 

1202 .. versionadded:: 0.19.0 

1203 

1204 Returns 

1205 ------- 

1206 convolve : array 

1207 An N-dimensional array containing a subset of the discrete linear 

1208 convolution of `in1` with `in2`. 

1209 

1210 See Also 

1211 -------- 

1212 numpy.polymul : performs polynomial multiplication (same operation, but 

1213 also accepts poly1d objects) 

1214 choose_conv_method : chooses the fastest appropriate convolution method 

1215 fftconvolve : Always uses the FFT method. 

1216 oaconvolve : Uses the overlap-add method to do convolution, which is 

1217 generally faster when the input arrays are large and 

1218 significantly different in size. 

1219 

1220 Notes 

1221 ----- 

1222 By default, `convolve` and `correlate` use ``method='auto'``, which calls 

1223 `choose_conv_method` to choose the fastest method using pre-computed 

1224 values (`choose_conv_method` can also measure real-world timing with a 

1225 keyword argument). Because `fftconvolve` relies on floating point numbers, 

1226 there are certain constraints that may force `method=direct` (more detail 

1227 in `choose_conv_method` docstring). 

1228 

1229 Examples 

1230 -------- 

1231 Smooth a square pulse using a Hann window: 

1232 

1233 >>> from scipy import signal 

1234 >>> sig = np.repeat([0., 1., 0.], 100) 

1235 >>> win = signal.hann(50) 

1236 >>> filtered = signal.convolve(sig, win, mode='same') / sum(win) 

1237 

1238 >>> import matplotlib.pyplot as plt 

1239 >>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True) 

1240 >>> ax_orig.plot(sig) 

1241 >>> ax_orig.set_title('Original pulse') 

1242 >>> ax_orig.margins(0, 0.1) 

1243 >>> ax_win.plot(win) 

1244 >>> ax_win.set_title('Filter impulse response') 

1245 >>> ax_win.margins(0, 0.1) 

1246 >>> ax_filt.plot(filtered) 

1247 >>> ax_filt.set_title('Filtered signal') 

1248 >>> ax_filt.margins(0, 0.1) 

1249 >>> fig.tight_layout() 

1250 >>> fig.show() 

1251 

1252 """ 

1253 volume = np.asarray(in1) 

1254 kernel = np.asarray(in2) 

1255 

1256 if volume.ndim == kernel.ndim == 0: 

1257 return volume * kernel 

1258 elif volume.ndim != kernel.ndim: 

1259 raise ValueError("volume and kernel should have the same " 

1260 "dimensionality") 

1261 

1262 if _inputs_swap_needed(mode, volume.shape, kernel.shape): 

1263 # Convolution is commutative; order doesn't have any effect on output 

1264 volume, kernel = kernel, volume 

1265 

1266 if method == 'auto': 

1267 method = choose_conv_method(volume, kernel, mode=mode) 

1268 

1269 if method == 'fft': 

1270 out = fftconvolve(volume, kernel, mode=mode) 

1271 result_type = np.result_type(volume, kernel) 

1272 if result_type.kind in {'u', 'i'}: 

1273 out = np.around(out) 

1274 return out.astype(result_type) 

1275 elif method == 'direct': 

1276 # fastpath to faster numpy.convolve for 1d inputs when possible 

1277 if _np_conv_ok(volume, kernel, mode): 

1278 return np.convolve(volume, kernel, mode) 

1279 

1280 return correlate(volume, _reverse_and_conj(kernel), mode, 'direct') 

1281 else: 

1282 raise ValueError("Acceptable method flags are 'auto'," 

1283 " 'direct', or 'fft'.") 

1284 

1285 

1286def order_filter(a, domain, rank): 

1287 """ 

1288 Perform an order filter on an N-D array. 

1289 

1290 Perform an order filter on the array in. The domain argument acts as a 

1291 mask centered over each pixel. The non-zero elements of domain are 

1292 used to select elements surrounding each input pixel which are placed 

1293 in a list. The list is sorted, and the output for that pixel is the 

1294 element corresponding to rank in the sorted list. 

1295 

1296 Parameters 

1297 ---------- 

1298 a : ndarray 

1299 The N-dimensional input array. 

1300 domain : array_like 

1301 A mask array with the same number of dimensions as `a`. 

1302 Each dimension should have an odd number of elements. 

1303 rank : int 

1304 A non-negative integer which selects the element from the 

1305 sorted list (0 corresponds to the smallest element, 1 is the 

1306 next smallest element, etc.). 

1307 

1308 Returns 

1309 ------- 

1310 out : ndarray 

1311 The results of the order filter in an array with the same 

1312 shape as `a`. 

1313 

1314 Examples 

1315 -------- 

1316 >>> from scipy import signal 

1317 >>> x = np.arange(25).reshape(5, 5) 

1318 >>> domain = np.identity(3) 

1319 >>> x 

1320 array([[ 0, 1, 2, 3, 4], 

1321 [ 5, 6, 7, 8, 9], 

1322 [10, 11, 12, 13, 14], 

1323 [15, 16, 17, 18, 19], 

1324 [20, 21, 22, 23, 24]]) 

1325 >>> signal.order_filter(x, domain, 0) 

1326 array([[ 0., 0., 0., 0., 0.], 

1327 [ 0., 0., 1., 2., 0.], 

1328 [ 0., 5., 6., 7., 0.], 

1329 [ 0., 10., 11., 12., 0.], 

1330 [ 0., 0., 0., 0., 0.]]) 

1331 >>> signal.order_filter(x, domain, 2) 

1332 array([[ 6., 7., 8., 9., 4.], 

1333 [ 11., 12., 13., 14., 9.], 

1334 [ 16., 17., 18., 19., 14.], 

1335 [ 21., 22., 23., 24., 19.], 

1336 [ 20., 21., 22., 23., 24.]]) 

1337 

1338 """ 

1339 domain = np.asarray(domain) 

1340 size = domain.shape 

1341 for k in range(len(size)): 

1342 if (size[k] % 2) != 1: 

1343 raise ValueError("Each dimension of domain argument " 

1344 " should have an odd number of elements.") 

1345 return sigtools._order_filterND(a, domain, rank) 

1346 

1347 

1348def medfilt(volume, kernel_size=None): 

1349 """ 

1350 Perform a median filter on an N-dimensional array. 

1351 

1352 Apply a median filter to the input array using a local window-size 

1353 given by `kernel_size`. The array will automatically be zero-padded. 

1354 

1355 Parameters 

1356 ---------- 

1357 volume : array_like 

1358 An N-dimensional input array. 

1359 kernel_size : array_like, optional 

1360 A scalar or an N-length list giving the size of the median filter 

1361 window in each dimension. Elements of `kernel_size` should be odd. 

1362 If `kernel_size` is a scalar, then this scalar is used as the size in 

1363 each dimension. Default size is 3 for each dimension. 

1364 

1365 Returns 

1366 ------- 

1367 out : ndarray 

1368 An array the same size as input containing the median filtered 

1369 result. 

1370 

1371 Warns 

1372 ----- 

1373 UserWarning 

1374 If array size is smaller than kernel size along any dimension 

1375 

1376 See Also 

1377 -------- 

1378 scipy.ndimage.median_filter 

1379 

1380 Notes 

1381 ------- 

1382 The more general function `scipy.ndimage.median_filter` has a more 

1383 efficient implementation of a median filter and therefore runs much faster. 

1384 """ 

1385 volume = np.atleast_1d(volume) 

1386 if kernel_size is None: 

1387 kernel_size = [3] * volume.ndim 

1388 kernel_size = np.asarray(kernel_size) 

1389 if kernel_size.shape == (): 

1390 kernel_size = np.repeat(kernel_size.item(), volume.ndim) 

1391 

1392 for k in range(volume.ndim): 

1393 if (kernel_size[k] % 2) != 1: 

1394 raise ValueError("Each element of kernel_size should be odd.") 

1395 if any(k > s for k, s in zip(kernel_size, volume.shape)): 

1396 warnings.warn('kernel_size exceeds volume extent: the volume will be ' 

1397 'zero-padded.') 

1398 

1399 domain = np.ones(kernel_size) 

1400 

1401 numels = np.prod(kernel_size, axis=0) 

1402 order = numels // 2 

1403 return sigtools._order_filterND(volume, domain, order) 

1404 

1405 

1406def wiener(im, mysize=None, noise=None): 

1407 """ 

1408 Perform a Wiener filter on an N-dimensional array. 

1409 

1410 Apply a Wiener filter to the N-dimensional array `im`. 

1411 

1412 Parameters 

1413 ---------- 

1414 im : ndarray 

1415 An N-dimensional array. 

1416 mysize : int or array_like, optional 

1417 A scalar or an N-length list giving the size of the Wiener filter 

1418 window in each dimension. Elements of mysize should be odd. 

1419 If mysize is a scalar, then this scalar is used as the size 

1420 in each dimension. 

1421 noise : float, optional 

1422 The noise-power to use. If None, then noise is estimated as the 

1423 average of the local variance of the input. 

1424 

1425 Returns 

1426 ------- 

1427 out : ndarray 

1428 Wiener filtered result with the same shape as `im`. 

1429 

1430 Examples 

1431 -------- 

1432 

1433 >>> from scipy.misc import face 

1434 >>> from scipy.signal.signaltools import wiener 

1435 >>> import matplotlib.pyplot as plt 

1436 >>> import numpy as np 

1437 >>> img = np.random.random((40, 40)) #Create a random image 

1438 >>> filtered_img = wiener(img, (5, 5)) #Filter the image 

1439 >>> f, (plot1, plot2) = plt.subplots(1, 2) 

1440 >>> plot1.imshow(img) 

1441 >>> plot2.imshow(filtered_img) 

1442 >>> plt.show() 

1443 

1444 Notes 

1445 ----- 

1446 This implementation is similar to wiener2 in Matlab/Octave. 

1447 For more details see [1]_ 

1448 

1449 References 

1450 ---------- 

1451 .. [1] Lim, Jae S., Two-Dimensional Signal and Image Processing, 

1452 Englewood Cliffs, NJ, Prentice Hall, 1990, p. 548. 

1453 

1454 

1455 """ 

1456 im = np.asarray(im) 

1457 if mysize is None: 

1458 mysize = [3] * im.ndim 

1459 mysize = np.asarray(mysize) 

1460 if mysize.shape == (): 

1461 mysize = np.repeat(mysize.item(), im.ndim) 

1462 

1463 # Estimate the local mean 

1464 lMean = correlate(im, np.ones(mysize), 'same') / np.prod(mysize, axis=0) 

1465 

1466 # Estimate the local variance 

1467 lVar = (correlate(im ** 2, np.ones(mysize), 'same') / 

1468 np.prod(mysize, axis=0) - lMean ** 2) 

1469 

1470 # Estimate the noise power if needed. 

1471 if noise is None: 

1472 noise = np.mean(np.ravel(lVar), axis=0) 

1473 

1474 res = (im - lMean) 

1475 res *= (1 - noise / lVar) 

1476 res += lMean 

1477 out = np.where(lVar < noise, lMean, res) 

1478 

1479 return out 

1480 

1481 

1482def convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0): 

1483 """ 

1484 Convolve two 2-dimensional arrays. 

1485 

1486 Convolve `in1` and `in2` with output size determined by `mode`, and 

1487 boundary conditions determined by `boundary` and `fillvalue`. 

1488 

1489 Parameters 

1490 ---------- 

1491 in1 : array_like 

1492 First input. 

1493 in2 : array_like 

1494 Second input. Should have the same number of dimensions as `in1`. 

1495 mode : str {'full', 'valid', 'same'}, optional 

1496 A string indicating the size of the output: 

1497 

1498 ``full`` 

1499 The output is the full discrete linear convolution 

1500 of the inputs. (Default) 

1501 ``valid`` 

1502 The output consists only of those elements that do not 

1503 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 

1504 must be at least as large as the other in every dimension. 

1505 ``same`` 

1506 The output is the same size as `in1`, centered 

1507 with respect to the 'full' output. 

1508 boundary : str {'fill', 'wrap', 'symm'}, optional 

1509 A flag indicating how to handle boundaries: 

1510 

1511 ``fill`` 

1512 pad input arrays with fillvalue. (default) 

1513 ``wrap`` 

1514 circular boundary conditions. 

1515 ``symm`` 

1516 symmetrical boundary conditions. 

1517 

1518 fillvalue : scalar, optional 

1519 Value to fill pad input arrays with. Default is 0. 

1520 

1521 Returns 

1522 ------- 

1523 out : ndarray 

1524 A 2-dimensional array containing a subset of the discrete linear 

1525 convolution of `in1` with `in2`. 

1526 

1527 Examples 

1528 -------- 

1529 Compute the gradient of an image by 2D convolution with a complex Scharr 

1530 operator. (Horizontal operator is real, vertical is imaginary.) Use 

1531 symmetric boundary condition to avoid creating edges at the image 

1532 boundaries. 

1533 

1534 >>> from scipy import signal 

1535 >>> from scipy import misc 

1536 >>> ascent = misc.ascent() 

1537 >>> scharr = np.array([[ -3-3j, 0-10j, +3 -3j], 

1538 ... [-10+0j, 0+ 0j, +10 +0j], 

1539 ... [ -3+3j, 0+10j, +3 +3j]]) # Gx + j*Gy 

1540 >>> grad = signal.convolve2d(ascent, scharr, boundary='symm', mode='same') 

1541 

1542 >>> import matplotlib.pyplot as plt 

1543 >>> fig, (ax_orig, ax_mag, ax_ang) = plt.subplots(3, 1, figsize=(6, 15)) 

1544 >>> ax_orig.imshow(ascent, cmap='gray') 

1545 >>> ax_orig.set_title('Original') 

1546 >>> ax_orig.set_axis_off() 

1547 >>> ax_mag.imshow(np.absolute(grad), cmap='gray') 

1548 >>> ax_mag.set_title('Gradient magnitude') 

1549 >>> ax_mag.set_axis_off() 

1550 >>> ax_ang.imshow(np.angle(grad), cmap='hsv') # hsv is cyclic, like angles 

1551 >>> ax_ang.set_title('Gradient orientation') 

1552 >>> ax_ang.set_axis_off() 

1553 >>> fig.show() 

1554 

1555 """ 

1556 in1 = np.asarray(in1) 

1557 in2 = np.asarray(in2) 

1558 

1559 if not in1.ndim == in2.ndim == 2: 

1560 raise ValueError('convolve2d inputs must both be 2-D arrays') 

1561 

1562 if _inputs_swap_needed(mode, in1.shape, in2.shape): 

1563 in1, in2 = in2, in1 

1564 

1565 val = _valfrommode(mode) 

1566 bval = _bvalfromboundary(boundary) 

1567 out = sigtools._convolve2d(in1, in2, 1, val, bval, fillvalue) 

1568 return out 

1569 

1570 

1571def correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0): 

1572 """ 

1573 Cross-correlate two 2-dimensional arrays. 

1574 

1575 Cross correlate `in1` and `in2` with output size determined by `mode`, and 

1576 boundary conditions determined by `boundary` and `fillvalue`. 

1577 

1578 Parameters 

1579 ---------- 

1580 in1 : array_like 

1581 First input. 

1582 in2 : array_like 

1583 Second input. Should have the same number of dimensions as `in1`. 

1584 mode : str {'full', 'valid', 'same'}, optional 

1585 A string indicating the size of the output: 

1586 

1587 ``full`` 

1588 The output is the full discrete linear cross-correlation 

1589 of the inputs. (Default) 

1590 ``valid`` 

1591 The output consists only of those elements that do not 

1592 rely on the zero-padding. In 'valid' mode, either `in1` or `in2` 

1593 must be at least as large as the other in every dimension. 

1594 ``same`` 

1595 The output is the same size as `in1`, centered 

1596 with respect to the 'full' output. 

1597 boundary : str {'fill', 'wrap', 'symm'}, optional 

1598 A flag indicating how to handle boundaries: 

1599 

1600 ``fill`` 

1601 pad input arrays with fillvalue. (default) 

1602 ``wrap`` 

1603 circular boundary conditions. 

1604 ``symm`` 

1605 symmetrical boundary conditions. 

1606 

1607 fillvalue : scalar, optional 

1608 Value to fill pad input arrays with. Default is 0. 

1609 

1610 Returns 

1611 ------- 

1612 correlate2d : ndarray 

1613 A 2-dimensional array containing a subset of the discrete linear 

1614 cross-correlation of `in1` with `in2`. 

1615 

1616 Notes 

1617 ----- 

1618 When using "same" mode with even-length inputs, the outputs of `correlate` 

1619 and `correlate2d` differ: There is a 1-index offset between them. 

1620 

1621 Examples 

1622 -------- 

1623 Use 2D cross-correlation to find the location of a template in a noisy 

1624 image: 

1625 

1626 >>> from scipy import signal 

1627 >>> from scipy import misc 

1628 >>> face = misc.face(gray=True) - misc.face(gray=True).mean() 

1629 >>> template = np.copy(face[300:365, 670:750]) # right eye 

1630 >>> template -= template.mean() 

1631 >>> face = face + np.random.randn(*face.shape) * 50 # add noise 

1632 >>> corr = signal.correlate2d(face, template, boundary='symm', mode='same') 

1633 >>> y, x = np.unravel_index(np.argmax(corr), corr.shape) # find the match 

1634 

1635 >>> import matplotlib.pyplot as plt 

1636 >>> fig, (ax_orig, ax_template, ax_corr) = plt.subplots(3, 1, 

1637 ... figsize=(6, 15)) 

1638 >>> ax_orig.imshow(face, cmap='gray') 

1639 >>> ax_orig.set_title('Original') 

1640 >>> ax_orig.set_axis_off() 

1641 >>> ax_template.imshow(template, cmap='gray') 

1642 >>> ax_template.set_title('Template') 

1643 >>> ax_template.set_axis_off() 

1644 >>> ax_corr.imshow(corr, cmap='gray') 

1645 >>> ax_corr.set_title('Cross-correlation') 

1646 >>> ax_corr.set_axis_off() 

1647 >>> ax_orig.plot(x, y, 'ro') 

1648 >>> fig.show() 

1649 

1650 """ 

1651 in1 = np.asarray(in1) 

1652 in2 = np.asarray(in2) 

1653 

1654 if not in1.ndim == in2.ndim == 2: 

1655 raise ValueError('correlate2d inputs must both be 2-D arrays') 

1656 

1657 swapped_inputs = _inputs_swap_needed(mode, in1.shape, in2.shape) 

1658 if swapped_inputs: 

1659 in1, in2 = in2, in1 

1660 

1661 val = _valfrommode(mode) 

1662 bval = _bvalfromboundary(boundary) 

1663 out = sigtools._convolve2d(in1, in2.conj(), 0, val, bval, fillvalue) 

1664 

1665 if swapped_inputs: 

1666 out = out[::-1, ::-1] 

1667 

1668 return out 

1669 

1670 

1671def medfilt2d(input, kernel_size=3): 

1672 """ 

1673 Median filter a 2-dimensional array. 

1674 

1675 Apply a median filter to the `input` array using a local window-size 

1676 given by `kernel_size` (must be odd). The array is zero-padded 

1677 automatically. 

1678 

1679 Parameters 

1680 ---------- 

1681 input : array_like 

1682 A 2-dimensional input array. 

1683 kernel_size : array_like, optional 

1684 A scalar or a list of length 2, giving the size of the 

1685 median filter window in each dimension. Elements of 

1686 `kernel_size` should be odd. If `kernel_size` is a scalar, 

1687 then this scalar is used as the size in each dimension. 

1688 Default is a kernel of size (3, 3). 

1689 

1690 Returns 

1691 ------- 

1692 out : ndarray 

1693 An array the same size as input containing the median filtered 

1694 result. 

1695 

1696 See also 

1697 -------- 

1698 scipy.ndimage.median_filter 

1699 

1700 Notes 

1701 ------- 

1702 The more general function `scipy.ndimage.median_filter` has a more 

1703 efficient implementation of a median filter and therefore runs much faster. 

1704 """ 

1705 image = np.asarray(input) 

1706 if kernel_size is None: 

1707 kernel_size = [3] * 2 

1708 kernel_size = np.asarray(kernel_size) 

1709 if kernel_size.shape == (): 

1710 kernel_size = np.repeat(kernel_size.item(), 2) 

1711 

1712 for size in kernel_size: 

1713 if (size % 2) != 1: 

1714 raise ValueError("Each element of kernel_size should be odd.") 

1715 

1716 return sigtools._medfilt2d(image, kernel_size) 

1717 

1718 

1719def lfilter(b, a, x, axis=-1, zi=None): 

1720 """ 

1721 Filter data along one-dimension with an IIR or FIR filter. 

1722 

1723 Filter a data sequence, `x`, using a digital filter. This works for many 

1724 fundamental data types (including Object type). The filter is a direct 

1725 form II transposed implementation of the standard difference equation 

1726 (see Notes). 

1727 

1728 The function `sosfilt` (and filter design using ``output='sos'``) should be 

1729 preferred over `lfilter` for most filtering tasks, as second-order sections 

1730 have fewer numerical problems. 

1731 

1732 Parameters 

1733 ---------- 

1734 b : array_like 

1735 The numerator coefficient vector in a 1-D sequence. 

1736 a : array_like 

1737 The denominator coefficient vector in a 1-D sequence. If ``a[0]`` 

1738 is not 1, then both `a` and `b` are normalized by ``a[0]``. 

1739 x : array_like 

1740 An N-dimensional input array. 

1741 axis : int, optional 

1742 The axis of the input data array along which to apply the 

1743 linear filter. The filter is applied to each subarray along 

1744 this axis. Default is -1. 

1745 zi : array_like, optional 

1746 Initial conditions for the filter delays. It is a vector 

1747 (or array of vectors for an N-dimensional input) of length 

1748 ``max(len(a), len(b)) - 1``. If `zi` is None or is not given then 

1749 initial rest is assumed. See `lfiltic` for more information. 

1750 

1751 Returns 

1752 ------- 

1753 y : array 

1754 The output of the digital filter. 

1755 zf : array, optional 

1756 If `zi` is None, this is not returned, otherwise, `zf` holds the 

1757 final filter delay values. 

1758 

1759 See Also 

1760 -------- 

1761 lfiltic : Construct initial conditions for `lfilter`. 

1762 lfilter_zi : Compute initial state (steady state of step response) for 

1763 `lfilter`. 

1764 filtfilt : A forward-backward filter, to obtain a filter with linear phase. 

1765 savgol_filter : A Savitzky-Golay filter. 

1766 sosfilt: Filter data using cascaded second-order sections. 

1767 sosfiltfilt: A forward-backward filter using second-order sections. 

1768 

1769 Notes 

1770 ----- 

1771 The filter function is implemented as a direct II transposed structure. 

1772 This means that the filter implements:: 

1773 

1774 a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M] 

1775 - a[1]*y[n-1] - ... - a[N]*y[n-N] 

1776 

1777 where `M` is the degree of the numerator, `N` is the degree of the 

1778 denominator, and `n` is the sample number. It is implemented using 

1779 the following difference equations (assuming M = N):: 

1780 

1781 a[0]*y[n] = b[0] * x[n] + d[0][n-1] 

1782 d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n-1] 

1783 d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n-1] 

1784 ... 

1785 d[N-2][n] = b[N-1]*x[n] - a[N-1]*y[n] + d[N-1][n-1] 

1786 d[N-1][n] = b[N] * x[n] - a[N] * y[n] 

1787 

1788 where `d` are the state variables. 

1789 

1790 The rational transfer function describing this filter in the 

1791 z-transform domain is:: 

1792 

1793 -1 -M 

1794 b[0] + b[1]z + ... + b[M] z 

1795 Y(z) = -------------------------------- X(z) 

1796 -1 -N 

1797 a[0] + a[1]z + ... + a[N] z 

1798 

1799 Examples 

1800 -------- 

1801 Generate a noisy signal to be filtered: 

1802 

1803 >>> from scipy import signal 

1804 >>> import matplotlib.pyplot as plt 

1805 >>> t = np.linspace(-1, 1, 201) 

1806 >>> x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) + 

1807 ... 0.1*np.sin(2*np.pi*1.25*t + 1) + 

1808 ... 0.18*np.cos(2*np.pi*3.85*t)) 

1809 >>> xn = x + np.random.randn(len(t)) * 0.08 

1810 

1811 Create an order 3 lowpass butterworth filter: 

1812 

1813 >>> b, a = signal.butter(3, 0.05) 

1814 

1815 Apply the filter to xn. Use lfilter_zi to choose the initial condition of 

1816 the filter: 

1817 

1818 >>> zi = signal.lfilter_zi(b, a) 

1819 >>> z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0]) 

1820 

1821 Apply the filter again, to have a result filtered at an order the same as 

1822 filtfilt: 

1823 

1824 >>> z2, _ = signal.lfilter(b, a, z, zi=zi*z[0]) 

1825 

1826 Use filtfilt to apply the filter: 

1827 

1828 >>> y = signal.filtfilt(b, a, xn) 

1829 

1830 Plot the original signal and the various filtered versions: 

1831 

1832 >>> plt.figure 

1833 >>> plt.plot(t, xn, 'b', alpha=0.75) 

1834 >>> plt.plot(t, z, 'r--', t, z2, 'r', t, y, 'k') 

1835 >>> plt.legend(('noisy signal', 'lfilter, once', 'lfilter, twice', 

1836 ... 'filtfilt'), loc='best') 

1837 >>> plt.grid(True) 

1838 >>> plt.show() 

1839 

1840 """ 

1841 a = np.atleast_1d(a) 

1842 if len(a) == 1: 

1843 # This path only supports types fdgFDGO to mirror _linear_filter below. 

1844 # Any of b, a, x, or zi can set the dtype, but there is no default 

1845 # casting of other types; instead a NotImplementedError is raised. 

1846 b = np.asarray(b) 

1847 a = np.asarray(a) 

1848 if b.ndim != 1 and a.ndim != 1: 

1849 raise ValueError('object of too small depth for desired array') 

1850 x = _validate_x(x) 

1851 inputs = [b, a, x] 

1852 if zi is not None: 

1853 # _linear_filter does not broadcast zi, but does do expansion of 

1854 # singleton dims. 

1855 zi = np.asarray(zi) 

1856 if zi.ndim != x.ndim: 

1857 raise ValueError('object of too small depth for desired array') 

1858 expected_shape = list(x.shape) 

1859 expected_shape[axis] = b.shape[0] - 1 

1860 expected_shape = tuple(expected_shape) 

1861 # check the trivial case where zi is the right shape first 

1862 if zi.shape != expected_shape: 

1863 strides = zi.ndim * [None] 

1864 if axis < 0: 

1865 axis += zi.ndim 

1866 for k in range(zi.ndim): 

1867 if k == axis and zi.shape[k] == expected_shape[k]: 

1868 strides[k] = zi.strides[k] 

1869 elif k != axis and zi.shape[k] == expected_shape[k]: 

1870 strides[k] = zi.strides[k] 

1871 elif k != axis and zi.shape[k] == 1: 

1872 strides[k] = 0 

1873 else: 

1874 raise ValueError('Unexpected shape for zi: expected ' 

1875 '%s, found %s.' % 

1876 (expected_shape, zi.shape)) 

1877 zi = np.lib.stride_tricks.as_strided(zi, expected_shape, 

1878 strides) 

1879 inputs.append(zi) 

1880 dtype = np.result_type(*inputs) 

1881 

1882 if dtype.char not in 'fdgFDGO': 

1883 raise NotImplementedError("input type '%s' not supported" % dtype) 

1884 

1885 b = np.array(b, dtype=dtype) 

1886 a = np.array(a, dtype=dtype, copy=False) 

1887 b /= a[0] 

1888 x = np.array(x, dtype=dtype, copy=False) 

1889 

1890 out_full = np.apply_along_axis(lambda y: np.convolve(b, y), axis, x) 

1891 ind = out_full.ndim * [slice(None)] 

1892 if zi is not None: 

1893 ind[axis] = slice(zi.shape[axis]) 

1894 out_full[tuple(ind)] += zi 

1895 

1896 ind[axis] = slice(out_full.shape[axis] - len(b) + 1) 

1897 out = out_full[tuple(ind)] 

1898 

1899 if zi is None: 

1900 return out 

1901 else: 

1902 ind[axis] = slice(out_full.shape[axis] - len(b) + 1, None) 

1903 zf = out_full[tuple(ind)] 

1904 return out, zf 

1905 else: 

1906 if zi is None: 

1907 return sigtools._linear_filter(b, a, x, axis) 

1908 else: 

1909 return sigtools._linear_filter(b, a, x, axis, zi) 

1910 

1911 

1912def lfiltic(b, a, y, x=None): 

1913 """ 

1914 Construct initial conditions for lfilter given input and output vectors. 

1915 

1916 Given a linear filter (b, a) and initial conditions on the output `y` 

1917 and the input `x`, return the initial conditions on the state vector zi 

1918 which is used by `lfilter` to generate the output given the input. 

1919 

1920 Parameters 

1921 ---------- 

1922 b : array_like 

1923 Linear filter term. 

1924 a : array_like 

1925 Linear filter term. 

1926 y : array_like 

1927 Initial conditions. 

1928 

1929 If ``N = len(a) - 1``, then ``y = {y[-1], y[-2], ..., y[-N]}``. 

1930 

1931 If `y` is too short, it is padded with zeros. 

1932 x : array_like, optional 

1933 Initial conditions. 

1934 

1935 If ``M = len(b) - 1``, then ``x = {x[-1], x[-2], ..., x[-M]}``. 

1936 

1937 If `x` is not given, its initial conditions are assumed zero. 

1938 

1939 If `x` is too short, it is padded with zeros. 

1940 

1941 Returns 

1942 ------- 

1943 zi : ndarray 

1944 The state vector ``zi = {z_0[-1], z_1[-1], ..., z_K-1[-1]}``, 

1945 where ``K = max(M, N)``. 

1946 

1947 See Also 

1948 -------- 

1949 lfilter, lfilter_zi 

1950 

1951 """ 

1952 N = np.size(a) - 1 

1953 M = np.size(b) - 1 

1954 K = max(M, N) 

1955 y = np.asarray(y) 

1956 if y.dtype.kind in 'bui': 

1957 # ensure calculations are floating point 

1958 y = y.astype(np.float64) 

1959 zi = np.zeros(K, y.dtype) 

1960 if x is None: 

1961 x = np.zeros(M, y.dtype) 

1962 else: 

1963 x = np.asarray(x) 

1964 L = np.size(x) 

1965 if L < M: 

1966 x = np.r_[x, np.zeros(M - L)] 

1967 L = np.size(y) 

1968 if L < N: 

1969 y = np.r_[y, np.zeros(N - L)] 

1970 

1971 for m in range(M): 

1972 zi[m] = np.sum(b[m + 1:] * x[:M - m], axis=0) 

1973 

1974 for m in range(N): 

1975 zi[m] -= np.sum(a[m + 1:] * y[:N - m], axis=0) 

1976 

1977 return zi 

1978 

1979 

1980def deconvolve(signal, divisor): 

1981 """Deconvolves ``divisor`` out of ``signal`` using inverse filtering. 

1982 

1983 Returns the quotient and remainder such that 

1984 ``signal = convolve(divisor, quotient) + remainder`` 

1985 

1986 Parameters 

1987 ---------- 

1988 signal : array_like 

1989 Signal data, typically a recorded signal 

1990 divisor : array_like 

1991 Divisor data, typically an impulse response or filter that was 

1992 applied to the original signal 

1993 

1994 Returns 

1995 ------- 

1996 quotient : ndarray 

1997 Quotient, typically the recovered original signal 

1998 remainder : ndarray 

1999 Remainder 

2000 

2001 Examples 

2002 -------- 

2003 Deconvolve a signal that's been filtered: 

2004 

2005 >>> from scipy import signal 

2006 >>> original = [0, 1, 0, 0, 1, 1, 0, 0] 

2007 >>> impulse_response = [2, 1] 

2008 >>> recorded = signal.convolve(impulse_response, original) 

2009 >>> recorded 

2010 array([0, 2, 1, 0, 2, 3, 1, 0, 0]) 

2011 >>> recovered, remainder = signal.deconvolve(recorded, impulse_response) 

2012 >>> recovered 

2013 array([ 0., 1., 0., 0., 1., 1., 0., 0.]) 

2014 

2015 See Also 

2016 -------- 

2017 numpy.polydiv : performs polynomial division (same operation, but 

2018 also accepts poly1d objects) 

2019 

2020 """ 

2021 num = np.atleast_1d(signal) 

2022 den = np.atleast_1d(divisor) 

2023 N = len(num) 

2024 D = len(den) 

2025 if D > N: 

2026 quot = [] 

2027 rem = num 

2028 else: 

2029 input = np.zeros(N - D + 1, float) 

2030 input[0] = 1 

2031 quot = lfilter(num, den, input) 

2032 rem = num - convolve(den, quot, mode='full') 

2033 return quot, rem 

2034 

2035 

2036def hilbert(x, N=None, axis=-1): 

2037 """ 

2038 Compute the analytic signal, using the Hilbert transform. 

2039 

2040 The transformation is done along the last axis by default. 

2041 

2042 Parameters 

2043 ---------- 

2044 x : array_like 

2045 Signal data. Must be real. 

2046 N : int, optional 

2047 Number of Fourier components. Default: ``x.shape[axis]`` 

2048 axis : int, optional 

2049 Axis along which to do the transformation. Default: -1. 

2050 

2051 Returns 

2052 ------- 

2053 xa : ndarray 

2054 Analytic signal of `x`, of each 1-D array along `axis` 

2055 

2056 Notes 

2057 ----- 

2058 The analytic signal ``x_a(t)`` of signal ``x(t)`` is: 

2059 

2060 .. math:: x_a = F^{-1}(F(x) 2U) = x + i y 

2061 

2062 where `F` is the Fourier transform, `U` the unit step function, 

2063 and `y` the Hilbert transform of `x`. [1]_ 

2064 

2065 In other words, the negative half of the frequency spectrum is zeroed 

2066 out, turning the real-valued signal into a complex signal. The Hilbert 

2067 transformed signal can be obtained from ``np.imag(hilbert(x))``, and the 

2068 original signal from ``np.real(hilbert(x))``. 

2069 

2070 Examples 

2071 --------- 

2072 In this example we use the Hilbert transform to determine the amplitude 

2073 envelope and instantaneous frequency of an amplitude-modulated signal. 

2074 

2075 >>> import numpy as np 

2076 >>> import matplotlib.pyplot as plt 

2077 >>> from scipy.signal import hilbert, chirp 

2078 

2079 >>> duration = 1.0 

2080 >>> fs = 400.0 

2081 >>> samples = int(fs*duration) 

2082 >>> t = np.arange(samples) / fs 

2083 

2084 We create a chirp of which the frequency increases from 20 Hz to 100 Hz and 

2085 apply an amplitude modulation. 

2086 

2087 >>> signal = chirp(t, 20.0, t[-1], 100.0) 

2088 >>> signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) ) 

2089 

2090 The amplitude envelope is given by magnitude of the analytic signal. The 

2091 instantaneous frequency can be obtained by differentiating the 

2092 instantaneous phase in respect to time. The instantaneous phase corresponds 

2093 to the phase angle of the analytic signal. 

2094 

2095 >>> analytic_signal = hilbert(signal) 

2096 >>> amplitude_envelope = np.abs(analytic_signal) 

2097 >>> instantaneous_phase = np.unwrap(np.angle(analytic_signal)) 

2098 >>> instantaneous_frequency = (np.diff(instantaneous_phase) / 

2099 ... (2.0*np.pi) * fs) 

2100 

2101 >>> fig = plt.figure() 

2102 >>> ax0 = fig.add_subplot(211) 

2103 >>> ax0.plot(t, signal, label='signal') 

2104 >>> ax0.plot(t, amplitude_envelope, label='envelope') 

2105 >>> ax0.set_xlabel("time in seconds") 

2106 >>> ax0.legend() 

2107 >>> ax1 = fig.add_subplot(212) 

2108 >>> ax1.plot(t[1:], instantaneous_frequency) 

2109 >>> ax1.set_xlabel("time in seconds") 

2110 >>> ax1.set_ylim(0.0, 120.0) 

2111 

2112 References 

2113 ---------- 

2114 .. [1] Wikipedia, "Analytic signal". 

2115 https://en.wikipedia.org/wiki/Analytic_signal 

2116 .. [2] Leon Cohen, "Time-Frequency Analysis", 1995. Chapter 2. 

2117 .. [3] Alan V. Oppenheim, Ronald W. Schafer. Discrete-Time Signal 

2118 Processing, Third Edition, 2009. Chapter 12. 

2119 ISBN 13: 978-1292-02572-8 

2120 

2121 """ 

2122 x = np.asarray(x) 

2123 if np.iscomplexobj(x): 

2124 raise ValueError("x must be real.") 

2125 if N is None: 

2126 N = x.shape[axis] 

2127 if N <= 0: 

2128 raise ValueError("N must be positive.") 

2129 

2130 Xf = sp_fft.fft(x, N, axis=axis) 

2131 h = np.zeros(N) 

2132 if N % 2 == 0: 

2133 h[0] = h[N // 2] = 1 

2134 h[1:N // 2] = 2 

2135 else: 

2136 h[0] = 1 

2137 h[1:(N + 1) // 2] = 2 

2138 

2139 if x.ndim > 1: 

2140 ind = [np.newaxis] * x.ndim 

2141 ind[axis] = slice(None) 

2142 h = h[tuple(ind)] 

2143 x = sp_fft.ifft(Xf * h, axis=axis) 

2144 return x 

2145 

2146 

2147def hilbert2(x, N=None): 

2148 """ 

2149 Compute the '2-D' analytic signal of `x` 

2150 

2151 Parameters 

2152 ---------- 

2153 x : array_like 

2154 2-D signal data. 

2155 N : int or tuple of two ints, optional 

2156 Number of Fourier components. Default is ``x.shape`` 

2157 

2158 Returns 

2159 ------- 

2160 xa : ndarray 

2161 Analytic signal of `x` taken along axes (0,1). 

2162 

2163 References 

2164 ---------- 

2165 .. [1] Wikipedia, "Analytic signal", 

2166 https://en.wikipedia.org/wiki/Analytic_signal 

2167 

2168 """ 

2169 x = np.atleast_2d(x) 

2170 if x.ndim > 2: 

2171 raise ValueError("x must be 2-D.") 

2172 if np.iscomplexobj(x): 

2173 raise ValueError("x must be real.") 

2174 if N is None: 

2175 N = x.shape 

2176 elif isinstance(N, int): 

2177 if N <= 0: 

2178 raise ValueError("N must be positive.") 

2179 N = (N, N) 

2180 elif len(N) != 2 or np.any(np.asarray(N) <= 0): 

2181 raise ValueError("When given as a tuple, N must hold exactly " 

2182 "two positive integers") 

2183 

2184 Xf = sp_fft.fft2(x, N, axes=(0, 1)) 

2185 h1 = np.zeros(N[0], 'd') 

2186 h2 = np.zeros(N[1], 'd') 

2187 for p in range(2): 

2188 h = eval("h%d" % (p + 1)) 

2189 N1 = N[p] 

2190 if N1 % 2 == 0: 

2191 h[0] = h[N1 // 2] = 1 

2192 h[1:N1 // 2] = 2 

2193 else: 

2194 h[0] = 1 

2195 h[1:(N1 + 1) // 2] = 2 

2196 exec("h%d = h" % (p + 1), globals(), locals()) 

2197 

2198 h = h1[:, np.newaxis] * h2[np.newaxis, :] 

2199 k = x.ndim 

2200 while k > 2: 

2201 h = h[:, np.newaxis] 

2202 k -= 1 

2203 x = sp_fft.ifft2(Xf * h, axes=(0, 1)) 

2204 return x 

2205 

2206 

2207def cmplx_sort(p): 

2208 """Sort roots based on magnitude. 

2209 

2210 Parameters 

2211 ---------- 

2212 p : array_like 

2213 The roots to sort, as a 1-D array. 

2214 

2215 Returns 

2216 ------- 

2217 p_sorted : ndarray 

2218 Sorted roots. 

2219 indx : ndarray 

2220 Array of indices needed to sort the input `p`. 

2221 

2222 Examples 

2223 -------- 

2224 >>> from scipy import signal 

2225 >>> vals = [1, 4, 1+1.j, 3] 

2226 >>> p_sorted, indx = signal.cmplx_sort(vals) 

2227 >>> p_sorted 

2228 array([1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j]) 

2229 >>> indx 

2230 array([0, 2, 3, 1]) 

2231 """ 

2232 p = np.asarray(p) 

2233 indx = np.argsort(abs(p)) 

2234 return np.take(p, indx, 0), indx 

2235 

2236 

2237def unique_roots(p, tol=1e-3, rtype='min'): 

2238 """Determine unique roots and their multiplicities from a list of roots. 

2239 

2240 Parameters 

2241 ---------- 

2242 p : array_like 

2243 The list of roots. 

2244 tol : float, optional 

2245 The tolerance for two roots to be considered equal in terms of 

2246 the distance between them. Default is 1e-3. Refer to Notes about 

2247 the details on roots grouping. 

2248 rtype : {'max', 'maximum', 'min', 'minimum', 'avg', 'mean'}, optional 

2249 How to determine the returned root if multiple roots are within 

2250 `tol` of each other. 

2251 

2252 - 'max', 'maximum': pick the maximum of those roots 

2253 - 'min', 'minimum': pick the minimum of those roots 

2254 - 'avg', 'mean': take the average of those roots 

2255 

2256 When finding minimum or maximum among complex roots they are compared 

2257 first by the real part and then by the imaginary part. 

2258 

2259 Returns 

2260 ------- 

2261 unique : ndarray 

2262 The list of unique roots. 

2263 multiplicity : ndarray 

2264 The multiplicity of each root. 

2265 

2266 Notes 

2267 ----- 

2268 If we have 3 roots ``a``, ``b`` and ``c``, such that ``a`` is close to 

2269 ``b`` and ``b`` is close to ``c`` (distance is less than `tol`), then it 

2270 doesn't necessarily mean that ``a`` is close to ``c``. It means that roots 

2271 grouping is not unique. In this function we use "greedy" grouping going 

2272 through the roots in the order they are given in the input `p`. 

2273 

2274 This utility function is not specific to roots but can be used for any 

2275 sequence of values for which uniqueness and multiplicity has to be 

2276 determined. For a more general routine, see `numpy.unique`. 

2277 

2278 Examples 

2279 -------- 

2280 >>> from scipy import signal 

2281 >>> vals = [0, 1.3, 1.31, 2.8, 1.25, 2.2, 10.3] 

2282 >>> uniq, mult = signal.unique_roots(vals, tol=2e-2, rtype='avg') 

2283 

2284 Check which roots have multiplicity larger than 1: 

2285 

2286 >>> uniq[mult > 1] 

2287 array([ 1.305]) 

2288 """ 

2289 if rtype in ['max', 'maximum']: 

2290 reduce = np.max 

2291 elif rtype in ['min', 'minimum']: 

2292 reduce = np.min 

2293 elif rtype in ['avg', 'mean']: 

2294 reduce = np.mean 

2295 else: 

2296 raise ValueError("`rtype` must be one of " 

2297 "{'max', 'maximum', 'min', 'minimum', 'avg', 'mean'}") 

2298 

2299 p = np.asarray(p) 

2300 

2301 points = np.empty((len(p), 2)) 

2302 points[:, 0] = np.real(p) 

2303 points[:, 1] = np.imag(p) 

2304 tree = cKDTree(points) 

2305 

2306 p_unique = [] 

2307 p_multiplicity = [] 

2308 used = np.zeros(len(p), dtype=bool) 

2309 for i in range(len(p)): 

2310 if used[i]: 

2311 continue 

2312 

2313 group = tree.query_ball_point(points[i], tol) 

2314 group = [x for x in group if not used[x]] 

2315 

2316 p_unique.append(reduce(p[group])) 

2317 p_multiplicity.append(len(group)) 

2318 

2319 used[group] = True 

2320 

2321 return np.asarray(p_unique), np.asarray(p_multiplicity) 

2322 

2323 

2324def invres(r, p, k, tol=1e-3, rtype='avg'): 

2325 """Compute b(s) and a(s) from partial fraction expansion. 

2326 

2327 If `M` is the degree of numerator `b` and `N` the degree of denominator 

2328 `a`:: 

2329 

2330 b(s) b[0] s**(M) + b[1] s**(M-1) + ... + b[M] 

2331 H(s) = ------ = ------------------------------------------ 

2332 a(s) a[0] s**(N) + a[1] s**(N-1) + ... + a[N] 

2333 

2334 then the partial-fraction expansion H(s) is defined as:: 

2335 

2336 r[0] r[1] r[-1] 

2337 = -------- + -------- + ... + --------- + k(s) 

2338 (s-p[0]) (s-p[1]) (s-p[-1]) 

2339 

2340 If there are any repeated roots (closer together than `tol`), then H(s) 

2341 has terms like:: 

2342 

2343 r[i] r[i+1] r[i+n-1] 

2344 -------- + ----------- + ... + ----------- 

2345 (s-p[i]) (s-p[i])**2 (s-p[i])**n 

2346 

2347 This function is used for polynomials in positive powers of s or z, 

2348 such as analog filters or digital filters in controls engineering. For 

2349 negative powers of z (typical for digital filters in DSP), use `invresz`. 

2350 

2351 Parameters 

2352 ---------- 

2353 r : array_like 

2354 Residues corresponding to the poles. For repeated poles, the residues 

2355 must be ordered to correspond to ascending by power fractions. 

2356 p : array_like 

2357 Poles. Equal poles must be adjacent. 

2358 k : array_like 

2359 Coefficients of the direct polynomial term. 

2360 tol : float, optional 

2361 The tolerance for two roots to be considered equal in terms of 

2362 the distance between them. Default is 1e-3. See `unique_roots` 

2363 for further details. 

2364 rtype : {'avg', 'min', 'max'}, optional 

2365 Method for computing a root to represent a group of identical roots. 

2366 Default is 'avg'. See `unique_roots` for further details. 

2367 

2368 Returns 

2369 ------- 

2370 b : ndarray 

2371 Numerator polynomial coefficients. 

2372 a : ndarray 

2373 Denominator polynomial coefficients. 

2374 

2375 See Also 

2376 -------- 

2377 residue, invresz, unique_roots 

2378 

2379 """ 

2380 r = np.atleast_1d(r) 

2381 p = np.atleast_1d(p) 

2382 k = np.trim_zeros(np.atleast_1d(k), 'f') 

2383 

2384 unique_poles, multiplicity = _group_poles(p, tol, rtype) 

2385 factors, denominator = _compute_factors(unique_poles, multiplicity, 

2386 include_powers=True) 

2387 

2388 if len(k) == 0: 

2389 numerator = 0 

2390 else: 

2391 numerator = np.polymul(k, denominator) 

2392 

2393 for residue, factor in zip(r, factors): 

2394 numerator = np.polyadd(numerator, residue * factor) 

2395 

2396 return numerator, denominator 

2397 

2398 

2399def _compute_factors(roots, multiplicity, include_powers=False): 

2400 """Compute the total polynomial divided by factors for each root.""" 

2401 current = np.array([1]) 

2402 suffixes = [current] 

2403 for pole, mult in zip(roots[-1:0:-1], multiplicity[-1:0:-1]): 

2404 monomial = np.array([1, -pole]) 

2405 for _ in range(mult): 

2406 current = np.polymul(current, monomial) 

2407 suffixes.append(current) 

2408 suffixes = suffixes[::-1] 

2409 

2410 factors = [] 

2411 current = np.array([1]) 

2412 for pole, mult, suffix in zip(roots, multiplicity, suffixes): 

2413 monomial = np.array([1, -pole]) 

2414 block = [] 

2415 for i in range(mult): 

2416 if i == 0 or include_powers: 

2417 block.append(np.polymul(current, suffix)) 

2418 current = np.polymul(current, monomial) 

2419 factors.extend(reversed(block)) 

2420 

2421 return factors, current 

2422 

2423 

2424def _compute_residues(poles, multiplicity, numerator): 

2425 denominator_factors, _ = _compute_factors(poles, multiplicity) 

2426 numerator = numerator.astype(poles.dtype) 

2427 

2428 residues = [] 

2429 for pole, mult, factor in zip(poles, multiplicity, 

2430 denominator_factors): 

2431 if mult == 1: 

2432 residues.append(np.polyval(numerator, pole) / 

2433 np.polyval(factor, pole)) 

2434 else: 

2435 numer = numerator.copy() 

2436 monomial = np.array([1, -pole]) 

2437 factor, d = np.polydiv(factor, monomial) 

2438 

2439 block = [] 

2440 for _ in range(mult): 

2441 numer, n = np.polydiv(numer, monomial) 

2442 r = n[0] / d[0] 

2443 numer = np.polysub(numer, r * factor) 

2444 block.append(r) 

2445 

2446 residues.extend(reversed(block)) 

2447 

2448 return np.asarray(residues) 

2449 

2450 

2451def residue(b, a, tol=1e-3, rtype='avg'): 

2452 """Compute partial-fraction expansion of b(s) / a(s). 

2453 

2454 If `M` is the degree of numerator `b` and `N` the degree of denominator 

2455 `a`:: 

2456 

2457 b(s) b[0] s**(M) + b[1] s**(M-1) + ... + b[M] 

2458 H(s) = ------ = ------------------------------------------ 

2459 a(s) a[0] s**(N) + a[1] s**(N-1) + ... + a[N] 

2460 

2461 then the partial-fraction expansion H(s) is defined as:: 

2462 

2463 r[0] r[1] r[-1] 

2464 = -------- + -------- + ... + --------- + k(s) 

2465 (s-p[0]) (s-p[1]) (s-p[-1]) 

2466 

2467 If there are any repeated roots (closer together than `tol`), then H(s) 

2468 has terms like:: 

2469 

2470 r[i] r[i+1] r[i+n-1] 

2471 -------- + ----------- + ... + ----------- 

2472 (s-p[i]) (s-p[i])**2 (s-p[i])**n 

2473 

2474 This function is used for polynomials in positive powers of s or z, 

2475 such as analog filters or digital filters in controls engineering. For 

2476 negative powers of z (typical for digital filters in DSP), use `residuez`. 

2477 

2478 See Notes for details about the algorithm. 

2479 

2480 Parameters 

2481 ---------- 

2482 b : array_like 

2483 Numerator polynomial coefficients. 

2484 a : array_like 

2485 Denominator polynomial coefficients. 

2486 tol : float, optional 

2487 The tolerance for two roots to be considered equal in terms of 

2488 the distance between them. Default is 1e-3. See `unique_roots` 

2489 for further details. 

2490 rtype : {'avg', 'min', 'max'}, optional 

2491 Method for computing a root to represent a group of identical roots. 

2492 Default is 'avg'. See `unique_roots` for further details. 

2493 

2494 Returns 

2495 ------- 

2496 r : ndarray 

2497 Residues corresponding to the poles. For repeated poles, the residues 

2498 are ordered to correspond to ascending by power fractions. 

2499 p : ndarray 

2500 Poles ordered by magnitude in ascending order. 

2501 k : ndarray 

2502 Coefficients of the direct polynomial term. 

2503 

2504 See Also 

2505 -------- 

2506 invres, residuez, numpy.poly, unique_roots 

2507 

2508 Notes 

2509 ----- 

2510 The "deflation through subtraction" algorithm is used for 

2511 computations --- method 6 in [1]_. 

2512 

2513 The form of partial fraction expansion depends on poles multiplicity in 

2514 the exact mathematical sense. However there is no way to exactly 

2515 determine multiplicity of roots of a polynomial in numerical computing. 

2516 Thus you should think of the result of `residue` with given `tol` as 

2517 partial fraction expansion computed for the denominator composed of the 

2518 computed poles with empirically determined multiplicity. The choice of 

2519 `tol` can drastically change the result if there are close poles. 

2520 

2521 References 

2522 ---------- 

2523 .. [1] J. F. Mahoney, B. D. Sivazlian, "Partial fractions expansion: a 

2524 review of computational methodology and efficiency", Journal of 

2525 Computational and Applied Mathematics, Vol. 9, 1983. 

2526 """ 

2527 b = np.asarray(b) 

2528 a = np.asarray(a) 

2529 if (np.issubdtype(b.dtype, np.complexfloating) 

2530 or np.issubdtype(a.dtype, np.complexfloating)): 

2531 b = b.astype(complex) 

2532 a = a.astype(complex) 

2533 else: 

2534 b = b.astype(float) 

2535 a = a.astype(float) 

2536 

2537 b = np.trim_zeros(np.atleast_1d(b), 'f') 

2538 a = np.trim_zeros(np.atleast_1d(a), 'f') 

2539 

2540 if a.size == 0: 

2541 raise ValueError("Denominator `a` is zero.") 

2542 

2543 poles = np.roots(a) 

2544 if b.size == 0: 

2545 return np.zeros(poles.shape), cmplx_sort(poles)[0], np.array([]) 

2546 

2547 if len(b) < len(a): 

2548 k = np.empty(0) 

2549 else: 

2550 k, b = np.polydiv(b, a) 

2551 

2552 unique_poles, multiplicity = unique_roots(poles, tol=tol, rtype=rtype) 

2553 unique_poles, order = cmplx_sort(unique_poles) 

2554 multiplicity = multiplicity[order] 

2555 

2556 residues = _compute_residues(unique_poles, multiplicity, b) 

2557 

2558 index = 0 

2559 for pole, mult in zip(unique_poles, multiplicity): 

2560 poles[index:index + mult] = pole 

2561 index += mult 

2562 

2563 return residues / a[0], poles, k 

2564 

2565 

2566def residuez(b, a, tol=1e-3, rtype='avg'): 

2567 """Compute partial-fraction expansion of b(z) / a(z). 

2568 

2569 If `M` is the degree of numerator `b` and `N` the degree of denominator 

2570 `a`:: 

2571 

2572 b(z) b[0] + b[1] z**(-1) + ... + b[M] z**(-M) 

2573 H(z) = ------ = ------------------------------------------ 

2574 a(z) a[0] + a[1] z**(-1) + ... + a[N] z**(-N) 

2575 

2576 then the partial-fraction expansion H(z) is defined as:: 

2577 

2578 r[0] r[-1] 

2579 = --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ... 

2580 (1-p[0]z**(-1)) (1-p[-1]z**(-1)) 

2581 

2582 If there are any repeated roots (closer than `tol`), then the partial 

2583 fraction expansion has terms like:: 

2584 

2585 r[i] r[i+1] r[i+n-1] 

2586 -------------- + ------------------ + ... + ------------------ 

2587 (1-p[i]z**(-1)) (1-p[i]z**(-1))**2 (1-p[i]z**(-1))**n 

2588 

2589 This function is used for polynomials in negative powers of z, 

2590 such as digital filters in DSP. For positive powers, use `residue`. 

2591 

2592 See Notes of `residue` for details about the algorithm. 

2593 

2594 Parameters 

2595 ---------- 

2596 b : array_like 

2597 Numerator polynomial coefficients. 

2598 a : array_like 

2599 Denominator polynomial coefficients. 

2600 tol : float, optional 

2601 The tolerance for two roots to be considered equal in terms of 

2602 the distance between them. Default is 1e-3. See `unique_roots` 

2603 for further details. 

2604 rtype : {'avg', 'min', 'max'}, optional 

2605 Method for computing a root to represent a group of identical roots. 

2606 Default is 'avg'. See `unique_roots` for further details. 

2607 

2608 Returns 

2609 ------- 

2610 r : ndarray 

2611 Residues corresponding to the poles. For repeated poles, the residues 

2612 are ordered to correspond to ascending by power fractions. 

2613 p : ndarray 

2614 Poles ordered by magnitude in ascending order. 

2615 k : ndarray 

2616 Coefficients of the direct polynomial term. 

2617 

2618 See Also 

2619 -------- 

2620 invresz, residue, unique_roots 

2621 """ 

2622 b = np.asarray(b) 

2623 a = np.asarray(a) 

2624 if (np.issubdtype(b.dtype, np.complexfloating) 

2625 or np.issubdtype(a.dtype, np.complexfloating)): 

2626 b = b.astype(complex) 

2627 a = a.astype(complex) 

2628 else: 

2629 b = b.astype(float) 

2630 a = a.astype(float) 

2631 

2632 b = np.trim_zeros(np.atleast_1d(b), 'b') 

2633 a = np.trim_zeros(np.atleast_1d(a), 'b') 

2634 

2635 if a.size == 0: 

2636 raise ValueError("Denominator `a` is zero.") 

2637 elif a[0] == 0: 

2638 raise ValueError("First coefficient of determinant `a` must be " 

2639 "non-zero.") 

2640 

2641 poles = np.roots(a) 

2642 if b.size == 0: 

2643 return np.zeros(poles.shape), cmplx_sort(poles)[0], np.array([]) 

2644 

2645 b_rev = b[::-1] 

2646 a_rev = a[::-1] 

2647 

2648 if len(b_rev) < len(a_rev): 

2649 k_rev = np.empty(0) 

2650 else: 

2651 k_rev, b_rev = np.polydiv(b_rev, a_rev) 

2652 

2653 unique_poles, multiplicity = unique_roots(poles, tol=tol, rtype=rtype) 

2654 unique_poles, order = cmplx_sort(unique_poles) 

2655 multiplicity = multiplicity[order] 

2656 

2657 residues = _compute_residues(1 / unique_poles, multiplicity, b_rev) 

2658 

2659 index = 0 

2660 powers = np.empty(len(residues), dtype=int) 

2661 for pole, mult in zip(unique_poles, multiplicity): 

2662 poles[index:index + mult] = pole 

2663 powers[index:index + mult] = 1 + np.arange(mult) 

2664 index += mult 

2665 

2666 residues *= (-poles) ** powers / a_rev[0] 

2667 

2668 return residues, poles, k_rev[::-1] 

2669 

2670 

2671def _group_poles(poles, tol, rtype): 

2672 if rtype in ['max', 'maximum']: 

2673 reduce = np.max 

2674 elif rtype in ['min', 'minimum']: 

2675 reduce = np.min 

2676 elif rtype in ['avg', 'mean']: 

2677 reduce = np.mean 

2678 else: 

2679 raise ValueError("`rtype` must be one of " 

2680 "{'max', 'maximum', 'min', 'minimum', 'avg', 'mean'}") 

2681 

2682 unique = [] 

2683 multiplicity = [] 

2684 

2685 pole = poles[0] 

2686 block = [pole] 

2687 for i in range(1, len(poles)): 

2688 if abs(poles[i] - pole) <= tol: 

2689 block.append(pole) 

2690 else: 

2691 unique.append(reduce(block)) 

2692 multiplicity.append(len(block)) 

2693 pole = poles[i] 

2694 block = [pole] 

2695 

2696 unique.append(reduce(block)) 

2697 multiplicity.append(len(block)) 

2698 

2699 return np.asarray(unique), np.asarray(multiplicity) 

2700 

2701 

2702def invresz(r, p, k, tol=1e-3, rtype='avg'): 

2703 """Compute b(z) and a(z) from partial fraction expansion. 

2704 

2705 If `M` is the degree of numerator `b` and `N` the degree of denominator 

2706 `a`:: 

2707 

2708 b(z) b[0] + b[1] z**(-1) + ... + b[M] z**(-M) 

2709 H(z) = ------ = ------------------------------------------ 

2710 a(z) a[0] + a[1] z**(-1) + ... + a[N] z**(-N) 

2711 

2712 then the partial-fraction expansion H(z) is defined as:: 

2713 

2714 r[0] r[-1] 

2715 = --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ... 

2716 (1-p[0]z**(-1)) (1-p[-1]z**(-1)) 

2717 

2718 If there are any repeated roots (closer than `tol`), then the partial 

2719 fraction expansion has terms like:: 

2720 

2721 r[i] r[i+1] r[i+n-1] 

2722 -------------- + ------------------ + ... + ------------------ 

2723 (1-p[i]z**(-1)) (1-p[i]z**(-1))**2 (1-p[i]z**(-1))**n 

2724 

2725 This function is used for polynomials in negative powers of z, 

2726 such as digital filters in DSP. For positive powers, use `invres`. 

2727 

2728 Parameters 

2729 ---------- 

2730 r : array_like 

2731 Residues corresponding to the poles. For repeated poles, the residues 

2732 must be ordered to correspond to ascending by power fractions. 

2733 p : array_like 

2734 Poles. Equal poles must be adjacent. 

2735 k : array_like 

2736 Coefficients of the direct polynomial term. 

2737 tol : float, optional 

2738 The tolerance for two roots to be considered equal in terms of 

2739 the distance between them. Default is 1e-3. See `unique_roots` 

2740 for further details. 

2741 rtype : {'avg', 'min', 'max'}, optional 

2742 Method for computing a root to represent a group of identical roots. 

2743 Default is 'avg'. See `unique_roots` for further details. 

2744 

2745 Returns 

2746 ------- 

2747 b : ndarray 

2748 Numerator polynomial coefficients. 

2749 a : ndarray 

2750 Denominator polynomial coefficients. 

2751 

2752 See Also 

2753 -------- 

2754 residuez, unique_roots, invres 

2755 

2756 """ 

2757 r = np.atleast_1d(r) 

2758 p = np.atleast_1d(p) 

2759 k = np.trim_zeros(np.atleast_1d(k), 'b') 

2760 

2761 unique_poles, multiplicity = _group_poles(p, tol, rtype) 

2762 factors, denominator = _compute_factors(unique_poles, multiplicity, 

2763 include_powers=True) 

2764 

2765 if len(k) == 0: 

2766 numerator = 0 

2767 else: 

2768 numerator = np.polymul(k[::-1], denominator[::-1]) 

2769 

2770 for residue, factor in zip(r, factors): 

2771 numerator = np.polyadd(numerator, residue * factor[::-1]) 

2772 

2773 return numerator[::-1], denominator 

2774 

2775 

2776def resample(x, num, t=None, axis=0, window=None, domain='time'): 

2777 """ 

2778 Resample `x` to `num` samples using Fourier method along the given axis. 

2779 

2780 The resampled signal starts at the same value as `x` but is sampled 

2781 with a spacing of ``len(x) / num * (spacing of x)``. Because a 

2782 Fourier method is used, the signal is assumed to be periodic. 

2783 

2784 Parameters 

2785 ---------- 

2786 x : array_like 

2787 The data to be resampled. 

2788 num : int 

2789 The number of samples in the resampled signal. 

2790 t : array_like, optional 

2791 If `t` is given, it is assumed to be the equally spaced sample 

2792 positions associated with the signal data in `x`. 

2793 axis : int, optional 

2794 The axis of `x` that is resampled. Default is 0. 

2795 window : array_like, callable, string, float, or tuple, optional 

2796 Specifies the window applied to the signal in the Fourier 

2797 domain. See below for details. 

2798 domain : string, optional 

2799 A string indicating the domain of the input `x`: 

2800 ``time`` Consider the input `x` as time-domain (Default), 

2801 ``freq`` Consider the input `x` as frequency-domain. 

2802 

2803 Returns 

2804 ------- 

2805 resampled_x or (resampled_x, resampled_t) 

2806 Either the resampled array, or, if `t` was given, a tuple 

2807 containing the resampled array and the corresponding resampled 

2808 positions. 

2809 

2810 See Also 

2811 -------- 

2812 decimate : Downsample the signal after applying an FIR or IIR filter. 

2813 resample_poly : Resample using polyphase filtering and an FIR filter. 

2814 

2815 Notes 

2816 ----- 

2817 The argument `window` controls a Fourier-domain window that tapers 

2818 the Fourier spectrum before zero-padding to alleviate ringing in 

2819 the resampled values for sampled signals you didn't intend to be 

2820 interpreted as band-limited. 

2821 

2822 If `window` is a function, then it is called with a vector of inputs 

2823 indicating the frequency bins (i.e. fftfreq(x.shape[axis]) ). 

2824 

2825 If `window` is an array of the same length as `x.shape[axis]` it is 

2826 assumed to be the window to be applied directly in the Fourier 

2827 domain (with dc and low-frequency first). 

2828 

2829 For any other type of `window`, the function `scipy.signal.get_window` 

2830 is called to generate the window. 

2831 

2832 The first sample of the returned vector is the same as the first 

2833 sample of the input vector. The spacing between samples is changed 

2834 from ``dx`` to ``dx * len(x) / num``. 

2835 

2836 If `t` is not None, then it is used solely to calculate the resampled 

2837 positions `resampled_t` 

2838 

2839 As noted, `resample` uses FFT transformations, which can be very 

2840 slow if the number of input or output samples is large and prime; 

2841 see `scipy.fft.fft`. 

2842 

2843 Examples 

2844 -------- 

2845 Note that the end of the resampled data rises to meet the first 

2846 sample of the next cycle: 

2847 

2848 >>> from scipy import signal 

2849 

2850 >>> x = np.linspace(0, 10, 20, endpoint=False) 

2851 >>> y = np.cos(-x**2/6.0) 

2852 >>> f = signal.resample(y, 100) 

2853 >>> xnew = np.linspace(0, 10, 100, endpoint=False) 

2854 

2855 >>> import matplotlib.pyplot as plt 

2856 >>> plt.plot(x, y, 'go-', xnew, f, '.-', 10, y[0], 'ro') 

2857 >>> plt.legend(['data', 'resampled'], loc='best') 

2858 >>> plt.show() 

2859 """ 

2860 

2861 if domain not in ('time', 'freq'): 

2862 raise ValueError("Acceptable domain flags are 'time' or" 

2863 " 'freq', not domain={}".format(domain)) 

2864 

2865 x = np.asarray(x) 

2866 Nx = x.shape[axis] 

2867 

2868 # Check if we can use faster real FFT 

2869 real_input = np.isrealobj(x) 

2870 

2871 if domain == 'time': 

2872 # Forward transform 

2873 if real_input: 

2874 X = sp_fft.rfft(x, axis=axis) 

2875 else: # Full complex FFT 

2876 X = sp_fft.fft(x, axis=axis) 

2877 else: # domain == 'freq' 

2878 X = x 

2879 

2880 # Apply window to spectrum 

2881 if window is not None: 

2882 if callable(window): 

2883 W = window(sp_fft.fftfreq(Nx)) 

2884 elif isinstance(window, np.ndarray): 

2885 if window.shape != (Nx,): 

2886 raise ValueError('window must have the same length as data') 

2887 W = window 

2888 else: 

2889 W = sp_fft.ifftshift(get_window(window, Nx)) 

2890 

2891 newshape_W = [1] * x.ndim 

2892 newshape_W[axis] = X.shape[axis] 

2893 if real_input: 

2894 # Fold the window back on itself to mimic complex behavior 

2895 W_real = W.copy() 

2896 W_real[1:] += W_real[-1:0:-1] 

2897 W_real[1:] *= 0.5 

2898 X *= W_real[:newshape_W[axis]].reshape(newshape_W) 

2899 else: 

2900 X *= W.reshape(newshape_W) 

2901 

2902 # Copy each half of the original spectrum to the output spectrum, either 

2903 # truncating high frequences (downsampling) or zero-padding them 

2904 # (upsampling) 

2905 

2906 # Placeholder array for output spectrum 

2907 newshape = list(x.shape) 

2908 if real_input: 

2909 newshape[axis] = num // 2 + 1 

2910 else: 

2911 newshape[axis] = num 

2912 Y = np.zeros(newshape, X.dtype) 

2913 

2914 # Copy positive frequency components (and Nyquist, if present) 

2915 N = min(num, Nx) 

2916 nyq = N // 2 + 1 # Slice index that includes Nyquist if present 

2917 sl = [slice(None)] * x.ndim 

2918 sl[axis] = slice(0, nyq) 

2919 Y[tuple(sl)] = X[tuple(sl)] 

2920 if not real_input: 

2921 # Copy negative frequency components 

2922 if N > 2: # (slice expression doesn't collapse to empty array) 

2923 sl[axis] = slice(nyq - N, None) 

2924 Y[tuple(sl)] = X[tuple(sl)] 

2925 

2926 # Split/join Nyquist component(s) if present 

2927 # So far we have set Y[+N/2]=X[+N/2] 

2928 if N % 2 == 0: 

2929 if num < Nx: # downsampling 

2930 if real_input: 

2931 sl[axis] = slice(N//2, N//2 + 1) 

2932 Y[tuple(sl)] *= 2. 

2933 else: 

2934 # select the component of Y at frequency +N/2, 

2935 # add the component of X at -N/2 

2936 sl[axis] = slice(-N//2, -N//2 + 1) 

2937 Y[tuple(sl)] += X[tuple(sl)] 

2938 elif Nx < num: # upsampling 

2939 # select the component at frequency +N/2 and halve it 

2940 sl[axis] = slice(N//2, N//2 + 1) 

2941 Y[tuple(sl)] *= 0.5 

2942 if not real_input: 

2943 temp = Y[tuple(sl)] 

2944 # set the component at -N/2 equal to the component at +N/2 

2945 sl[axis] = slice(num-N//2, num-N//2 + 1) 

2946 Y[tuple(sl)] = temp 

2947 

2948 # Inverse transform 

2949 if real_input: 

2950 y = sp_fft.irfft(Y, num, axis=axis) 

2951 else: 

2952 y = sp_fft.ifft(Y, axis=axis, overwrite_x=True) 

2953 

2954 y *= (float(num) / float(Nx)) 

2955 

2956 if t is None: 

2957 return y 

2958 else: 

2959 new_t = np.arange(0, num) * (t[1] - t[0]) * Nx / float(num) + t[0] 

2960 return y, new_t 

2961 

2962 

2963def resample_poly(x, up, down, axis=0, window=('kaiser', 5.0), 

2964 padtype='constant', cval=None): 

2965 """ 

2966 Resample `x` along the given axis using polyphase filtering. 

2967 

2968 The signal `x` is upsampled by the factor `up`, a zero-phase low-pass 

2969 FIR filter is applied, and then it is downsampled by the factor `down`. 

2970 The resulting sample rate is ``up / down`` times the original sample 

2971 rate. By default, values beyond the boundary of the signal are assumed 

2972 to be zero during the filtering step. 

2973 

2974 Parameters 

2975 ---------- 

2976 x : array_like 

2977 The data to be resampled. 

2978 up : int 

2979 The upsampling factor. 

2980 down : int 

2981 The downsampling factor. 

2982 axis : int, optional 

2983 The axis of `x` that is resampled. Default is 0. 

2984 window : string, tuple, or array_like, optional 

2985 Desired window to use to design the low-pass filter, or the FIR filter 

2986 coefficients to employ. See below for details. 

2987 padtype : string, optional 

2988 `constant`, `line`, `mean`, `median`, `maximum`, `minimum` or any of 

2989 the other signal extension modes supported by `scipy.signal.upfirdn`. 

2990 Changes assumptions on values beyond the boundary. If `constant`, 

2991 assumed to be `cval` (default zero). If `line` assumed to continue a 

2992 linear trend defined by the first and last points. `mean`, `median`, 

2993 `maximum` and `minimum` work as in `np.pad` and assume that the values 

2994 beyond the boundary are the mean, median, maximum or minimum 

2995 respectively of the array along the axis. 

2996 

2997 .. versionadded:: 1.4.0 

2998 cval : float, optional 

2999 Value to use if `padtype='constant'`. Default is zero. 

3000 

3001 .. versionadded:: 1.4.0 

3002 

3003 Returns 

3004 ------- 

3005 resampled_x : array 

3006 The resampled array. 

3007 

3008 See Also 

3009 -------- 

3010 decimate : Downsample the signal after applying an FIR or IIR filter. 

3011 resample : Resample up or down using the FFT method. 

3012 

3013 Notes 

3014 ----- 

3015 This polyphase method will likely be faster than the Fourier method 

3016 in `scipy.signal.resample` when the number of samples is large and 

3017 prime, or when the number of samples is large and `up` and `down` 

3018 share a large greatest common denominator. The length of the FIR 

3019 filter used will depend on ``max(up, down) // gcd(up, down)``, and 

3020 the number of operations during polyphase filtering will depend on 

3021 the filter length and `down` (see `scipy.signal.upfirdn` for details). 

3022 

3023 The argument `window` specifies the FIR low-pass filter design. 

3024 

3025 If `window` is an array_like it is assumed to be the FIR filter 

3026 coefficients. Note that the FIR filter is applied after the upsampling 

3027 step, so it should be designed to operate on a signal at a sampling 

3028 frequency higher than the original by a factor of `up//gcd(up, down)`. 

3029 This function's output will be centered with respect to this array, so it 

3030 is best to pass a symmetric filter with an odd number of samples if, as 

3031 is usually the case, a zero-phase filter is desired. 

3032 

3033 For any other type of `window`, the functions `scipy.signal.get_window` 

3034 and `scipy.signal.firwin` are called to generate the appropriate filter 

3035 coefficients. 

3036 

3037 The first sample of the returned vector is the same as the first 

3038 sample of the input vector. The spacing between samples is changed 

3039 from ``dx`` to ``dx * down / float(up)``. 

3040 

3041 Examples 

3042 -------- 

3043 By default, the end of the resampled data rises to meet the first 

3044 sample of the next cycle for the FFT method, and gets closer to zero 

3045 for the polyphase method: 

3046 

3047 >>> from scipy import signal 

3048 

3049 >>> x = np.linspace(0, 10, 20, endpoint=False) 

3050 >>> y = np.cos(-x**2/6.0) 

3051 >>> f_fft = signal.resample(y, 100) 

3052 >>> f_poly = signal.resample_poly(y, 100, 20) 

3053 >>> xnew = np.linspace(0, 10, 100, endpoint=False) 

3054 

3055 >>> import matplotlib.pyplot as plt 

3056 >>> plt.plot(xnew, f_fft, 'b.-', xnew, f_poly, 'r.-') 

3057 >>> plt.plot(x, y, 'ko-') 

3058 >>> plt.plot(10, y[0], 'bo', 10, 0., 'ro') # boundaries 

3059 >>> plt.legend(['resample', 'resamp_poly', 'data'], loc='best') 

3060 >>> plt.show() 

3061 

3062 This default behaviour can be changed by using the padtype option: 

3063 

3064 >>> import numpy as np 

3065 >>> from scipy import signal 

3066 

3067 >>> N = 5 

3068 >>> x = np.linspace(0, 1, N, endpoint=False) 

3069 >>> y = 2 + x**2 - 1.7*np.sin(x) + .2*np.cos(11*x) 

3070 >>> y2 = 1 + x**3 + 0.1*np.sin(x) + .1*np.cos(11*x) 

3071 >>> Y = np.stack([y, y2], axis=-1) 

3072 >>> up = 4 

3073 >>> xr = np.linspace(0, 1, N*up, endpoint=False) 

3074 

3075 >>> y2 = signal.resample_poly(Y, up, 1, padtype='constant') 

3076 >>> y3 = signal.resample_poly(Y, up, 1, padtype='mean') 

3077 >>> y4 = signal.resample_poly(Y, up, 1, padtype='line') 

3078 

3079 >>> import matplotlib.pyplot as plt 

3080 >>> for i in [0,1]: 

3081 ... plt.figure() 

3082 ... plt.plot(xr, y4[:,i], 'g.', label='line') 

3083 ... plt.plot(xr, y3[:,i], 'y.', label='mean') 

3084 ... plt.plot(xr, y2[:,i], 'r.', label='constant') 

3085 ... plt.plot(x, Y[:,i], 'k-') 

3086 ... plt.legend() 

3087 >>> plt.show() 

3088 

3089 """ 

3090 x = np.asarray(x) 

3091 if up != int(up): 

3092 raise ValueError("up must be an integer") 

3093 if down != int(down): 

3094 raise ValueError("down must be an integer") 

3095 up = int(up) 

3096 down = int(down) 

3097 if up < 1 or down < 1: 

3098 raise ValueError('up and down must be >= 1') 

3099 if cval is not None and padtype != 'constant': 

3100 raise ValueError('cval has no effect when padtype is ', padtype) 

3101 

3102 # Determine our up and down factors 

3103 # Use a rational approximation to save computation time on really long 

3104 # signals 

3105 g_ = math.gcd(up, down) 

3106 up //= g_ 

3107 down //= g_ 

3108 if up == down == 1: 

3109 return x.copy() 

3110 n_in = x.shape[axis] 

3111 n_out = n_in * up 

3112 n_out = n_out // down + bool(n_out % down) 

3113 

3114 if isinstance(window, (list, np.ndarray)): 

3115 window = np.array(window) # use array to force a copy (we modify it) 

3116 if window.ndim > 1: 

3117 raise ValueError('window must be 1-D') 

3118 half_len = (window.size - 1) // 2 

3119 h = window 

3120 else: 

3121 # Design a linear-phase low-pass FIR filter 

3122 max_rate = max(up, down) 

3123 f_c = 1. / max_rate # cutoff of FIR filter (rel. to Nyquist) 

3124 half_len = 10 * max_rate # reasonable cutoff for our sinc-like function 

3125 h = firwin(2 * half_len + 1, f_c, window=window) 

3126 h *= up 

3127 

3128 # Zero-pad our filter to put the output samples at the center 

3129 n_pre_pad = (down - half_len % down) 

3130 n_post_pad = 0 

3131 n_pre_remove = (half_len + n_pre_pad) // down 

3132 # We should rarely need to do this given our filter lengths... 

3133 while _output_len(len(h) + n_pre_pad + n_post_pad, n_in, 

3134 up, down) < n_out + n_pre_remove: 

3135 n_post_pad += 1 

3136 h = np.concatenate((np.zeros(n_pre_pad, dtype=h.dtype), h, 

3137 np.zeros(n_post_pad, dtype=h.dtype))) 

3138 n_pre_remove_end = n_pre_remove + n_out 

3139 

3140 # Remove background depending on the padtype option 

3141 funcs = {'mean': np.mean, 'median': np.median, 

3142 'minimum': np.amin, 'maximum': np.amax} 

3143 upfirdn_kwargs = {'mode': 'constant', 'cval': 0} 

3144 if padtype in funcs: 

3145 background_values = funcs[padtype](x, axis=axis, keepdims=True) 

3146 elif padtype in _upfirdn_modes: 

3147 upfirdn_kwargs = {'mode': padtype} 

3148 if padtype == 'constant': 

3149 if cval is None: 

3150 cval = 0 

3151 upfirdn_kwargs['cval'] = cval 

3152 else: 

3153 raise ValueError( 

3154 'padtype must be one of: maximum, mean, median, minimum, ' + 

3155 ', '.join(_upfirdn_modes)) 

3156 

3157 if padtype in funcs: 

3158 x = x - background_values 

3159 

3160 # filter then remove excess 

3161 y = upfirdn(h, x, up, down, axis=axis, **upfirdn_kwargs) 

3162 keep = [slice(None), ]*x.ndim 

3163 keep[axis] = slice(n_pre_remove, n_pre_remove_end) 

3164 y_keep = y[tuple(keep)] 

3165 

3166 # Add background back 

3167 if padtype in funcs: 

3168 y_keep += background_values 

3169 

3170 return y_keep 

3171 

3172 

3173def vectorstrength(events, period): 

3174 ''' 

3175 Determine the vector strength of the events corresponding to the given 

3176 period. 

3177 

3178 The vector strength is a measure of phase synchrony, how well the 

3179 timing of the events is synchronized to a single period of a periodic 

3180 signal. 

3181 

3182 If multiple periods are used, calculate the vector strength of each. 

3183 This is called the "resonating vector strength". 

3184 

3185 Parameters 

3186 ---------- 

3187 events : 1D array_like 

3188 An array of time points containing the timing of the events. 

3189 period : float or array_like 

3190 The period of the signal that the events should synchronize to. 

3191 The period is in the same units as `events`. It can also be an array 

3192 of periods, in which case the outputs are arrays of the same length. 

3193 

3194 Returns 

3195 ------- 

3196 strength : float or 1D array 

3197 The strength of the synchronization. 1.0 is perfect synchronization 

3198 and 0.0 is no synchronization. If `period` is an array, this is also 

3199 an array with each element containing the vector strength at the 

3200 corresponding period. 

3201 phase : float or array 

3202 The phase that the events are most strongly synchronized to in radians. 

3203 If `period` is an array, this is also an array with each element 

3204 containing the phase for the corresponding period. 

3205 

3206 References 

3207 ---------- 

3208 van Hemmen, JL, Longtin, A, and Vollmayr, AN. Testing resonating vector 

3209 strength: Auditory system, electric fish, and noise. 

3210 Chaos 21, 047508 (2011); 

3211 :doi:`10.1063/1.3670512`. 

3212 van Hemmen, JL. Vector strength after Goldberg, Brown, and von Mises: 

3213 biological and mathematical perspectives. Biol Cybern. 

3214 2013 Aug;107(4):385-96. :doi:`10.1007/s00422-013-0561-7`. 

3215 van Hemmen, JL and Vollmayr, AN. Resonating vector strength: what happens 

3216 when we vary the "probing" frequency while keeping the spike times 

3217 fixed. Biol Cybern. 2013 Aug;107(4):491-94. 

3218 :doi:`10.1007/s00422-013-0560-8`. 

3219 ''' 

3220 events = np.asarray(events) 

3221 period = np.asarray(period) 

3222 if events.ndim > 1: 

3223 raise ValueError('events cannot have dimensions more than 1') 

3224 if period.ndim > 1: 

3225 raise ValueError('period cannot have dimensions more than 1') 

3226 

3227 # we need to know later if period was originally a scalar 

3228 scalarperiod = not period.ndim 

3229 

3230 events = np.atleast_2d(events) 

3231 period = np.atleast_2d(period) 

3232 if (period <= 0).any(): 

3233 raise ValueError('periods must be positive') 

3234 

3235 # this converts the times to vectors 

3236 vectors = np.exp(np.dot(2j*np.pi/period.T, events)) 

3237 

3238 # the vector strength is just the magnitude of the mean of the vectors 

3239 # the vector phase is the angle of the mean of the vectors 

3240 vectormean = np.mean(vectors, axis=1) 

3241 strength = abs(vectormean) 

3242 phase = np.angle(vectormean) 

3243 

3244 # if the original period was a scalar, return scalars 

3245 if scalarperiod: 

3246 strength = strength[0] 

3247 phase = phase[0] 

3248 return strength, phase 

3249 

3250 

3251def detrend(data, axis=-1, type='linear', bp=0, overwrite_data=False): 

3252 """ 

3253 Remove linear trend along axis from data. 

3254 

3255 Parameters 

3256 ---------- 

3257 data : array_like 

3258 The input data. 

3259 axis : int, optional 

3260 The axis along which to detrend the data. By default this is the 

3261 last axis (-1). 

3262 type : {'linear', 'constant'}, optional 

3263 The type of detrending. If ``type == 'linear'`` (default), 

3264 the result of a linear least-squares fit to `data` is subtracted 

3265 from `data`. 

3266 If ``type == 'constant'``, only the mean of `data` is subtracted. 

3267 bp : array_like of ints, optional 

3268 A sequence of break points. If given, an individual linear fit is 

3269 performed for each part of `data` between two break points. 

3270 Break points are specified as indices into `data`. This parameter 

3271 only has an effect when ``type == 'linear'``. 

3272 overwrite_data : bool, optional 

3273 If True, perform in place detrending and avoid a copy. Default is False 

3274 

3275 Returns 

3276 ------- 

3277 ret : ndarray 

3278 The detrended input data. 

3279 

3280 Examples 

3281 -------- 

3282 >>> from scipy import signal 

3283 >>> randgen = np.random.RandomState(9) 

3284 >>> npoints = 1000 

3285 >>> noise = randgen.randn(npoints) 

3286 >>> x = 3 + 2*np.linspace(0, 1, npoints) + noise 

3287 >>> (signal.detrend(x) - noise).max() < 0.01 

3288 True 

3289 

3290 """ 

3291 if type not in ['linear', 'l', 'constant', 'c']: 

3292 raise ValueError("Trend type must be 'linear' or 'constant'.") 

3293 data = np.asarray(data) 

3294 dtype = data.dtype.char 

3295 if dtype not in 'dfDF': 

3296 dtype = 'd' 

3297 if type in ['constant', 'c']: 

3298 ret = data - np.expand_dims(np.mean(data, axis), axis) 

3299 return ret 

3300 else: 

3301 dshape = data.shape 

3302 N = dshape[axis] 

3303 bp = np.sort(np.unique(np.r_[0, bp, N])) 

3304 if np.any(bp > N): 

3305 raise ValueError("Breakpoints must be less than length " 

3306 "of data along given axis.") 

3307 Nreg = len(bp) - 1 

3308 # Restructure data so that axis is along first dimension and 

3309 # all other dimensions are collapsed into second dimension 

3310 rnk = len(dshape) 

3311 if axis < 0: 

3312 axis = axis + rnk 

3313 newdims = np.r_[axis, 0:axis, axis + 1:rnk] 

3314 newdata = np.reshape(np.transpose(data, tuple(newdims)), 

3315 (N, _prod(dshape) // N)) 

3316 if not overwrite_data: 

3317 newdata = newdata.copy() # make sure we have a copy 

3318 if newdata.dtype.char not in 'dfDF': 

3319 newdata = newdata.astype(dtype) 

3320 # Find leastsq fit and remove it for each piece 

3321 for m in range(Nreg): 

3322 Npts = bp[m + 1] - bp[m] 

3323 A = np.ones((Npts, 2), dtype) 

3324 A[:, 0] = np.cast[dtype](np.arange(1, Npts + 1) * 1.0 / Npts) 

3325 sl = slice(bp[m], bp[m + 1]) 

3326 coef, resids, rank, s = linalg.lstsq(A, newdata[sl]) 

3327 newdata[sl] = newdata[sl] - np.dot(A, coef) 

3328 # Put data back in original shape. 

3329 tdshape = np.take(dshape, newdims, 0) 

3330 ret = np.reshape(newdata, tuple(tdshape)) 

3331 vals = list(range(1, rnk)) 

3332 olddims = vals[:axis] + [0] + vals[axis:] 

3333 ret = np.transpose(ret, tuple(olddims)) 

3334 return ret 

3335 

3336 

3337def lfilter_zi(b, a): 

3338 """ 

3339 Construct initial conditions for lfilter for step response steady-state. 

3340 

3341 Compute an initial state `zi` for the `lfilter` function that corresponds 

3342 to the steady state of the step response. 

3343 

3344 A typical use of this function is to set the initial state so that the 

3345 output of the filter starts at the same value as the first element of 

3346 the signal to be filtered. 

3347 

3348 Parameters 

3349 ---------- 

3350 b, a : array_like (1-D) 

3351 The IIR filter coefficients. See `lfilter` for more 

3352 information. 

3353 

3354 Returns 

3355 ------- 

3356 zi : 1-D ndarray 

3357 The initial state for the filter. 

3358 

3359 See Also 

3360 -------- 

3361 lfilter, lfiltic, filtfilt 

3362 

3363 Notes 

3364 ----- 

3365 A linear filter with order m has a state space representation (A, B, C, D), 

3366 for which the output y of the filter can be expressed as:: 

3367 

3368 z(n+1) = A*z(n) + B*x(n) 

3369 y(n) = C*z(n) + D*x(n) 

3370 

3371 where z(n) is a vector of length m, A has shape (m, m), B has shape 

3372 (m, 1), C has shape (1, m) and D has shape (1, 1) (assuming x(n) is 

3373 a scalar). lfilter_zi solves:: 

3374 

3375 zi = A*zi + B 

3376 

3377 In other words, it finds the initial condition for which the response 

3378 to an input of all ones is a constant. 

3379 

3380 Given the filter coefficients `a` and `b`, the state space matrices 

3381 for the transposed direct form II implementation of the linear filter, 

3382 which is the implementation used by scipy.signal.lfilter, are:: 

3383 

3384 A = scipy.linalg.companion(a).T 

3385 B = b[1:] - a[1:]*b[0] 

3386 

3387 assuming `a[0]` is 1.0; if `a[0]` is not 1, `a` and `b` are first 

3388 divided by a[0]. 

3389 

3390 Examples 

3391 -------- 

3392 The following code creates a lowpass Butterworth filter. Then it 

3393 applies that filter to an array whose values are all 1.0; the 

3394 output is also all 1.0, as expected for a lowpass filter. If the 

3395 `zi` argument of `lfilter` had not been given, the output would have 

3396 shown the transient signal. 

3397 

3398 >>> from numpy import array, ones 

3399 >>> from scipy.signal import lfilter, lfilter_zi, butter 

3400 >>> b, a = butter(5, 0.25) 

3401 >>> zi = lfilter_zi(b, a) 

3402 >>> y, zo = lfilter(b, a, ones(10), zi=zi) 

3403 >>> y 

3404 array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]) 

3405 

3406 Another example: 

3407 

3408 >>> x = array([0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0]) 

3409 >>> y, zf = lfilter(b, a, x, zi=zi*x[0]) 

3410 >>> y 

3411 array([ 0.5 , 0.5 , 0.5 , 0.49836039, 0.48610528, 

3412 0.44399389, 0.35505241]) 

3413 

3414 Note that the `zi` argument to `lfilter` was computed using 

3415 `lfilter_zi` and scaled by `x[0]`. Then the output `y` has no 

3416 transient until the input drops from 0.5 to 0.0. 

3417 

3418 """ 

3419 

3420 # FIXME: Can this function be replaced with an appropriate 

3421 # use of lfiltic? For example, when b,a = butter(N,Wn), 

3422 # lfiltic(b, a, y=numpy.ones_like(a), x=numpy.ones_like(b)). 

3423 # 

3424 

3425 # We could use scipy.signal.normalize, but it uses warnings in 

3426 # cases where a ValueError is more appropriate, and it allows 

3427 # b to be 2D. 

3428 b = np.atleast_1d(b) 

3429 if b.ndim != 1: 

3430 raise ValueError("Numerator b must be 1-D.") 

3431 a = np.atleast_1d(a) 

3432 if a.ndim != 1: 

3433 raise ValueError("Denominator a must be 1-D.") 

3434 

3435 while len(a) > 1 and a[0] == 0.0: 

3436 a = a[1:] 

3437 if a.size < 1: 

3438 raise ValueError("There must be at least one nonzero `a` coefficient.") 

3439 

3440 if a[0] != 1.0: 

3441 # Normalize the coefficients so a[0] == 1. 

3442 b = b / a[0] 

3443 a = a / a[0] 

3444 

3445 n = max(len(a), len(b)) 

3446 

3447 # Pad a or b with zeros so they are the same length. 

3448 if len(a) < n: 

3449 a = np.r_[a, np.zeros(n - len(a))] 

3450 elif len(b) < n: 

3451 b = np.r_[b, np.zeros(n - len(b))] 

3452 

3453 IminusA = np.eye(n - 1) - linalg.companion(a).T 

3454 B = b[1:] - a[1:] * b[0] 

3455 # Solve zi = A*zi + B 

3456 zi = np.linalg.solve(IminusA, B) 

3457 

3458 # For future reference: we could also use the following 

3459 # explicit formulas to solve the linear system: 

3460 # 

3461 # zi = np.zeros(n - 1) 

3462 # zi[0] = B.sum() / IminusA[:,0].sum() 

3463 # asum = 1.0 

3464 # csum = 0.0 

3465 # for k in range(1,n-1): 

3466 # asum += a[k] 

3467 # csum += b[k] - a[k]*b[0] 

3468 # zi[k] = asum*zi[0] - csum 

3469 

3470 return zi 

3471 

3472 

3473def sosfilt_zi(sos): 

3474 """ 

3475 Construct initial conditions for sosfilt for step response steady-state. 

3476 

3477 Compute an initial state `zi` for the `sosfilt` function that corresponds 

3478 to the steady state of the step response. 

3479 

3480 A typical use of this function is to set the initial state so that the 

3481 output of the filter starts at the same value as the first element of 

3482 the signal to be filtered. 

3483 

3484 Parameters 

3485 ---------- 

3486 sos : array_like 

3487 Array of second-order filter coefficients, must have shape 

3488 ``(n_sections, 6)``. See `sosfilt` for the SOS filter format 

3489 specification. 

3490 

3491 Returns 

3492 ------- 

3493 zi : ndarray 

3494 Initial conditions suitable for use with ``sosfilt``, shape 

3495 ``(n_sections, 2)``. 

3496 

3497 See Also 

3498 -------- 

3499 sosfilt, zpk2sos 

3500 

3501 Notes 

3502 ----- 

3503 .. versionadded:: 0.16.0 

3504 

3505 Examples 

3506 -------- 

3507 Filter a rectangular pulse that begins at time 0, with and without 

3508 the use of the `zi` argument of `scipy.signal.sosfilt`. 

3509 

3510 >>> from scipy import signal 

3511 >>> import matplotlib.pyplot as plt 

3512 

3513 >>> sos = signal.butter(9, 0.125, output='sos') 

3514 >>> zi = signal.sosfilt_zi(sos) 

3515 >>> x = (np.arange(250) < 100).astype(int) 

3516 >>> f1 = signal.sosfilt(sos, x) 

3517 >>> f2, zo = signal.sosfilt(sos, x, zi=zi) 

3518 

3519 >>> plt.plot(x, 'k--', label='x') 

3520 >>> plt.plot(f1, 'b', alpha=0.5, linewidth=2, label='filtered') 

3521 >>> plt.plot(f2, 'g', alpha=0.25, linewidth=4, label='filtered with zi') 

3522 >>> plt.legend(loc='best') 

3523 >>> plt.show() 

3524 

3525 """ 

3526 sos = np.asarray(sos) 

3527 if sos.ndim != 2 or sos.shape[1] != 6: 

3528 raise ValueError('sos must be shape (n_sections, 6)') 

3529 

3530 n_sections = sos.shape[0] 

3531 zi = np.empty((n_sections, 2)) 

3532 scale = 1.0 

3533 for section in range(n_sections): 

3534 b = sos[section, :3] 

3535 a = sos[section, 3:] 

3536 zi[section] = scale * lfilter_zi(b, a) 

3537 # If H(z) = B(z)/A(z) is this section's transfer function, then 

3538 # b.sum()/a.sum() is H(1), the gain at omega=0. That's the steady 

3539 # state value of this section's step response. 

3540 scale *= b.sum() / a.sum() 

3541 

3542 return zi 

3543 

3544 

3545def _filtfilt_gust(b, a, x, axis=-1, irlen=None): 

3546 """Forward-backward IIR filter that uses Gustafsson's method. 

3547 

3548 Apply the IIR filter defined by `(b,a)` to `x` twice, first forward 

3549 then backward, using Gustafsson's initial conditions [1]_. 

3550 

3551 Let ``y_fb`` be the result of filtering first forward and then backward, 

3552 and let ``y_bf`` be the result of filtering first backward then forward. 

3553 Gustafsson's method is to compute initial conditions for the forward 

3554 pass and the backward pass such that ``y_fb == y_bf``. 

3555 

3556 Parameters 

3557 ---------- 

3558 b : scalar or 1-D ndarray 

3559 Numerator coefficients of the filter. 

3560 a : scalar or 1-D ndarray 

3561 Denominator coefficients of the filter. 

3562 x : ndarray 

3563 Data to be filtered. 

3564 axis : int, optional 

3565 Axis of `x` to be filtered. Default is -1. 

3566 irlen : int or None, optional 

3567 The length of the nonnegligible part of the impulse response. 

3568 If `irlen` is None, or if the length of the signal is less than 

3569 ``2 * irlen``, then no part of the impulse response is ignored. 

3570 

3571 Returns 

3572 ------- 

3573 y : ndarray 

3574 The filtered data. 

3575 x0 : ndarray 

3576 Initial condition for the forward filter. 

3577 x1 : ndarray 

3578 Initial condition for the backward filter. 

3579 

3580 Notes 

3581 ----- 

3582 Typically the return values `x0` and `x1` are not needed by the 

3583 caller. The intended use of these return values is in unit tests. 

3584 

3585 References 

3586 ---------- 

3587 .. [1] F. Gustaffson. Determining the initial states in forward-backward 

3588 filtering. Transactions on Signal Processing, 46(4):988-992, 1996. 

3589 

3590 """ 

3591 # In the comments, "Gustafsson's paper" and [1] refer to the 

3592 # paper referenced in the docstring. 

3593 

3594 b = np.atleast_1d(b) 

3595 a = np.atleast_1d(a) 

3596 

3597 order = max(len(b), len(a)) - 1 

3598 if order == 0: 

3599 # The filter is just scalar multiplication, with no state. 

3600 scale = (b[0] / a[0])**2 

3601 y = scale * x 

3602 return y, np.array([]), np.array([]) 

3603 

3604 if axis != -1 or axis != x.ndim - 1: 

3605 # Move the axis containing the data to the end. 

3606 x = np.swapaxes(x, axis, x.ndim - 1) 

3607 

3608 # n is the number of samples in the data to be filtered. 

3609 n = x.shape[-1] 

3610 

3611 if irlen is None or n <= 2*irlen: 

3612 m = n 

3613 else: 

3614 m = irlen 

3615 

3616 # Create Obs, the observability matrix (called O in the paper). 

3617 # This matrix can be interpreted as the operator that propagates 

3618 # an arbitrary initial state to the output, assuming the input is 

3619 # zero. 

3620 # In Gustafsson's paper, the forward and backward filters are not 

3621 # necessarily the same, so he has both O_f and O_b. We use the same 

3622 # filter in both directions, so we only need O. The same comment 

3623 # applies to S below. 

3624 Obs = np.zeros((m, order)) 

3625 zi = np.zeros(order) 

3626 zi[0] = 1 

3627 Obs[:, 0] = lfilter(b, a, np.zeros(m), zi=zi)[0] 

3628 for k in range(1, order): 

3629 Obs[k:, k] = Obs[:-k, 0] 

3630 

3631 # Obsr is O^R (Gustafsson's notation for row-reversed O) 

3632 Obsr = Obs[::-1] 

3633 

3634 # Create S. S is the matrix that applies the filter to the reversed 

3635 # propagated initial conditions. That is, 

3636 # out = S.dot(zi) 

3637 # is the same as 

3638 # tmp, _ = lfilter(b, a, zeros(), zi=zi) # Propagate ICs. 

3639 # out = lfilter(b, a, tmp[::-1]) # Reverse and filter. 

3640 

3641 # Equations (5) & (6) of [1] 

3642 S = lfilter(b, a, Obs[::-1], axis=0) 

3643 

3644 # Sr is S^R (row-reversed S) 

3645 Sr = S[::-1] 

3646 

3647 # M is [(S^R - O), (O^R - S)] 

3648 if m == n: 

3649 M = np.hstack((Sr - Obs, Obsr - S)) 

3650 else: 

3651 # Matrix described in section IV of [1]. 

3652 M = np.zeros((2*m, 2*order)) 

3653 M[:m, :order] = Sr - Obs 

3654 M[m:, order:] = Obsr - S 

3655 

3656 # Naive forward-backward and backward-forward filters. 

3657 # These have large transients because the filters use zero initial 

3658 # conditions. 

3659 y_f = lfilter(b, a, x) 

3660 y_fb = lfilter(b, a, y_f[..., ::-1])[..., ::-1] 

3661 

3662 y_b = lfilter(b, a, x[..., ::-1])[..., ::-1] 

3663 y_bf = lfilter(b, a, y_b) 

3664 

3665 delta_y_bf_fb = y_bf - y_fb 

3666 if m == n: 

3667 delta = delta_y_bf_fb 

3668 else: 

3669 start_m = delta_y_bf_fb[..., :m] 

3670 end_m = delta_y_bf_fb[..., -m:] 

3671 delta = np.concatenate((start_m, end_m), axis=-1) 

3672 

3673 # ic_opt holds the "optimal" initial conditions. 

3674 # The following code computes the result shown in the formula 

3675 # of the paper between equations (6) and (7). 

3676 if delta.ndim == 1: 

3677 ic_opt = linalg.lstsq(M, delta)[0] 

3678 else: 

3679 # Reshape delta so it can be used as an array of multiple 

3680 # right-hand-sides in linalg.lstsq. 

3681 delta2d = delta.reshape(-1, delta.shape[-1]).T 

3682 ic_opt0 = linalg.lstsq(M, delta2d)[0].T 

3683 ic_opt = ic_opt0.reshape(delta.shape[:-1] + (M.shape[-1],)) 

3684 

3685 # Now compute the filtered signal using equation (7) of [1]. 

3686 # First, form [S^R, O^R] and call it W. 

3687 if m == n: 

3688 W = np.hstack((Sr, Obsr)) 

3689 else: 

3690 W = np.zeros((2*m, 2*order)) 

3691 W[:m, :order] = Sr 

3692 W[m:, order:] = Obsr 

3693 

3694 # Equation (7) of [1] says 

3695 # Y_fb^opt = Y_fb^0 + W * [x_0^opt; x_{N-1}^opt] 

3696 # `wic` is (almost) the product on the right. 

3697 # W has shape (m, 2*order), and ic_opt has shape (..., 2*order), 

3698 # so we can't use W.dot(ic_opt). Instead, we dot ic_opt with W.T, 

3699 # so wic has shape (..., m). 

3700 wic = ic_opt.dot(W.T) 

3701 

3702 # `wic` is "almost" the product of W and the optimal ICs in equation 

3703 # (7)--if we're using a truncated impulse response (m < n), `wic` 

3704 # contains only the adjustments required for the ends of the signal. 

3705 # Here we form y_opt, taking this into account if necessary. 

3706 y_opt = y_fb 

3707 if m == n: 

3708 y_opt += wic 

3709 else: 

3710 y_opt[..., :m] += wic[..., :m] 

3711 y_opt[..., -m:] += wic[..., -m:] 

3712 

3713 x0 = ic_opt[..., :order] 

3714 x1 = ic_opt[..., -order:] 

3715 if axis != -1 or axis != x.ndim - 1: 

3716 # Restore the data axis to its original position. 

3717 x0 = np.swapaxes(x0, axis, x.ndim - 1) 

3718 x1 = np.swapaxes(x1, axis, x.ndim - 1) 

3719 y_opt = np.swapaxes(y_opt, axis, x.ndim - 1) 

3720 

3721 return y_opt, x0, x1 

3722 

3723 

3724def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad', 

3725 irlen=None): 

3726 """ 

3727 Apply a digital filter forward and backward to a signal. 

3728 

3729 This function applies a linear digital filter twice, once forward and 

3730 once backwards. The combined filter has zero phase and a filter order 

3731 twice that of the original. 

3732 

3733 The function provides options for handling the edges of the signal. 

3734 

3735 The function `sosfiltfilt` (and filter design using ``output='sos'``) 

3736 should be preferred over `filtfilt` for most filtering tasks, as 

3737 second-order sections have fewer numerical problems. 

3738 

3739 Parameters 

3740 ---------- 

3741 b : (N,) array_like 

3742 The numerator coefficient vector of the filter. 

3743 a : (N,) array_like 

3744 The denominator coefficient vector of the filter. If ``a[0]`` 

3745 is not 1, then both `a` and `b` are normalized by ``a[0]``. 

3746 x : array_like 

3747 The array of data to be filtered. 

3748 axis : int, optional 

3749 The axis of `x` to which the filter is applied. 

3750 Default is -1. 

3751 padtype : str or None, optional 

3752 Must be 'odd', 'even', 'constant', or None. This determines the 

3753 type of extension to use for the padded signal to which the filter 

3754 is applied. If `padtype` is None, no padding is used. The default 

3755 is 'odd'. 

3756 padlen : int or None, optional 

3757 The number of elements by which to extend `x` at both ends of 

3758 `axis` before applying the filter. This value must be less than 

3759 ``x.shape[axis] - 1``. ``padlen=0`` implies no padding. 

3760 The default value is ``3 * max(len(a), len(b))``. 

3761 method : str, optional 

3762 Determines the method for handling the edges of the signal, either 

3763 "pad" or "gust". When `method` is "pad", the signal is padded; the 

3764 type of padding is determined by `padtype` and `padlen`, and `irlen` 

3765 is ignored. When `method` is "gust", Gustafsson's method is used, 

3766 and `padtype` and `padlen` are ignored. 

3767 irlen : int or None, optional 

3768 When `method` is "gust", `irlen` specifies the length of the 

3769 impulse response of the filter. If `irlen` is None, no part 

3770 of the impulse response is ignored. For a long signal, specifying 

3771 `irlen` can significantly improve the performance of the filter. 

3772 

3773 Returns 

3774 ------- 

3775 y : ndarray 

3776 The filtered output with the same shape as `x`. 

3777 

3778 See Also 

3779 -------- 

3780 sosfiltfilt, lfilter_zi, lfilter, lfiltic, savgol_filter, sosfilt 

3781 

3782 Notes 

3783 ----- 

3784 When `method` is "pad", the function pads the data along the given axis 

3785 in one of three ways: odd, even or constant. The odd and even extensions 

3786 have the corresponding symmetry about the end point of the data. The 

3787 constant extension extends the data with the values at the end points. On 

3788 both the forward and backward passes, the initial condition of the 

3789 filter is found by using `lfilter_zi` and scaling it by the end point of 

3790 the extended data. 

3791 

3792 When `method` is "gust", Gustafsson's method [1]_ is used. Initial 

3793 conditions are chosen for the forward and backward passes so that the 

3794 forward-backward filter gives the same result as the backward-forward 

3795 filter. 

3796 

3797 The option to use Gustaffson's method was added in scipy version 0.16.0. 

3798 

3799 References 

3800 ---------- 

3801 .. [1] F. Gustaffson, "Determining the initial states in forward-backward 

3802 filtering", Transactions on Signal Processing, Vol. 46, pp. 988-992, 

3803 1996. 

3804 

3805 Examples 

3806 -------- 

3807 The examples will use several functions from `scipy.signal`. 

3808 

3809 >>> from scipy import signal 

3810 >>> import matplotlib.pyplot as plt 

3811 

3812 First we create a one second signal that is the sum of two pure sine 

3813 waves, with frequencies 5 Hz and 250 Hz, sampled at 2000 Hz. 

3814 

3815 >>> t = np.linspace(0, 1.0, 2001) 

3816 >>> xlow = np.sin(2 * np.pi * 5 * t) 

3817 >>> xhigh = np.sin(2 * np.pi * 250 * t) 

3818 >>> x = xlow + xhigh 

3819 

3820 Now create a lowpass Butterworth filter with a cutoff of 0.125 times 

3821 the Nyquist frequency, or 125 Hz, and apply it to ``x`` with `filtfilt`. 

3822 The result should be approximately ``xlow``, with no phase shift. 

3823 

3824 >>> b, a = signal.butter(8, 0.125) 

3825 >>> y = signal.filtfilt(b, a, x, padlen=150) 

3826 >>> np.abs(y - xlow).max() 

3827 9.1086182074789912e-06 

3828 

3829 We get a fairly clean result for this artificial example because 

3830 the odd extension is exact, and with the moderately long padding, 

3831 the filter's transients have dissipated by the time the actual data 

3832 is reached. In general, transient effects at the edges are 

3833 unavoidable. 

3834 

3835 The following example demonstrates the option ``method="gust"``. 

3836 

3837 First, create a filter. 

3838 

3839 >>> b, a = signal.ellip(4, 0.01, 120, 0.125) # Filter to be applied. 

3840 >>> np.random.seed(123456) 

3841 

3842 `sig` is a random input signal to be filtered. 

3843 

3844 >>> n = 60 

3845 >>> sig = np.random.randn(n)**3 + 3*np.random.randn(n).cumsum() 

3846 

3847 Apply `filtfilt` to `sig`, once using the Gustafsson method, and 

3848 once using padding, and plot the results for comparison. 

3849 

3850 >>> fgust = signal.filtfilt(b, a, sig, method="gust") 

3851 >>> fpad = signal.filtfilt(b, a, sig, padlen=50) 

3852 >>> plt.plot(sig, 'k-', label='input') 

3853 >>> plt.plot(fgust, 'b-', linewidth=4, label='gust') 

3854 >>> plt.plot(fpad, 'c-', linewidth=1.5, label='pad') 

3855 >>> plt.legend(loc='best') 

3856 >>> plt.show() 

3857 

3858 The `irlen` argument can be used to improve the performance 

3859 of Gustafsson's method. 

3860 

3861 Estimate the impulse response length of the filter. 

3862 

3863 >>> z, p, k = signal.tf2zpk(b, a) 

3864 >>> eps = 1e-9 

3865 >>> r = np.max(np.abs(p)) 

3866 >>> approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r))) 

3867 >>> approx_impulse_len 

3868 137 

3869 

3870 Apply the filter to a longer signal, with and without the `irlen` 

3871 argument. The difference between `y1` and `y2` is small. For long 

3872 signals, using `irlen` gives a significant performance improvement. 

3873 

3874 >>> x = np.random.randn(5000) 

3875 >>> y1 = signal.filtfilt(b, a, x, method='gust') 

3876 >>> y2 = signal.filtfilt(b, a, x, method='gust', irlen=approx_impulse_len) 

3877 >>> print(np.max(np.abs(y1 - y2))) 

3878 1.80056858312e-10 

3879 

3880 """ 

3881 b = np.atleast_1d(b) 

3882 a = np.atleast_1d(a) 

3883 x = np.asarray(x) 

3884 

3885 if method not in ["pad", "gust"]: 

3886 raise ValueError("method must be 'pad' or 'gust'.") 

3887 

3888 if method == "gust": 

3889 y, z1, z2 = _filtfilt_gust(b, a, x, axis=axis, irlen=irlen) 

3890 return y 

3891 

3892 # method == "pad" 

3893 edge, ext = _validate_pad(padtype, padlen, x, axis, 

3894 ntaps=max(len(a), len(b))) 

3895 

3896 # Get the steady state of the filter's step response. 

3897 zi = lfilter_zi(b, a) 

3898 

3899 # Reshape zi and create x0 so that zi*x0 broadcasts 

3900 # to the correct value for the 'zi' keyword argument 

3901 # to lfilter. 

3902 zi_shape = [1] * x.ndim 

3903 zi_shape[axis] = zi.size 

3904 zi = np.reshape(zi, zi_shape) 

3905 x0 = axis_slice(ext, stop=1, axis=axis) 

3906 

3907 # Forward filter. 

3908 (y, zf) = lfilter(b, a, ext, axis=axis, zi=zi * x0) 

3909 

3910 # Backward filter. 

3911 # Create y0 so zi*y0 broadcasts appropriately. 

3912 y0 = axis_slice(y, start=-1, axis=axis) 

3913 (y, zf) = lfilter(b, a, axis_reverse(y, axis=axis), axis=axis, zi=zi * y0) 

3914 

3915 # Reverse y. 

3916 y = axis_reverse(y, axis=axis) 

3917 

3918 if edge > 0: 

3919 # Slice the actual signal from the extended signal. 

3920 y = axis_slice(y, start=edge, stop=-edge, axis=axis) 

3921 

3922 return y 

3923 

3924 

3925def _validate_pad(padtype, padlen, x, axis, ntaps): 

3926 """Helper to validate padding for filtfilt""" 

3927 if padtype not in ['even', 'odd', 'constant', None]: 

3928 raise ValueError(("Unknown value '%s' given to padtype. padtype " 

3929 "must be 'even', 'odd', 'constant', or None.") % 

3930 padtype) 

3931 

3932 if padtype is None: 

3933 padlen = 0 

3934 

3935 if padlen is None: 

3936 # Original padding; preserved for backwards compatibility. 

3937 edge = ntaps * 3 

3938 else: 

3939 edge = padlen 

3940 

3941 # x's 'axis' dimension must be bigger than edge. 

3942 if x.shape[axis] <= edge: 

3943 raise ValueError("The length of the input vector x must be greater " 

3944 "than padlen, which is %d." % edge) 

3945 

3946 if padtype is not None and edge > 0: 

3947 # Make an extension of length `edge` at each 

3948 # end of the input array. 

3949 if padtype == 'even': 

3950 ext = even_ext(x, edge, axis=axis) 

3951 elif padtype == 'odd': 

3952 ext = odd_ext(x, edge, axis=axis) 

3953 else: 

3954 ext = const_ext(x, edge, axis=axis) 

3955 else: 

3956 ext = x 

3957 return edge, ext 

3958 

3959 

3960def _validate_x(x): 

3961 x = np.asarray(x) 

3962 if x.ndim == 0: 

3963 raise ValueError('x must be at least 1-D') 

3964 return x 

3965 

3966 

3967def sosfilt(sos, x, axis=-1, zi=None): 

3968 """ 

3969 Filter data along one dimension using cascaded second-order sections. 

3970 

3971 Filter a data sequence, `x`, using a digital IIR filter defined by 

3972 `sos`. 

3973 

3974 Parameters 

3975 ---------- 

3976 sos : array_like 

3977 Array of second-order filter coefficients, must have shape 

3978 ``(n_sections, 6)``. Each row corresponds to a second-order 

3979 section, with the first three columns providing the numerator 

3980 coefficients and the last three providing the denominator 

3981 coefficients. 

3982 x : array_like 

3983 An N-dimensional input array. 

3984 axis : int, optional 

3985 The axis of the input data array along which to apply the 

3986 linear filter. The filter is applied to each subarray along 

3987 this axis. Default is -1. 

3988 zi : array_like, optional 

3989 Initial conditions for the cascaded filter delays. It is a (at 

3990 least 2D) vector of shape ``(n_sections, ..., 2, ...)``, where 

3991 ``..., 2, ...`` denotes the shape of `x`, but with ``x.shape[axis]`` 

3992 replaced by 2. If `zi` is None or is not given then initial rest 

3993 (i.e. all zeros) is assumed. 

3994 Note that these initial conditions are *not* the same as the initial 

3995 conditions given by `lfiltic` or `lfilter_zi`. 

3996 

3997 Returns 

3998 ------- 

3999 y : ndarray 

4000 The output of the digital filter. 

4001 zf : ndarray, optional 

4002 If `zi` is None, this is not returned, otherwise, `zf` holds the 

4003 final filter delay values. 

4004 

4005 See Also 

4006 -------- 

4007 zpk2sos, sos2zpk, sosfilt_zi, sosfiltfilt, sosfreqz 

4008 

4009 Notes 

4010 ----- 

4011 The filter function is implemented as a series of second-order filters 

4012 with direct-form II transposed structure. It is designed to minimize 

4013 numerical precision errors for high-order filters. 

4014 

4015 .. versionadded:: 0.16.0 

4016 

4017 Examples 

4018 -------- 

4019 Plot a 13th-order filter's impulse response using both `lfilter` and 

4020 `sosfilt`, showing the instability that results from trying to do a 

4021 13th-order filter in a single stage (the numerical error pushes some poles 

4022 outside of the unit circle): 

4023 

4024 >>> import matplotlib.pyplot as plt 

4025 >>> from scipy import signal 

4026 >>> b, a = signal.ellip(13, 0.009, 80, 0.05, output='ba') 

4027 >>> sos = signal.ellip(13, 0.009, 80, 0.05, output='sos') 

4028 >>> x = signal.unit_impulse(700) 

4029 >>> y_tf = signal.lfilter(b, a, x) 

4030 >>> y_sos = signal.sosfilt(sos, x) 

4031 >>> plt.plot(y_tf, 'r', label='TF') 

4032 >>> plt.plot(y_sos, 'k', label='SOS') 

4033 >>> plt.legend(loc='best') 

4034 >>> plt.show() 

4035 

4036 """ 

4037 x = _validate_x(x) 

4038 sos, n_sections = _validate_sos(sos) 

4039 x_zi_shape = list(x.shape) 

4040 x_zi_shape[axis] = 2 

4041 x_zi_shape = tuple([n_sections] + x_zi_shape) 

4042 inputs = [sos, x] 

4043 if zi is not None: 

4044 inputs.append(np.asarray(zi)) 

4045 dtype = np.result_type(*inputs) 

4046 if dtype.char not in 'fdgFDGO': 

4047 raise NotImplementedError("input type '%s' not supported" % dtype) 

4048 if zi is not None: 

4049 zi = np.array(zi, dtype) # make a copy so that we can operate in place 

4050 if zi.shape != x_zi_shape: 

4051 raise ValueError('Invalid zi shape. With axis=%r, an input with ' 

4052 'shape %r, and an sos array with %d sections, zi ' 

4053 'must have shape %r, got %r.' % 

4054 (axis, x.shape, n_sections, x_zi_shape, zi.shape)) 

4055 return_zi = True 

4056 else: 

4057 zi = np.zeros(x_zi_shape, dtype=dtype) 

4058 return_zi = False 

4059 axis = axis % x.ndim # make positive 

4060 x = np.moveaxis(x, axis, -1) 

4061 zi = np.moveaxis(zi, [0, axis + 1], [-2, -1]) 

4062 x_shape, zi_shape = x.shape, zi.shape 

4063 x = np.reshape(x, (-1, x.shape[-1])) 

4064 x = np.array(x, dtype, order='C') # make a copy, can modify in place 

4065 zi = np.ascontiguousarray(np.reshape(zi, (-1, n_sections, 2))) 

4066 sos = sos.astype(dtype, copy=False) 

4067 _sosfilt(sos, x, zi) 

4068 x.shape = x_shape 

4069 x = np.moveaxis(x, -1, axis) 

4070 if return_zi: 

4071 zi.shape = zi_shape 

4072 zi = np.moveaxis(zi, [-2, -1], [0, axis + 1]) 

4073 out = (x, zi) 

4074 else: 

4075 out = x 

4076 return out 

4077 

4078 

4079def sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None): 

4080 """ 

4081 A forward-backward digital filter using cascaded second-order sections. 

4082 

4083 See `filtfilt` for more complete information about this method. 

4084 

4085 Parameters 

4086 ---------- 

4087 sos : array_like 

4088 Array of second-order filter coefficients, must have shape 

4089 ``(n_sections, 6)``. Each row corresponds to a second-order 

4090 section, with the first three columns providing the numerator 

4091 coefficients and the last three providing the denominator 

4092 coefficients. 

4093 x : array_like 

4094 The array of data to be filtered. 

4095 axis : int, optional 

4096 The axis of `x` to which the filter is applied. 

4097 Default is -1. 

4098 padtype : str or None, optional 

4099 Must be 'odd', 'even', 'constant', or None. This determines the 

4100 type of extension to use for the padded signal to which the filter 

4101 is applied. If `padtype` is None, no padding is used. The default 

4102 is 'odd'. 

4103 padlen : int or None, optional 

4104 The number of elements by which to extend `x` at both ends of 

4105 `axis` before applying the filter. This value must be less than 

4106 ``x.shape[axis] - 1``. ``padlen=0`` implies no padding. 

4107 The default value is:: 

4108 

4109 3 * (2 * len(sos) + 1 - min((sos[:, 2] == 0).sum(), 

4110 (sos[:, 5] == 0).sum())) 

4111 

4112 The extra subtraction at the end attempts to compensate for poles 

4113 and zeros at the origin (e.g. for odd-order filters) to yield 

4114 equivalent estimates of `padlen` to those of `filtfilt` for 

4115 second-order section filters built with `scipy.signal` functions. 

4116 

4117 Returns 

4118 ------- 

4119 y : ndarray 

4120 The filtered output with the same shape as `x`. 

4121 

4122 See Also 

4123 -------- 

4124 filtfilt, sosfilt, sosfilt_zi, sosfreqz 

4125 

4126 Notes 

4127 ----- 

4128 .. versionadded:: 0.18.0 

4129 

4130 Examples 

4131 -------- 

4132 >>> from scipy.signal import sosfiltfilt, butter 

4133 >>> import matplotlib.pyplot as plt 

4134 

4135 Create an interesting signal to filter. 

4136 

4137 >>> n = 201 

4138 >>> t = np.linspace(0, 1, n) 

4139 >>> np.random.seed(123) 

4140 >>> x = 1 + (t < 0.5) - 0.25*t**2 + 0.05*np.random.randn(n) 

4141 

4142 Create a lowpass Butterworth filter, and use it to filter `x`. 

4143 

4144 >>> sos = butter(4, 0.125, output='sos') 

4145 >>> y = sosfiltfilt(sos, x) 

4146 

4147 For comparison, apply an 8th order filter using `sosfilt`. The filter 

4148 is initialized using the mean of the first four values of `x`. 

4149 

4150 >>> from scipy.signal import sosfilt, sosfilt_zi 

4151 >>> sos8 = butter(8, 0.125, output='sos') 

4152 >>> zi = x[:4].mean() * sosfilt_zi(sos8) 

4153 >>> y2, zo = sosfilt(sos8, x, zi=zi) 

4154 

4155 Plot the results. Note that the phase of `y` matches the input, while 

4156 `y2` has a significant phase delay. 

4157 

4158 >>> plt.plot(t, x, alpha=0.5, label='x(t)') 

4159 >>> plt.plot(t, y, label='y(t)') 

4160 >>> plt.plot(t, y2, label='y2(t)') 

4161 >>> plt.legend(framealpha=1, shadow=True) 

4162 >>> plt.grid(alpha=0.25) 

4163 >>> plt.xlabel('t') 

4164 >>> plt.show() 

4165 

4166 """ 

4167 sos, n_sections = _validate_sos(sos) 

4168 x = _validate_x(x) 

4169 

4170 # `method` is "pad"... 

4171 ntaps = 2 * n_sections + 1 

4172 ntaps -= min((sos[:, 2] == 0).sum(), (sos[:, 5] == 0).sum()) 

4173 edge, ext = _validate_pad(padtype, padlen, x, axis, 

4174 ntaps=ntaps) 

4175 

4176 # These steps follow the same form as filtfilt with modifications 

4177 zi = sosfilt_zi(sos) # shape (n_sections, 2) --> (n_sections, ..., 2, ...) 

4178 zi_shape = [1] * x.ndim 

4179 zi_shape[axis] = 2 

4180 zi.shape = [n_sections] + zi_shape 

4181 x_0 = axis_slice(ext, stop=1, axis=axis) 

4182 (y, zf) = sosfilt(sos, ext, axis=axis, zi=zi * x_0) 

4183 y_0 = axis_slice(y, start=-1, axis=axis) 

4184 (y, zf) = sosfilt(sos, axis_reverse(y, axis=axis), axis=axis, zi=zi * y_0) 

4185 y = axis_reverse(y, axis=axis) 

4186 if edge > 0: 

4187 y = axis_slice(y, start=edge, stop=-edge, axis=axis) 

4188 return y 

4189 

4190 

4191def decimate(x, q, n=None, ftype='iir', axis=-1, zero_phase=True): 

4192 """ 

4193 Downsample the signal after applying an anti-aliasing filter. 

4194 

4195 By default, an order 8 Chebyshev type I filter is used. A 30 point FIR 

4196 filter with Hamming window is used if `ftype` is 'fir'. 

4197 

4198 Parameters 

4199 ---------- 

4200 x : array_like 

4201 The signal to be downsampled, as an N-dimensional array. 

4202 q : int 

4203 The downsampling factor. When using IIR downsampling, it is recommended 

4204 to call `decimate` multiple times for downsampling factors higher than 

4205 13. 

4206 n : int, optional 

4207 The order of the filter (1 less than the length for 'fir'). Defaults to 

4208 8 for 'iir' and 20 times the downsampling factor for 'fir'. 

4209 ftype : str {'iir', 'fir'} or ``dlti`` instance, optional 

4210 If 'iir' or 'fir', specifies the type of lowpass filter. If an instance 

4211 of an `dlti` object, uses that object to filter before downsampling. 

4212 axis : int, optional 

4213 The axis along which to decimate. 

4214 zero_phase : bool, optional 

4215 Prevent phase shift by filtering with `filtfilt` instead of `lfilter` 

4216 when using an IIR filter, and shifting the outputs back by the filter's 

4217 group delay when using an FIR filter. The default value of ``True`` is 

4218 recommended, since a phase shift is generally not desired. 

4219 

4220 .. versionadded:: 0.18.0 

4221 

4222 Returns 

4223 ------- 

4224 y : ndarray 

4225 The down-sampled signal. 

4226 

4227 See Also 

4228 -------- 

4229 resample : Resample up or down using the FFT method. 

4230 resample_poly : Resample using polyphase filtering and an FIR filter. 

4231 

4232 Notes 

4233 ----- 

4234 The ``zero_phase`` keyword was added in 0.18.0. 

4235 The possibility to use instances of ``dlti`` as ``ftype`` was added in 

4236 0.18.0. 

4237 """ 

4238 

4239 x = np.asarray(x) 

4240 q = operator.index(q) 

4241 

4242 if n is not None: 

4243 n = operator.index(n) 

4244 

4245 if ftype == 'fir': 

4246 if n is None: 

4247 half_len = 10 * q # reasonable cutoff for our sinc-like function 

4248 n = 2 * half_len 

4249 b, a = firwin(n+1, 1. / q, window='hamming'), 1. 

4250 elif ftype == 'iir': 

4251 if n is None: 

4252 n = 8 

4253 system = dlti(*cheby1(n, 0.05, 0.8 / q)) 

4254 b, a = system.num, system.den 

4255 elif isinstance(ftype, dlti): 

4256 system = ftype._as_tf() # Avoids copying if already in TF form 

4257 b, a = system.num, system.den 

4258 else: 

4259 raise ValueError('invalid ftype') 

4260 

4261 sl = [slice(None)] * x.ndim 

4262 a = np.asarray(a) 

4263 

4264 if a.size == 1: # FIR case 

4265 b = b / a 

4266 if zero_phase: 

4267 y = resample_poly(x, 1, q, axis=axis, window=b) 

4268 else: 

4269 # upfirdn is generally faster than lfilter by a factor equal to the 

4270 # downsampling factor, since it only calculates the needed outputs 

4271 n_out = x.shape[axis] // q + bool(x.shape[axis] % q) 

4272 y = upfirdn(b, x, up=1, down=q, axis=axis) 

4273 sl[axis] = slice(None, n_out, None) 

4274 

4275 else: # IIR case 

4276 if zero_phase: 

4277 y = filtfilt(b, a, x, axis=axis) 

4278 else: 

4279 y = lfilter(b, a, x, axis=axis) 

4280 sl[axis] = slice(None, None, q) 

4281 

4282 return y[tuple(sl)]