Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/signal/fir_filter_design.py : 7%

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
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)
13from . import sigtools
15__all__ = ['kaiser_beta', 'kaiser_atten', 'kaiserord',
16 'firwin', 'firwin2', 'remez', 'firls', 'minimum_phase']
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
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.
47def kaiser_beta(a):
48 """Compute the Kaiser parameter `beta`, given the attenuation `a`.
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.
56 Returns
57 -------
58 beta : float
59 The `beta` parameter to be used in the formula for a Kaiser window.
61 References
62 ----------
63 Oppenheim, Schafer, "Discrete-Time Signal Processing", p.475-476.
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)`:
71 >>> from scipy.signal import kaiser_beta
72 >>> kaiser_beta(65)
73 6.20426
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
85def kaiser_atten(numtaps, width):
86 """Compute the attenuation of a Kaiser FIR filter.
88 Given the number of taps `N` and the transition width `width`, compute the
89 attenuation `a` in dB, given by Kaiser's formula:
91 a = 2.285 * (N - 1) * pi * width + 7.95
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.
102 Returns
103 -------
104 a : float
105 The attenuation of the ripple, in dB.
107 See Also
108 --------
109 kaiserord, kaiser_beta
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:
119 >>> from scipy.signal import kaiser_atten
120 >>> kaiser_atten(211, 0.0375)
121 64.48099630593983
123 """
124 a = 2.285 * (numtaps - 1) * np.pi * width + 7.95
125 return a
128def kaiserord(ripple, width):
129 """
130 Determine the filter window parameters for the Kaiser window method.
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`.
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::
146 abs(A(w) - D(w))) < 10**(-ripple/20)
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.
154 Returns
155 -------
156 numtaps : int
157 The length of the Kaiser window.
158 beta : float
159 The beta parameter for the Kaiser window.
161 See Also
162 --------
163 kaiser_beta, kaiser_atten
165 Notes
166 -----
167 There are several ways to obtain the Kaiser window:
169 - ``signal.kaiser(numtaps, beta, sym=True)``
170 - ``signal.get_window(beta, numtaps)``
171 - ``signal.get_window(('kaiser', beta), numtaps)``
173 The empirical equations discovered by Kaiser are used.
175 References
176 ----------
177 Oppenheim, Schafer, "Discrete-Time Signal Processing", pp.475-476.
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.
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%.
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.
192 >>> from scipy.signal import kaiserord, firwin, freqz
193 >>> import matplotlib.pyplot as plt
194 >>> fs = 1000.0
195 >>> cutoff = 175
196 >>> width = 24
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.
203 Use `kaiserord` to determine the length of the filter and the
204 parameter for the Kaiser window.
206 >>> numtaps, beta = kaiserord(65, width/(0.5*fs))
207 >>> numtaps
208 167
209 >>> beta
210 6.20426
212 Use `firwin` to create the FIR filter.
214 >>> taps = firwin(numtaps, cutoff, window=('kaiser', beta),
215 ... scale=False, nyq=0.5*fs)
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.
221 >>> w, h = freqz(taps, worN=8000)
222 >>> w *= 0.5*fs/np.pi # Convert w to Hz.
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.
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
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.
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()
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)
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
258 return int(ceil(numtaps)), beta
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.
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.
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.
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).
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:
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
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.
320 Returns
321 -------
322 h : (numtaps,) ndarray
323 Coefficients of length `numtaps` FIR filter.
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.
333 See Also
334 --------
335 firwin2
336 firls
337 minimum_phase
338 remez
340 Examples
341 --------
342 Low-pass from 0 to f:
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])
350 Use a specific window function:
352 >>> signal.firwin(numtaps, f, window='nuttall')
353 array([ 3.56607041e-04, 9.99286786e-01, 3.56607041e-04])
355 High-pass ('stop' from 0 to f):
357 >>> signal.firwin(numtaps, f, pass_zero=False)
358 array([-0.00859313, 0.98281375, -0.00859313])
360 Band-pass:
362 >>> f1, f2 = 0.1, 0.2
363 >>> signal.firwin(numtaps, [f1, f2], pass_zero=False)
364 array([ 0.06301614, 0.88770441, 0.06301614])
366 Band-stop:
368 >>> signal.firwin(numtaps, [f1, f2])
369 array([-0.00801395, 1.0160279 , -0.00801395])
371 Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1]):
373 >>> f3, f4 = 0.3, 0.4
374 >>> signal.firwin(numtaps, [f1, f2, f3, f4])
375 array([-0.01376344, 1.02752689, -0.01376344])
377 Multi-band (passbands are [f1, f2] and [f3,f4]):
379 >>> signal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
380 array([ 0.04890915, 0.91284326, 0.04890915])
382 """ # noqa: E501
383 # The major enhancements to this function added in November 2010 were
384 # developed by Tom Krauss (see ticket #902).
386 nyq = 0.5 * _get_fs(fs, nyq)
388 cutoff = np.atleast_1d(cutoff) / float(nyq)
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.")
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)
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
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.")
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))
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)
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)
460 # Get and apply the window function.
461 from .signaltools import get_window
462 win = get_window(window, numtaps, fftbins=False)
463 h *= win
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
479 return h
482# Original version of firwin2 from scipy ticket #457, submitted by "tash".
483#
484# Rewritten by Warren Weckesser, 2010.
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.
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.
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.
531 Returns
532 -------
533 taps : ndarray
534 The filter coefficients of the FIR filter, as a 1-D array of length
535 `numtaps`.
537 See also
538 --------
539 firls
540 firwin
541 minimum_phase
542 remez
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.
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:
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
561 Magnitude response of all but type I filters are subjects to following
562 constraints:
564 - type II -- zero at the Nyquist frequency
565 - type III -- zero at zero and Nyquist frequencies
566 - type IV -- zero at zero frequency
568 .. versionadded:: 0.9.0
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.)
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
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:
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]
589 """
590 nyq = 0.5 * _get_fs(fs, nyq)
592 if len(freq) != len(gain):
593 raise ValueError('freq and gain must be of same length.')
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))
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')
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
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.")
634 if nfreqs is None:
635 nfreqs = 1 + 2 ** int(ceil(log(numtaps, 2)))
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))
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)
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
662 fx2 = fx * shift
664 # Use irfft to compute the inverse FFT.
665 out_full = irfft(fx2)
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
674 # Keep only the first `numtaps` coefficients in `out`, and multiply by
675 # the window.
676 out = out_full[:numtaps] * wind
678 if ftype == 3:
679 out[out.size // 2] = 0.0
681 return out
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.
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.
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:
715 * 'bandpass' : flat response in bands. This is the default.
717 * 'differentiator' : frequency proportional response in bands.
719 * 'hilbert' : filter with odd symmetry, that is, type III
720 (for even order) or type IV (for odd order)
721 linear phase filters.
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.
731 Returns
732 -------
733 out : ndarray
734 A rank-1 array containing the coefficients of the optimal
735 (in a minimax sense) filter.
737 See Also
738 --------
739 firls
740 firwin
741 firwin2
742 minimum_phase
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.
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.
762 >>> from scipy import signal
763 >>> import matplotlib.pyplot as plt
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)
777 This example shows a steep low pass transition according to the small
778 transition width and high filter order:
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")
788 This example shows a high pass filter:
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")
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:
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")
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.
818 The next example shows a bandstop filter. Because of the high filter order
819 the transition is quite steep:
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")
830 >>> plt.show()
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
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'")
847 # Convert weight
848 if weight is None:
849 weight = [1] * len(desired)
851 bands = np.asarray(bands).copy()
852 return sigtools._remez(numtaps, bands, desired, weight, tnum, fs,
853 maxiter, grid_density)
856def firls(numtaps, bands, desired, weight=None, nyq=None, fs=None):
857 """
858 FIR filter design using least-squares error minimization.
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).
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.
890 Returns
891 -------
892 coeffs : ndarray
893 Coefficients of the optimal (in a least squares sense) FIR filter.
895 See also
896 --------
897 firwin
898 firwin2
899 minimum_phase
900 remez
902 Notes
903 -----
904 This implementation follows the algorithm given in [1]_.
905 As noted there, least squares design has multiple advantages:
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.
913 This function constructs a Type I linear phase FIR filter, which
914 contains an odd number of `coeffs` satisfying for :math:`n < numtaps`:
916 .. math:: coeffs(n) = coeffs(numtaps - 1 - n)
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).
922 .. versionadded:: 0.18
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
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:
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()
964 """ # noqa
965 nyq = 0.5 * _get_fs(fs, nyq)
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
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)
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.")
1006 # Set up the linear matrix equation to be solved, Qa = b
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.
1011 # We omit the factor of 0.5 above, instead adding it during coefficient
1012 # calculation.
1014 # We also omit the 1/π from both Q and b equations, as they cancel
1015 # during solving.
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)
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
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)
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]
1063 # make coefficients symmetric (linear phase)
1064 coeffs = np.hstack((a[:0:-1], 2 * a[0], a[1:]))
1065 return coeffs
1068def _dhtm(mag):
1069 """Compute the modified 1-D discrete Hilbert transform
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
1090def minimum_phase(h, method='homomorphic', n_fft=None):
1091 """Convert a linear-phase FIR filter to minimum phase
1093 Parameters
1094 ----------
1095 h : array
1096 Linear-phase FIR filter coefficients.
1097 method : {'hilbert', 'homomorphic'}
1098 The method to use:
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.
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.
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).
1115 Returns
1116 -------
1117 h_minimum : array
1118 The minimum-phase version of the filter, with length
1119 ``(length(h) + 1) // 2``.
1121 See Also
1122 --------
1123 firwin
1124 firwin2
1125 remez
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.
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::
1136 epsilon = 2. * n_stop / n_fft
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::
1143 n_fft = 2 ** int(np.ceil(np.log2(2 * (len(h) - 1) / 0.01)))
1145 This gives reasonable results for both the Hilbert and homomorphic
1146 methods, and gives the value used when ``n_fft=None``.
1148 Alternative implementations exist for creating minimum-phase filters,
1149 including zero inversion [2]_ and spectral factorization [3]_ [4]_.
1150 For more information, see:
1152 http://dspguru.com/dsp/howtos/how-to-design-minimum-phase-fir-filters
1154 Examples
1155 --------
1156 Create an optimal linear-phase filter, then convert it to minimum phase:
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.)
1164 Convert it to minimum phase:
1166 >>> h_min_hom = minimum_phase(h_linear, method='homomorphic')
1167 >>> h_min_hil = minimum_phase(h_linear, method='hilbert')
1169 Compare the three filters:
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()
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]