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# -*- coding: utf-8 -*- 

2"""Functions for FIR filter design.""" 

3from math import ceil, log 

4import operator 

5import warnings 

6 

7import numpy as np 

8from numpy.fft import irfft, fft, ifft 

9from scipy.special import sinc 

10from scipy.linalg import (toeplitz, hankel, solve, LinAlgError, LinAlgWarning, 

11 lstsq) 

12 

13from . import sigtools 

14 

15__all__ = ['kaiser_beta', 'kaiser_atten', 'kaiserord', 

16 'firwin', 'firwin2', 'remez', 'firls', 'minimum_phase'] 

17 

18 

19def _get_fs(fs, nyq): 

20 """ 

21 Utility for replacing the argument 'nyq' (with default 1) with 'fs'. 

22 """ 

23 if nyq is None and fs is None: 

24 fs = 2 

25 elif nyq is not None: 

26 if fs is not None: 

27 raise ValueError("Values cannot be given for both 'nyq' and 'fs'.") 

28 fs = 2*nyq 

29 return fs 

30 

31 

32# Some notes on function parameters: 

33# 

34# `cutoff` and `width` are given as numbers between 0 and 1. These are 

35# relative frequencies, expressed as a fraction of the Nyquist frequency. 

36# For example, if the Nyquist frequency is 2 KHz, then width=0.15 is a width 

37# of 300 Hz. 

38# 

39# The `order` of a FIR filter is one less than the number of taps. 

40# This is a potential source of confusion, so in the following code, 

41# we will always use the number of taps as the parameterization of 

42# the 'size' of the filter. The "number of taps" means the number 

43# of coefficients, which is the same as the length of the impulse 

44# response of the filter. 

45 

46 

47def kaiser_beta(a): 

48 """Compute the Kaiser parameter `beta`, given the attenuation `a`. 

49 

50 Parameters 

51 ---------- 

52 a : float 

53 The desired attenuation in the stopband and maximum ripple in 

54 the passband, in dB. This should be a *positive* number. 

55 

56 Returns 

57 ------- 

58 beta : float 

59 The `beta` parameter to be used in the formula for a Kaiser window. 

60 

61 References 

62 ---------- 

63 Oppenheim, Schafer, "Discrete-Time Signal Processing", p.475-476. 

64 

65 Examples 

66 -------- 

67 Suppose we want to design a lowpass filter, with 65 dB attenuation 

68 in the stop band. The Kaiser window parameter to be used in the 

69 window method is computed by `kaiser_beta(65)`: 

70 

71 >>> from scipy.signal import kaiser_beta 

72 >>> kaiser_beta(65) 

73 6.20426 

74 

75 """ 

76 if a > 50: 

77 beta = 0.1102 * (a - 8.7) 

78 elif a > 21: 

79 beta = 0.5842 * (a - 21) ** 0.4 + 0.07886 * (a - 21) 

80 else: 

81 beta = 0.0 

82 return beta 

83 

84 

85def kaiser_atten(numtaps, width): 

86 """Compute the attenuation of a Kaiser FIR filter. 

87 

88 Given the number of taps `N` and the transition width `width`, compute the 

89 attenuation `a` in dB, given by Kaiser's formula: 

90 

91 a = 2.285 * (N - 1) * pi * width + 7.95 

92 

93 Parameters 

94 ---------- 

95 numtaps : int 

96 The number of taps in the FIR filter. 

97 width : float 

98 The desired width of the transition region between passband and 

99 stopband (or, in general, at any discontinuity) for the filter, 

100 expressed as a fraction of the Nyquist frequency. 

101 

102 Returns 

103 ------- 

104 a : float 

105 The attenuation of the ripple, in dB. 

106 

107 See Also 

108 -------- 

109 kaiserord, kaiser_beta 

110 

111 Examples 

112 -------- 

113 Suppose we want to design a FIR filter using the Kaiser window method 

114 that will have 211 taps and a transition width of 9 Hz for a signal that 

115 is sampled at 480 Hz. Expressed as a fraction of the Nyquist frequency, 

116 the width is 9/(0.5*480) = 0.0375. The approximate attenuation (in dB) 

117 is computed as follows: 

118 

119 >>> from scipy.signal import kaiser_atten 

120 >>> kaiser_atten(211, 0.0375) 

121 64.48099630593983 

122 

123 """ 

124 a = 2.285 * (numtaps - 1) * np.pi * width + 7.95 

125 return a 

126 

127 

128def kaiserord(ripple, width): 

129 """ 

130 Determine the filter window parameters for the Kaiser window method. 

131 

132 The parameters returned by this function are generally used to create 

133 a finite impulse response filter using the window method, with either 

134 `firwin` or `firwin2`. 

135 

136 Parameters 

137 ---------- 

138 ripple : float 

139 Upper bound for the deviation (in dB) of the magnitude of the 

140 filter's frequency response from that of the desired filter (not 

141 including frequencies in any transition intervals). That is, if w 

142 is the frequency expressed as a fraction of the Nyquist frequency, 

143 A(w) is the actual frequency response of the filter and D(w) is the 

144 desired frequency response, the design requirement is that:: 

145 

146 abs(A(w) - D(w))) < 10**(-ripple/20) 

147 

148 for 0 <= w <= 1 and w not in a transition interval. 

149 width : float 

150 Width of transition region, normalized so that 1 corresponds to pi 

151 radians / sample. That is, the frequency is expressed as a fraction 

152 of the Nyquist frequency. 

153 

154 Returns 

155 ------- 

156 numtaps : int 

157 The length of the Kaiser window. 

158 beta : float 

159 The beta parameter for the Kaiser window. 

160 

161 See Also 

162 -------- 

163 kaiser_beta, kaiser_atten 

164 

165 Notes 

166 ----- 

167 There are several ways to obtain the Kaiser window: 

168 

169 - ``signal.kaiser(numtaps, beta, sym=True)`` 

170 - ``signal.get_window(beta, numtaps)`` 

171 - ``signal.get_window(('kaiser', beta), numtaps)`` 

172 

173 The empirical equations discovered by Kaiser are used. 

174 

175 References 

176 ---------- 

177 Oppenheim, Schafer, "Discrete-Time Signal Processing", pp.475-476. 

178 

179 Examples 

180 -------- 

181 We will use the Kaiser window method to design a lowpass FIR filter 

182 for a signal that is sampled at 1000 Hz. 

183 

184 We want at least 65 dB rejection in the stop band, and in the pass 

185 band the gain should vary no more than 0.5%. 

186 

187 We want a cutoff frequency of 175 Hz, with a transition between the 

188 pass band and the stop band of 24 Hz. That is, in the band [0, 163], 

189 the gain varies no more than 0.5%, and in the band [187, 500], the 

190 signal is attenuated by at least 65 dB. 

191 

192 >>> from scipy.signal import kaiserord, firwin, freqz 

193 >>> import matplotlib.pyplot as plt 

194 >>> fs = 1000.0 

195 >>> cutoff = 175 

196 >>> width = 24 

197 

198 The Kaiser method accepts just a single parameter to control the pass 

199 band ripple and the stop band rejection, so we use the more restrictive 

200 of the two. In this case, the pass band ripple is 0.005, or 46.02 dB, 

201 so we will use 65 dB as the design parameter. 

202 

203 Use `kaiserord` to determine the length of the filter and the 

204 parameter for the Kaiser window. 

205 

206 >>> numtaps, beta = kaiserord(65, width/(0.5*fs)) 

207 >>> numtaps 

208 167 

209 >>> beta 

210 6.20426 

211 

212 Use `firwin` to create the FIR filter. 

213 

214 >>> taps = firwin(numtaps, cutoff, window=('kaiser', beta), 

215 ... scale=False, nyq=0.5*fs) 

216 

217 Compute the frequency response of the filter. ``w`` is the array of 

218 frequencies, and ``h`` is the corresponding complex array of frequency 

219 responses. 

220 

221 >>> w, h = freqz(taps, worN=8000) 

222 >>> w *= 0.5*fs/np.pi # Convert w to Hz. 

223 

224 Compute the deviation of the magnitude of the filter's response from 

225 that of the ideal lowpass filter. Values in the transition region are 

226 set to ``nan``, so they won't appear in the plot. 

227 

228 >>> ideal = w < cutoff # The "ideal" frequency response. 

229 >>> deviation = np.abs(np.abs(h) - ideal) 

230 >>> deviation[(w > cutoff - 0.5*width) & (w < cutoff + 0.5*width)] = np.nan 

231 

232 Plot the deviation. A close look at the left end of the stop band shows 

233 that the requirement for 65 dB attenuation is violated in the first lobe 

234 by about 0.125 dB. This is not unusual for the Kaiser window method. 

235 

236 >>> plt.plot(w, 20*np.log10(np.abs(deviation))) 

237 >>> plt.xlim(0, 0.5*fs) 

238 >>> plt.ylim(-90, -60) 

239 >>> plt.grid(alpha=0.25) 

240 >>> plt.axhline(-65, color='r', ls='--', alpha=0.3) 

241 >>> plt.xlabel('Frequency (Hz)') 

242 >>> plt.ylabel('Deviation from ideal (dB)') 

243 >>> plt.title('Lowpass Filter Frequency Response') 

244 >>> plt.show() 

245 

246 """ 

247 A = abs(ripple) # in case somebody is confused as to what's meant 

248 if A < 8: 

249 # Formula for N is not valid in this range. 

250 raise ValueError("Requested maximum ripple attentuation %f is too " 

251 "small for the Kaiser formula." % A) 

252 beta = kaiser_beta(A) 

253 

254 # Kaiser's formula (as given in Oppenheim and Schafer) is for the filter 

255 # order, so we have to add 1 to get the number of taps. 

256 numtaps = (A - 7.95) / 2.285 / (np.pi * width) + 1 

257 

258 return int(ceil(numtaps)), beta 

259 

260 

261def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, 

262 scale=True, nyq=None, fs=None): 

263 """ 

264 FIR filter design using the window method. 

265 

266 This function computes the coefficients of a finite impulse response 

267 filter. The filter will have linear phase; it will be Type I if 

268 `numtaps` is odd and Type II if `numtaps` is even. 

269 

270 Type II filters always have zero response at the Nyquist frequency, so a 

271 ValueError exception is raised if firwin is called with `numtaps` even and 

272 having a passband whose right end is at the Nyquist frequency. 

273 

274 Parameters 

275 ---------- 

276 numtaps : int 

277 Length of the filter (number of coefficients, i.e. the filter 

278 order + 1). `numtaps` must be odd if a passband includes the 

279 Nyquist frequency. 

280 cutoff : float or 1-D array_like 

281 Cutoff frequency of filter (expressed in the same units as `fs`) 

282 OR an array of cutoff frequencies (that is, band edges). In the 

283 latter case, the frequencies in `cutoff` should be positive and 

284 monotonically increasing between 0 and `fs/2`. The values 0 and 

285 `fs/2` must not be included in `cutoff`. 

286 width : float or None, optional 

287 If `width` is not None, then assume it is the approximate width 

288 of the transition region (expressed in the same units as `fs`) 

289 for use in Kaiser FIR filter design. In this case, the `window` 

290 argument is ignored. 

291 window : string or tuple of string and parameter values, optional 

292 Desired window to use. See `scipy.signal.get_window` for a list 

293 of windows and required parameters. 

294 pass_zero : {True, False, 'bandpass', 'lowpass', 'highpass', 'bandstop'}, optional 

295 If True, the gain at the frequency 0 (i.e., the "DC gain") is 1. 

296 If False, the DC gain is 0. Can also be a string argument for the 

297 desired filter type (equivalent to ``btype`` in IIR design functions). 

298 

299 .. versionadded:: 1.3.0 

300 Support for string arguments. 

301 scale : bool, optional 

302 Set to True to scale the coefficients so that the frequency 

303 response is exactly unity at a certain frequency. 

304 That frequency is either: 

305 

306 - 0 (DC) if the first passband starts at 0 (i.e. pass_zero 

307 is True) 

308 - `fs/2` (the Nyquist frequency) if the first passband ends at 

309 `fs/2` (i.e the filter is a single band highpass filter); 

310 center of first passband otherwise 

311 

312 nyq : float, optional 

313 *Deprecated. Use `fs` instead.* This is the Nyquist frequency. 

314 Each frequency in `cutoff` must be between 0 and `nyq`. Default 

315 is 1. 

316 fs : float, optional 

317 The sampling frequency of the signal. Each frequency in `cutoff` 

318 must be between 0 and ``fs/2``. Default is 2. 

319 

320 Returns 

321 ------- 

322 h : (numtaps,) ndarray 

323 Coefficients of length `numtaps` FIR filter. 

324 

325 Raises 

326 ------ 

327 ValueError 

328 If any value in `cutoff` is less than or equal to 0 or greater 

329 than or equal to ``fs/2``, if the values in `cutoff` are not strictly 

330 monotonically increasing, or if `numtaps` is even but a passband 

331 includes the Nyquist frequency. 

332 

333 See Also 

334 -------- 

335 firwin2 

336 firls 

337 minimum_phase 

338 remez 

339 

340 Examples 

341 -------- 

342 Low-pass from 0 to f: 

343 

344 >>> from scipy import signal 

345 >>> numtaps = 3 

346 >>> f = 0.1 

347 >>> signal.firwin(numtaps, f) 

348 array([ 0.06799017, 0.86401967, 0.06799017]) 

349 

350 Use a specific window function: 

351 

352 >>> signal.firwin(numtaps, f, window='nuttall') 

353 array([ 3.56607041e-04, 9.99286786e-01, 3.56607041e-04]) 

354 

355 High-pass ('stop' from 0 to f): 

356 

357 >>> signal.firwin(numtaps, f, pass_zero=False) 

358 array([-0.00859313, 0.98281375, -0.00859313]) 

359 

360 Band-pass: 

361 

362 >>> f1, f2 = 0.1, 0.2 

363 >>> signal.firwin(numtaps, [f1, f2], pass_zero=False) 

364 array([ 0.06301614, 0.88770441, 0.06301614]) 

365 

366 Band-stop: 

367 

368 >>> signal.firwin(numtaps, [f1, f2]) 

369 array([-0.00801395, 1.0160279 , -0.00801395]) 

370 

371 Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1]): 

372 

373 >>> f3, f4 = 0.3, 0.4 

374 >>> signal.firwin(numtaps, [f1, f2, f3, f4]) 

375 array([-0.01376344, 1.02752689, -0.01376344]) 

376 

377 Multi-band (passbands are [f1, f2] and [f3,f4]): 

378 

379 >>> signal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False) 

380 array([ 0.04890915, 0.91284326, 0.04890915]) 

381 

382 """ # noqa: E501 

383 # The major enhancements to this function added in November 2010 were 

384 # developed by Tom Krauss (see ticket #902). 

385 

386 nyq = 0.5 * _get_fs(fs, nyq) 

387 

388 cutoff = np.atleast_1d(cutoff) / float(nyq) 

389 

390 # Check for invalid input. 

391 if cutoff.ndim > 1: 

392 raise ValueError("The cutoff argument must be at most " 

393 "one-dimensional.") 

394 if cutoff.size == 0: 

395 raise ValueError("At least one cutoff frequency must be given.") 

396 if cutoff.min() <= 0 or cutoff.max() >= 1: 

397 raise ValueError("Invalid cutoff frequency: frequencies must be " 

398 "greater than 0 and less than fs/2.") 

399 if np.any(np.diff(cutoff) <= 0): 

400 raise ValueError("Invalid cutoff frequencies: the frequencies " 

401 "must be strictly increasing.") 

402 

403 if width is not None: 

404 # A width was given. Find the beta parameter of the Kaiser window 

405 # and set `window`. This overrides the value of `window` passed in. 

406 atten = kaiser_atten(numtaps, float(width) / nyq) 

407 beta = kaiser_beta(atten) 

408 window = ('kaiser', beta) 

409 

410 if isinstance(pass_zero, str): 

411 if pass_zero in ('bandstop', 'lowpass'): 

412 if pass_zero == 'lowpass': 

413 if cutoff.size != 1: 

414 raise ValueError('cutoff must have one element if ' 

415 'pass_zero=="lowpass", got %s' 

416 % (cutoff.shape,)) 

417 elif cutoff.size <= 1: 

418 raise ValueError('cutoff must have at least two elements if ' 

419 'pass_zero=="bandstop", got %s' 

420 % (cutoff.shape,)) 

421 pass_zero = True 

422 elif pass_zero in ('bandpass', 'highpass'): 

423 if pass_zero == 'highpass': 

424 if cutoff.size != 1: 

425 raise ValueError('cutoff must have one element if ' 

426 'pass_zero=="highpass", got %s' 

427 % (cutoff.shape,)) 

428 elif cutoff.size <= 1: 

429 raise ValueError('cutoff must have at least two elements if ' 

430 'pass_zero=="bandpass", got %s' 

431 % (cutoff.shape,)) 

432 pass_zero = False 

433 else: 

434 raise ValueError('pass_zero must be True, False, "bandpass", ' 

435 '"lowpass", "highpass", or "bandstop", got ' 

436 '%s' % (pass_zero,)) 

437 pass_zero = bool(operator.index(pass_zero)) # ensure bool-like 

438 

439 pass_nyquist = bool(cutoff.size & 1) ^ pass_zero 

440 if pass_nyquist and numtaps % 2 == 0: 

441 raise ValueError("A filter with an even number of coefficients must " 

442 "have zero response at the Nyquist frequency.") 

443 

444 # Insert 0 and/or 1 at the ends of cutoff so that the length of cutoff 

445 # is even, and each pair in cutoff corresponds to passband. 

446 cutoff = np.hstack(([0.0] * pass_zero, cutoff, [1.0] * pass_nyquist)) 

447 

448 # `bands` is a 2-D array; each row gives the left and right edges of 

449 # a passband. 

450 bands = cutoff.reshape(-1, 2) 

451 

452 # Build up the coefficients. 

453 alpha = 0.5 * (numtaps - 1) 

454 m = np.arange(0, numtaps) - alpha 

455 h = 0 

456 for left, right in bands: 

457 h += right * sinc(right * m) 

458 h -= left * sinc(left * m) 

459 

460 # Get and apply the window function. 

461 from .signaltools import get_window 

462 win = get_window(window, numtaps, fftbins=False) 

463 h *= win 

464 

465 # Now handle scaling if desired. 

466 if scale: 

467 # Get the first passband. 

468 left, right = bands[0] 

469 if left == 0: 

470 scale_frequency = 0.0 

471 elif right == 1: 

472 scale_frequency = 1.0 

473 else: 

474 scale_frequency = 0.5 * (left + right) 

475 c = np.cos(np.pi * m * scale_frequency) 

476 s = np.sum(h * c) 

477 h /= s 

478 

479 return h 

480 

481 

482# Original version of firwin2 from scipy ticket #457, submitted by "tash". 

483# 

484# Rewritten by Warren Weckesser, 2010. 

485 

486def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None, 

487 antisymmetric=False, fs=None): 

488 """ 

489 FIR filter design using the window method. 

490 

491 From the given frequencies `freq` and corresponding gains `gain`, 

492 this function constructs an FIR filter with linear phase and 

493 (approximately) the given frequency response. 

494 

495 Parameters 

496 ---------- 

497 numtaps : int 

498 The number of taps in the FIR filter. `numtaps` must be less than 

499 `nfreqs`. 

500 freq : array_like, 1-D 

501 The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being 

502 Nyquist. The Nyquist frequency is half `fs`. 

503 The values in `freq` must be nondecreasing. A value can be repeated 

504 once to implement a discontinuity. The first value in `freq` must 

505 be 0, and the last value must be ``fs/2``. Values 0 and ``fs/2`` must 

506 not be repeated. 

507 gain : array_like 

508 The filter gains at the frequency sampling points. Certain 

509 constraints to gain values, depending on the filter type, are applied, 

510 see Notes for details. 

511 nfreqs : int, optional 

512 The size of the interpolation mesh used to construct the filter. 

513 For most efficient behavior, this should be a power of 2 plus 1 

514 (e.g, 129, 257, etc). The default is one more than the smallest 

515 power of 2 that is not less than `numtaps`. `nfreqs` must be greater 

516 than `numtaps`. 

517 window : string or (string, float) or float, or None, optional 

518 Window function to use. Default is "hamming". See 

519 `scipy.signal.get_window` for the complete list of possible values. 

520 If None, no window function is applied. 

521 nyq : float, optional 

522 *Deprecated. Use `fs` instead.* This is the Nyquist frequency. 

523 Each frequency in `freq` must be between 0 and `nyq`. Default is 1. 

524 antisymmetric : bool, optional 

525 Whether resulting impulse response is symmetric/antisymmetric. 

526 See Notes for more details. 

527 fs : float, optional 

528 The sampling frequency of the signal. Each frequency in `cutoff` 

529 must be between 0 and ``fs/2``. Default is 2. 

530 

531 Returns 

532 ------- 

533 taps : ndarray 

534 The filter coefficients of the FIR filter, as a 1-D array of length 

535 `numtaps`. 

536 

537 See also 

538 -------- 

539 firls 

540 firwin 

541 minimum_phase 

542 remez 

543 

544 Notes 

545 ----- 

546 From the given set of frequencies and gains, the desired response is 

547 constructed in the frequency domain. The inverse FFT is applied to the 

548 desired response to create the associated convolution kernel, and the 

549 first `numtaps` coefficients of this kernel, scaled by `window`, are 

550 returned. 

551 

552 The FIR filter will have linear phase. The type of filter is determined by 

553 the value of 'numtaps` and `antisymmetric` flag. 

554 There are four possible combinations: 

555 

556 - odd `numtaps`, `antisymmetric` is False, type I filter is produced 

557 - even `numtaps`, `antisymmetric` is False, type II filter is produced 

558 - odd `numtaps`, `antisymmetric` is True, type III filter is produced 

559 - even `numtaps`, `antisymmetric` is True, type IV filter is produced 

560 

561 Magnitude response of all but type I filters are subjects to following 

562 constraints: 

563 

564 - type II -- zero at the Nyquist frequency 

565 - type III -- zero at zero and Nyquist frequencies 

566 - type IV -- zero at zero frequency 

567 

568 .. versionadded:: 0.9.0 

569 

570 References 

571 ---------- 

572 .. [1] Oppenheim, A. V. and Schafer, R. W., "Discrete-Time Signal 

573 Processing", Prentice-Hall, Englewood Cliffs, New Jersey (1989). 

574 (See, for example, Section 7.4.) 

575 

576 .. [2] Smith, Steven W., "The Scientist and Engineer's Guide to Digital 

577 Signal Processing", Ch. 17. http://www.dspguide.com/ch17/1.htm 

578 

579 Examples 

580 -------- 

581 A lowpass FIR filter with a response that is 1 on [0.0, 0.5], and 

582 that decreases linearly on [0.5, 1.0] from 1 to 0: 

583 

584 >>> from scipy import signal 

585 >>> taps = signal.firwin2(150, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0]) 

586 >>> print(taps[72:78]) 

587 [-0.02286961 -0.06362756 0.57310236 0.57310236 -0.06362756 -0.02286961] 

588 

589 """ 

590 nyq = 0.5 * _get_fs(fs, nyq) 

591 

592 if len(freq) != len(gain): 

593 raise ValueError('freq and gain must be of same length.') 

594 

595 if nfreqs is not None and numtaps >= nfreqs: 

596 raise ValueError(('ntaps must be less than nfreqs, but firwin2 was ' 

597 'called with ntaps=%d and nfreqs=%s') % 

598 (numtaps, nfreqs)) 

599 

600 if freq[0] != 0 or freq[-1] != nyq: 

601 raise ValueError('freq must start with 0 and end with fs/2.') 

602 d = np.diff(freq) 

603 if (d < 0).any(): 

604 raise ValueError('The values in freq must be nondecreasing.') 

605 d2 = d[:-1] + d[1:] 

606 if (d2 == 0).any(): 

607 raise ValueError('A value in freq must not occur more than twice.') 

608 if freq[1] == 0: 

609 raise ValueError('Value 0 must not be repeated in freq') 

610 if freq[-2] == nyq: 

611 raise ValueError('Value fs/2 must not be repeated in freq') 

612 

613 if antisymmetric: 

614 if numtaps % 2 == 0: 

615 ftype = 4 

616 else: 

617 ftype = 3 

618 else: 

619 if numtaps % 2 == 0: 

620 ftype = 2 

621 else: 

622 ftype = 1 

623 

624 if ftype == 2 and gain[-1] != 0.0: 

625 raise ValueError("A Type II filter must have zero gain at the " 

626 "Nyquist frequency.") 

627 elif ftype == 3 and (gain[0] != 0.0 or gain[-1] != 0.0): 

628 raise ValueError("A Type III filter must have zero gain at zero " 

629 "and Nyquist frequencies.") 

630 elif ftype == 4 and gain[0] != 0.0: 

631 raise ValueError("A Type IV filter must have zero gain at zero " 

632 "frequency.") 

633 

634 if nfreqs is None: 

635 nfreqs = 1 + 2 ** int(ceil(log(numtaps, 2))) 

636 

637 if (d == 0).any(): 

638 # Tweak any repeated values in freq so that interp works. 

639 freq = np.array(freq, copy=True) 

640 eps = np.finfo(float).eps * nyq 

641 for k in range(len(freq) - 1): 

642 if freq[k] == freq[k + 1]: 

643 freq[k] = freq[k] - eps 

644 freq[k + 1] = freq[k + 1] + eps 

645 # Check if freq is strictly increasing after tweak 

646 d = np.diff(freq) 

647 if (d <= 0).any(): 

648 raise ValueError("freq cannot contain numbers that are too close " 

649 "(within eps * (fs/2): " 

650 "{}) to a repeated value".format(eps)) 

651 

652 # Linearly interpolate the desired response on a uniform mesh `x`. 

653 x = np.linspace(0.0, nyq, nfreqs) 

654 fx = np.interp(x, freq, gain) 

655 

656 # Adjust the phases of the coefficients so that the first `ntaps` of the 

657 # inverse FFT are the desired filter coefficients. 

658 shift = np.exp(-(numtaps - 1) / 2. * 1.j * np.pi * x / nyq) 

659 if ftype > 2: 

660 shift *= 1j 

661 

662 fx2 = fx * shift 

663 

664 # Use irfft to compute the inverse FFT. 

665 out_full = irfft(fx2) 

666 

667 if window is not None: 

668 # Create the window to apply to the filter coefficients. 

669 from .signaltools import get_window 

670 wind = get_window(window, numtaps, fftbins=False) 

671 else: 

672 wind = 1 

673 

674 # Keep only the first `numtaps` coefficients in `out`, and multiply by 

675 # the window. 

676 out = out_full[:numtaps] * wind 

677 

678 if ftype == 3: 

679 out[out.size // 2] = 0.0 

680 

681 return out 

682 

683 

684def remez(numtaps, bands, desired, weight=None, Hz=None, type='bandpass', 

685 maxiter=25, grid_density=16, fs=None): 

686 """ 

687 Calculate the minimax optimal filter using the Remez exchange algorithm. 

688 

689 Calculate the filter-coefficients for the finite impulse response 

690 (FIR) filter whose transfer function minimizes the maximum error 

691 between the desired gain and the realized gain in the specified 

692 frequency bands using the Remez exchange algorithm. 

693 

694 Parameters 

695 ---------- 

696 numtaps : int 

697 The desired number of taps in the filter. The number of taps is 

698 the number of terms in the filter, or the filter order plus one. 

699 bands : array_like 

700 A monotonic sequence containing the band edges. 

701 All elements must be non-negative and less than half the sampling 

702 frequency as given by `fs`. 

703 desired : array_like 

704 A sequence half the size of bands containing the desired gain 

705 in each of the specified bands. 

706 weight : array_like, optional 

707 A relative weighting to give to each band region. The length of 

708 `weight` has to be half the length of `bands`. 

709 Hz : scalar, optional 

710 *Deprecated. Use `fs` instead.* 

711 The sampling frequency in Hz. Default is 1. 

712 type : {'bandpass', 'differentiator', 'hilbert'}, optional 

713 The type of filter: 

714 

715 * 'bandpass' : flat response in bands. This is the default. 

716 

717 * 'differentiator' : frequency proportional response in bands. 

718 

719 * 'hilbert' : filter with odd symmetry, that is, type III 

720 (for even order) or type IV (for odd order) 

721 linear phase filters. 

722 

723 maxiter : int, optional 

724 Maximum number of iterations of the algorithm. Default is 25. 

725 grid_density : int, optional 

726 Grid density. The dense grid used in `remez` is of size 

727 ``(numtaps + 1) * grid_density``. Default is 16. 

728 fs : float, optional 

729 The sampling frequency of the signal. Default is 1. 

730 

731 Returns 

732 ------- 

733 out : ndarray 

734 A rank-1 array containing the coefficients of the optimal 

735 (in a minimax sense) filter. 

736 

737 See Also 

738 -------- 

739 firls 

740 firwin 

741 firwin2 

742 minimum_phase 

743 

744 References 

745 ---------- 

746 .. [1] J. H. McClellan and T. W. Parks, "A unified approach to the 

747 design of optimum FIR linear phase digital filters", 

748 IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973. 

749 .. [2] J. H. McClellan, T. W. Parks and L. R. Rabiner, "A Computer 

750 Program for Designing Optimum FIR Linear Phase Digital 

751 Filters", IEEE Trans. Audio Electroacoust., vol. AU-21, 

752 pp. 506-525, 1973. 

753 

754 Examples 

755 -------- 

756 In these examples `remez` gets used creating a bandpass, bandstop, lowpass 

757 and highpass filter. The used parameters are the filter order, an array 

758 with according frequency boundaries, the desired attenuation values and the 

759 sampling frequency. Using `freqz` the corresponding frequency response 

760 gets calculated and plotted. 

761 

762 >>> from scipy import signal 

763 >>> import matplotlib.pyplot as plt 

764 

765 >>> def plot_response(fs, w, h, title): 

766 ... "Utility function to plot response functions" 

767 ... fig = plt.figure() 

768 ... ax = fig.add_subplot(111) 

769 ... ax.plot(0.5*fs*w/np.pi, 20*np.log10(np.abs(h))) 

770 ... ax.set_ylim(-40, 5) 

771 ... ax.set_xlim(0, 0.5*fs) 

772 ... ax.grid(True) 

773 ... ax.set_xlabel('Frequency (Hz)') 

774 ... ax.set_ylabel('Gain (dB)') 

775 ... ax.set_title(title) 

776 

777 This example shows a steep low pass transition according to the small 

778 transition width and high filter order: 

779 

780 >>> fs = 22050.0 # Sample rate, Hz 

781 >>> cutoff = 8000.0 # Desired cutoff frequency, Hz 

782 >>> trans_width = 100 # Width of transition from pass band to stop band, Hz 

783 >>> numtaps = 400 # Size of the FIR filter. 

784 >>> taps = signal.remez(numtaps, [0, cutoff, cutoff + trans_width, 0.5*fs], [1, 0], Hz=fs) 

785 >>> w, h = signal.freqz(taps, [1], worN=2000) 

786 >>> plot_response(fs, w, h, "Low-pass Filter") 

787 

788 This example shows a high pass filter: 

789 

790 >>> fs = 22050.0 # Sample rate, Hz 

791 >>> cutoff = 2000.0 # Desired cutoff frequency, Hz 

792 >>> trans_width = 250 # Width of transition from pass band to stop band, Hz 

793 >>> numtaps = 125 # Size of the FIR filter. 

794 >>> taps = signal.remez(numtaps, [0, cutoff - trans_width, cutoff, 0.5*fs], 

795 ... [0, 1], Hz=fs) 

796 >>> w, h = signal.freqz(taps, [1], worN=2000) 

797 >>> plot_response(fs, w, h, "High-pass Filter") 

798 

799 For a signal sampled with 22 kHz a bandpass filter with a pass band of 2-5 

800 kHz gets calculated using the Remez algorithm. The transition width is 260 

801 Hz and the filter order 10: 

802 

803 >>> fs = 22000.0 # Sample rate, Hz 

804 >>> band = [2000, 5000] # Desired pass band, Hz 

805 >>> trans_width = 260 # Width of transition from pass band to stop band, Hz 

806 >>> numtaps = 10 # Size of the FIR filter. 

807 >>> edges = [0, band[0] - trans_width, band[0], band[1], 

808 ... band[1] + trans_width, 0.5*fs] 

809 >>> taps = signal.remez(numtaps, edges, [0, 1, 0], Hz=fs) 

810 >>> w, h = signal.freqz(taps, [1], worN=2000) 

811 >>> plot_response(fs, w, h, "Band-pass Filter") 

812 

813 It can be seen that for this bandpass filter, the low order leads to higher 

814 ripple and less steep transitions. There is very low attenuation in the 

815 stop band and little overshoot in the pass band. Of course the desired 

816 gain can be better approximated with a higher filter order. 

817 

818 The next example shows a bandstop filter. Because of the high filter order 

819 the transition is quite steep: 

820 

821 >>> fs = 20000.0 # Sample rate, Hz 

822 >>> band = [6000, 8000] # Desired stop band, Hz 

823 >>> trans_width = 200 # Width of transition from pass band to stop band, Hz 

824 >>> numtaps = 175 # Size of the FIR filter. 

825 >>> edges = [0, band[0] - trans_width, band[0], band[1], band[1] + trans_width, 0.5*fs] 

826 >>> taps = signal.remez(numtaps, edges, [1, 0, 1], Hz=fs) 

827 >>> w, h = signal.freqz(taps, [1], worN=2000) 

828 >>> plot_response(fs, w, h, "Band-stop Filter") 

829 

830 >>> plt.show() 

831 

832 """ 

833 if Hz is None and fs is None: 

834 fs = 1.0 

835 elif Hz is not None: 

836 if fs is not None: 

837 raise ValueError("Values cannot be given for both 'Hz' and 'fs'.") 

838 fs = Hz 

839 

840 # Convert type 

841 try: 

842 tnum = {'bandpass': 1, 'differentiator': 2, 'hilbert': 3}[type] 

843 except KeyError: 

844 raise ValueError("Type must be 'bandpass', 'differentiator', " 

845 "or 'hilbert'") 

846 

847 # Convert weight 

848 if weight is None: 

849 weight = [1] * len(desired) 

850 

851 bands = np.asarray(bands).copy() 

852 return sigtools._remez(numtaps, bands, desired, weight, tnum, fs, 

853 maxiter, grid_density) 

854 

855 

856def firls(numtaps, bands, desired, weight=None, nyq=None, fs=None): 

857 """ 

858 FIR filter design using least-squares error minimization. 

859 

860 Calculate the filter coefficients for the linear-phase finite 

861 impulse response (FIR) filter which has the best approximation 

862 to the desired frequency response described by `bands` and 

863 `desired` in the least squares sense (i.e., the integral of the 

864 weighted mean-squared error within the specified bands is 

865 minimized). 

866 

867 Parameters 

868 ---------- 

869 numtaps : int 

870 The number of taps in the FIR filter. `numtaps` must be odd. 

871 bands : array_like 

872 A monotonic nondecreasing sequence containing the band edges in 

873 Hz. All elements must be non-negative and less than or equal to 

874 the Nyquist frequency given by `nyq`. 

875 desired : array_like 

876 A sequence the same size as `bands` containing the desired gain 

877 at the start and end point of each band. 

878 weight : array_like, optional 

879 A relative weighting to give to each band region when solving 

880 the least squares problem. `weight` has to be half the size of 

881 `bands`. 

882 nyq : float, optional 

883 *Deprecated. Use `fs` instead.* 

884 Nyquist frequency. Each frequency in `bands` must be between 0 

885 and `nyq` (inclusive). Default is 1. 

886 fs : float, optional 

887 The sampling frequency of the signal. Each frequency in `bands` 

888 must be between 0 and ``fs/2`` (inclusive). Default is 2. 

889 

890 Returns 

891 ------- 

892 coeffs : ndarray 

893 Coefficients of the optimal (in a least squares sense) FIR filter. 

894 

895 See also 

896 -------- 

897 firwin 

898 firwin2 

899 minimum_phase 

900 remez 

901 

902 Notes 

903 ----- 

904 This implementation follows the algorithm given in [1]_. 

905 As noted there, least squares design has multiple advantages: 

906 

907 1. Optimal in a least-squares sense. 

908 2. Simple, non-iterative method. 

909 3. The general solution can obtained by solving a linear 

910 system of equations. 

911 4. Allows the use of a frequency dependent weighting function. 

912 

913 This function constructs a Type I linear phase FIR filter, which 

914 contains an odd number of `coeffs` satisfying for :math:`n < numtaps`: 

915 

916 .. math:: coeffs(n) = coeffs(numtaps - 1 - n) 

917 

918 The odd number of coefficients and filter symmetry avoid boundary 

919 conditions that could otherwise occur at the Nyquist and 0 frequencies 

920 (e.g., for Type II, III, or IV variants). 

921 

922 .. versionadded:: 0.18 

923 

924 References 

925 ---------- 

926 .. [1] Ivan Selesnick, Linear-Phase Fir Filter Design By Least Squares. 

927 OpenStax CNX. Aug 9, 2005. 

928 http://cnx.org/contents/eb1ecb35-03a9-4610-ba87-41cd771c95f2@7 

929 

930 Examples 

931 -------- 

932 We want to construct a band-pass filter. Note that the behavior in the 

933 frequency ranges between our stop bands and pass bands is unspecified, 

934 and thus may overshoot depending on the parameters of our filter: 

935 

936 >>> from scipy import signal 

937 >>> import matplotlib.pyplot as plt 

938 >>> fig, axs = plt.subplots(2) 

939 >>> fs = 10.0 # Hz 

940 >>> desired = (0, 0, 1, 1, 0, 0) 

941 >>> for bi, bands in enumerate(((0, 1, 2, 3, 4, 5), (0, 1, 2, 4, 4.5, 5))): 

942 ... fir_firls = signal.firls(73, bands, desired, fs=fs) 

943 ... fir_remez = signal.remez(73, bands, desired[::2], fs=fs) 

944 ... fir_firwin2 = signal.firwin2(73, bands, desired, fs=fs) 

945 ... hs = list() 

946 ... ax = axs[bi] 

947 ... for fir in (fir_firls, fir_remez, fir_firwin2): 

948 ... freq, response = signal.freqz(fir) 

949 ... hs.append(ax.semilogy(0.5*fs*freq/np.pi, np.abs(response))[0]) 

950 ... for band, gains in zip(zip(bands[::2], bands[1::2]), 

951 ... zip(desired[::2], desired[1::2])): 

952 ... ax.semilogy(band, np.maximum(gains, 1e-7), 'k--', linewidth=2) 

953 ... if bi == 0: 

954 ... ax.legend(hs, ('firls', 'remez', 'firwin2'), 

955 ... loc='lower center', frameon=False) 

956 ... else: 

957 ... ax.set_xlabel('Frequency (Hz)') 

958 ... ax.grid(True) 

959 ... ax.set(title='Band-pass %d-%d Hz' % bands[2:4], ylabel='Magnitude') 

960 ... 

961 >>> fig.tight_layout() 

962 >>> plt.show() 

963 

964 """ # noqa 

965 nyq = 0.5 * _get_fs(fs, nyq) 

966 

967 numtaps = int(numtaps) 

968 if numtaps % 2 == 0 or numtaps < 1: 

969 raise ValueError("numtaps must be odd and >= 1") 

970 M = (numtaps-1) // 2 

971 

972 # normalize bands 0->1 and make it 2 columns 

973 nyq = float(nyq) 

974 if nyq <= 0: 

975 raise ValueError('nyq must be positive, got %s <= 0.' % nyq) 

976 bands = np.asarray(bands).flatten() / nyq 

977 if len(bands) % 2 != 0: 

978 raise ValueError("bands must contain frequency pairs.") 

979 if (bands < 0).any() or (bands > 1).any(): 

980 raise ValueError("bands must be between 0 and 1 relative to Nyquist") 

981 bands.shape = (-1, 2) 

982 

983 # check remaining params 

984 desired = np.asarray(desired).flatten() 

985 if bands.size != desired.size: 

986 raise ValueError("desired must have one entry per frequency, got %s " 

987 "gains for %s frequencies." 

988 % (desired.size, bands.size)) 

989 desired.shape = (-1, 2) 

990 if (np.diff(bands) <= 0).any() or (np.diff(bands[:, 0]) < 0).any(): 

991 raise ValueError("bands must be monotonically nondecreasing and have " 

992 "width > 0.") 

993 if (bands[:-1, 1] > bands[1:, 0]).any(): 

994 raise ValueError("bands must not overlap.") 

995 if (desired < 0).any(): 

996 raise ValueError("desired must be non-negative.") 

997 if weight is None: 

998 weight = np.ones(len(desired)) 

999 weight = np.asarray(weight).flatten() 

1000 if len(weight) != len(desired): 

1001 raise ValueError("weight must be the same size as the number of " 

1002 "band pairs (%s)." % (len(bands),)) 

1003 if (weight < 0).any(): 

1004 raise ValueError("weight must be non-negative.") 

1005 

1006 # Set up the linear matrix equation to be solved, Qa = b 

1007 

1008 # We can express Q(k,n) = 0.5 Q1(k,n) + 0.5 Q2(k,n) 

1009 # where Q1(k,n)=q(k−n) and Q2(k,n)=q(k+n), i.e. a Toeplitz plus Hankel. 

1010 

1011 # We omit the factor of 0.5 above, instead adding it during coefficient 

1012 # calculation. 

1013 

1014 # We also omit the 1/π from both Q and b equations, as they cancel 

1015 # during solving. 

1016 

1017 # We have that: 

1018 # q(n) = 1/π ∫W(ω)cos(nω)dω (over 0->π) 

1019 # Using our nomalization ω=πf and with a constant weight W over each 

1020 # interval f1->f2 we get: 

1021 # q(n) = W∫cos(πnf)df (0->1) = Wf sin(πnf)/πnf 

1022 # integrated over each f1->f2 pair (i.e., value at f2 - value at f1). 

1023 n = np.arange(numtaps)[:, np.newaxis, np.newaxis] 

1024 q = np.dot(np.diff(np.sinc(bands * n) * bands, axis=2)[:, :, 0], weight) 

1025 

1026 # Now we assemble our sum of Toeplitz and Hankel 

1027 Q1 = toeplitz(q[:M+1]) 

1028 Q2 = hankel(q[:M+1], q[M:]) 

1029 Q = Q1 + Q2 

1030 

1031 # Now for b(n) we have that: 

1032 # b(n) = 1/π ∫ W(ω)D(ω)cos(nω)dω (over 0->π) 

1033 # Using our normalization ω=πf and with a constant weight W over each 

1034 # interval and a linear term for D(ω) we get (over each f1->f2 interval): 

1035 # b(n) = W ∫ (mf+c)cos(πnf)df 

1036 # = f(mf+c)sin(πnf)/πnf + mf**2 cos(nπf)/(πnf)**2 

1037 # integrated over each f1->f2 pair (i.e., value at f2 - value at f1). 

1038 n = n[:M + 1] # only need this many coefficients here 

1039 # Choose m and c such that we are at the start and end weights 

1040 m = (np.diff(desired, axis=1) / np.diff(bands, axis=1)) 

1041 c = desired[:, [0]] - bands[:, [0]] * m 

1042 b = bands * (m*bands + c) * np.sinc(bands * n) 

1043 # Use L'Hospital's rule here for cos(nπf)/(πnf)**2 @ n=0 

1044 b[0] -= m * bands * bands / 2. 

1045 b[1:] += m * np.cos(n[1:] * np.pi * bands) / (np.pi * n[1:]) ** 2 

1046 b = np.dot(np.diff(b, axis=2)[:, :, 0], weight) 

1047 

1048 # Now we can solve the equation 

1049 try: # try the fast way 

1050 with warnings.catch_warnings(record=True) as w: 

1051 warnings.simplefilter('always') 

1052 a = solve(Q, b, sym_pos=True, check_finite=False) 

1053 for ww in w: 

1054 if (ww.category == LinAlgWarning and 

1055 str(ww.message).startswith('Ill-conditioned matrix')): 

1056 raise LinAlgError(str(ww.message)) 

1057 except LinAlgError: # in case Q is rank deficient 

1058 # This is faster than pinvh, even though we don't explicitly use 

1059 # the symmetry here. gelsy was faster than gelsd and gelss in 

1060 # some non-exhaustive tests. 

1061 a = lstsq(Q, b, lapack_driver='gelsy')[0] 

1062 

1063 # make coefficients symmetric (linear phase) 

1064 coeffs = np.hstack((a[:0:-1], 2 * a[0], a[1:])) 

1065 return coeffs 

1066 

1067 

1068def _dhtm(mag): 

1069 """Compute the modified 1-D discrete Hilbert transform 

1070 

1071 Parameters 

1072 ---------- 

1073 mag : ndarray 

1074 The magnitude spectrum. Should be 1-D with an even length, and 

1075 preferably a fast length for FFT/IFFT. 

1076 """ 

1077 # Adapted based on code by Niranjan Damera-Venkata, 

1078 # Brian L. Evans and Shawn R. McCaslin (see refs for `minimum_phase`) 

1079 sig = np.zeros(len(mag)) 

1080 # Leave Nyquist and DC at 0, knowing np.abs(fftfreq(N)[midpt]) == 0.5 

1081 midpt = len(mag) // 2 

1082 sig[1:midpt] = 1 

1083 sig[midpt+1:] = -1 

1084 # eventually if we want to support complex filters, we will need a 

1085 # np.abs() on the mag inside the log, and should remove the .real 

1086 recon = ifft(mag * np.exp(fft(sig * ifft(np.log(mag))))).real 

1087 return recon 

1088 

1089 

1090def minimum_phase(h, method='homomorphic', n_fft=None): 

1091 """Convert a linear-phase FIR filter to minimum phase 

1092 

1093 Parameters 

1094 ---------- 

1095 h : array 

1096 Linear-phase FIR filter coefficients. 

1097 method : {'hilbert', 'homomorphic'} 

1098 The method to use: 

1099 

1100 'homomorphic' (default) 

1101 This method [4]_ [5]_ works best with filters with an 

1102 odd number of taps, and the resulting minimum phase filter 

1103 will have a magnitude response that approximates the square 

1104 root of the the original filter's magnitude response. 

1105 

1106 'hilbert' 

1107 This method [1]_ is designed to be used with equiripple 

1108 filters (e.g., from `remez`) with unity or zero gain 

1109 regions. 

1110 

1111 n_fft : int 

1112 The number of points to use for the FFT. Should be at least a 

1113 few times larger than the signal length (see Notes). 

1114 

1115 Returns 

1116 ------- 

1117 h_minimum : array 

1118 The minimum-phase version of the filter, with length 

1119 ``(length(h) + 1) // 2``. 

1120 

1121 See Also 

1122 -------- 

1123 firwin 

1124 firwin2 

1125 remez 

1126 

1127 Notes 

1128 ----- 

1129 Both the Hilbert [1]_ or homomorphic [4]_ [5]_ methods require selection 

1130 of an FFT length to estimate the complex cepstrum of the filter. 

1131 

1132 In the case of the Hilbert method, the deviation from the ideal 

1133 spectrum ``epsilon`` is related to the number of stopband zeros 

1134 ``n_stop`` and FFT length ``n_fft`` as:: 

1135 

1136 epsilon = 2. * n_stop / n_fft 

1137 

1138 For example, with 100 stopband zeros and a FFT length of 2048, 

1139 ``epsilon = 0.0976``. If we conservatively assume that the number of 

1140 stopband zeros is one less than the filter length, we can take the FFT 

1141 length to be the next power of 2 that satisfies ``epsilon=0.01`` as:: 

1142 

1143 n_fft = 2 ** int(np.ceil(np.log2(2 * (len(h) - 1) / 0.01))) 

1144 

1145 This gives reasonable results for both the Hilbert and homomorphic 

1146 methods, and gives the value used when ``n_fft=None``. 

1147 

1148 Alternative implementations exist for creating minimum-phase filters, 

1149 including zero inversion [2]_ and spectral factorization [3]_ [4]_. 

1150 For more information, see: 

1151 

1152 http://dspguru.com/dsp/howtos/how-to-design-minimum-phase-fir-filters 

1153 

1154 Examples 

1155 -------- 

1156 Create an optimal linear-phase filter, then convert it to minimum phase: 

1157 

1158 >>> from scipy.signal import remez, minimum_phase, freqz, group_delay 

1159 >>> import matplotlib.pyplot as plt 

1160 >>> freq = [0, 0.2, 0.3, 1.0] 

1161 >>> desired = [1, 0] 

1162 >>> h_linear = remez(151, freq, desired, Hz=2.) 

1163 

1164 Convert it to minimum phase: 

1165 

1166 >>> h_min_hom = minimum_phase(h_linear, method='homomorphic') 

1167 >>> h_min_hil = minimum_phase(h_linear, method='hilbert') 

1168 

1169 Compare the three filters: 

1170 

1171 >>> fig, axs = plt.subplots(4, figsize=(4, 8)) 

1172 >>> for h, style, color in zip((h_linear, h_min_hom, h_min_hil), 

1173 ... ('-', '-', '--'), ('k', 'r', 'c')): 

1174 ... w, H = freqz(h) 

1175 ... w, gd = group_delay((h, 1)) 

1176 ... w /= np.pi 

1177 ... axs[0].plot(h, color=color, linestyle=style) 

1178 ... axs[1].plot(w, np.abs(H), color=color, linestyle=style) 

1179 ... axs[2].plot(w, 20 * np.log10(np.abs(H)), color=color, linestyle=style) 

1180 ... axs[3].plot(w, gd, color=color, linestyle=style) 

1181 >>> for ax in axs: 

1182 ... ax.grid(True, color='0.5') 

1183 ... ax.fill_between(freq[1:3], *ax.get_ylim(), color='#ffeeaa', zorder=1) 

1184 >>> axs[0].set(xlim=[0, len(h_linear) - 1], ylabel='Amplitude', xlabel='Samples') 

1185 >>> axs[1].legend(['Linear', 'Min-Hom', 'Min-Hil'], title='Phase') 

1186 >>> for ax, ylim in zip(axs[1:], ([0, 1.1], [-150, 10], [-60, 60])): 

1187 ... ax.set(xlim=[0, 1], ylim=ylim, xlabel='Frequency') 

1188 >>> axs[1].set(ylabel='Magnitude') 

1189 >>> axs[2].set(ylabel='Magnitude (dB)') 

1190 >>> axs[3].set(ylabel='Group delay') 

1191 >>> plt.tight_layout() 

1192 

1193 References 

1194 ---------- 

1195 .. [1] N. Damera-Venkata and B. L. Evans, "Optimal design of real and 

1196 complex minimum phase digital FIR filters," Acoustics, Speech, 

1197 and Signal Processing, 1999. Proceedings., 1999 IEEE International 

1198 Conference on, Phoenix, AZ, 1999, pp. 1145-1148 vol.3. 

1199 doi: 10.1109/ICASSP.1999.756179 

1200 .. [2] X. Chen and T. W. Parks, "Design of optimal minimum phase FIR 

1201 filters by direct factorization," Signal Processing, 

1202 vol. 10, no. 4, pp. 369-383, Jun. 1986. 

1203 .. [3] T. Saramaki, "Finite Impulse Response Filter Design," in 

1204 Handbook for Digital Signal Processing, chapter 4, 

1205 New York: Wiley-Interscience, 1993. 

1206 .. [4] J. S. Lim, Advanced Topics in Signal Processing. 

1207 Englewood Cliffs, N.J.: Prentice Hall, 1988. 

1208 .. [5] A. V. Oppenheim, R. W. Schafer, and J. R. Buck, 

1209 "Discrete-Time Signal Processing," 2nd edition. 

1210 Upper Saddle River, N.J.: Prentice Hall, 1999. 

1211 """ # noqa 

1212 h = np.asarray(h) 

1213 if np.iscomplexobj(h): 

1214 raise ValueError('Complex filters not supported') 

1215 if h.ndim != 1 or h.size <= 2: 

1216 raise ValueError('h must be 1-D and at least 2 samples long') 

1217 n_half = len(h) // 2 

1218 if not np.allclose(h[-n_half:][::-1], h[:n_half]): 

1219 warnings.warn('h does not appear to by symmetric, conversion may ' 

1220 'fail', RuntimeWarning) 

1221 if not isinstance(method, str) or method not in \ 

1222 ('homomorphic', 'hilbert',): 

1223 raise ValueError('method must be "homomorphic" or "hilbert", got %r' 

1224 % (method,)) 

1225 if n_fft is None: 

1226 n_fft = 2 ** int(np.ceil(np.log2(2 * (len(h) - 1) / 0.01))) 

1227 n_fft = int(n_fft) 

1228 if n_fft < len(h): 

1229 raise ValueError('n_fft must be at least len(h)==%s' % len(h)) 

1230 if method == 'hilbert': 

1231 w = np.arange(n_fft) * (2 * np.pi / n_fft * n_half) 

1232 H = np.real(fft(h, n_fft) * np.exp(1j * w)) 

1233 dp = max(H) - 1 

1234 ds = 0 - min(H) 

1235 S = 4. / (np.sqrt(1+dp+ds) + np.sqrt(1-dp+ds)) ** 2 

1236 H += ds 

1237 H *= S 

1238 H = np.sqrt(H, out=H) 

1239 H += 1e-10 # ensure that the log does not explode 

1240 h_minimum = _dhtm(H) 

1241 else: # method == 'homomorphic' 

1242 # zero-pad; calculate the DFT 

1243 h_temp = np.abs(fft(h, n_fft)) 

1244 # take 0.25*log(|H|**2) = 0.5*log(|H|) 

1245 h_temp += 1e-7 * h_temp[h_temp > 0].min() # don't let log blow up 

1246 np.log(h_temp, out=h_temp) 

1247 h_temp *= 0.5 

1248 # IDFT 

1249 h_temp = ifft(h_temp).real 

1250 # multiply pointwise by the homomorphic filter 

1251 # lmin[n] = 2u[n] - d[n] 

1252 win = np.zeros(n_fft) 

1253 win[0] = 1 

1254 stop = (len(h) + 1) // 2 

1255 win[1:stop] = 2 

1256 if len(h) % 2: 

1257 win[stop] = 1 

1258 h_temp *= win 

1259 h_temp = ifft(np.exp(fft(h_temp))) 

1260 h_minimum = h_temp.real 

1261 n_out = n_half + len(h) % 2 

1262 return h_minimum[:n_out]