signal.fir_filter_design

Functions for FIR filter design.

Module Contents

Functions

_get_fs(fs,nyq) Utility for replacing the argument ‘nyq’ (with default 1) with ‘fs’.
kaiser_beta(a) Compute the Kaiser parameter beta, given the attenuation a.
kaiser_atten(numtaps,width) Compute the attenuation of a Kaiser FIR filter.
kaiserord(ripple,width) Determine the filter window parameters for the Kaiser window method.
firwin(numtaps,cutoff,width=None,window=”hamming”,pass_zero=True,scale=True,nyq=None,fs=None) FIR filter design using the window method.
firwin2(numtaps,freq,gain,nfreqs=None,window=”hamming”,nyq=None,antisymmetric=False,fs=None) FIR filter design using the window method.
remez(numtaps,bands,desired,weight=None,Hz=None,type=”bandpass”,maxiter=25,grid_density=16,fs=None) Calculate the minimax optimal filter using the Remez exchange algorithm.
firls(numtaps,bands,desired,weight=None,nyq=None,fs=None) FIR filter design using least-squares error minimization.
_dhtm(mag) Compute the modified 1D discrete Hilbert transform
minimum_phase(h,method=”homomorphic”,n_fft=None) Convert a linear-phase FIR filter to minimum phase
_get_fs(fs, nyq)

Utility for replacing the argument ‘nyq’ (with default 1) with ‘fs’.

kaiser_beta(a)

Compute the Kaiser parameter beta, given the attenuation a.

a : float
The desired attenuation in the stopband and maximum ripple in the passband, in dB. This should be a positive number.
beta : float
The beta parameter to be used in the formula for a Kaiser window.

Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.

Suppose we want to design a lowpass filter, with 65 dB attenuation in the stop band. The Kaiser window parameter to be used in the window method is computed by kaiser_beta(65):

>>> from scipy.signal import kaiser_beta
>>> kaiser_beta(65)
6.20426
kaiser_atten(numtaps, width)

Compute the attenuation of a Kaiser FIR filter.

Given the number of taps N and the transition width width, compute the attenuation a in dB, given by Kaiser’s formula:

a = 2.285 * (N - 1) * pi * width + 7.95
numtaps : int
The number of taps in the FIR filter.
width : float
The desired width of the transition region between passband and stopband (or, in general, at any discontinuity) for the filter, expressed as a fraction of the Nyquist frequency.
a : float
The attenuation of the ripple, in dB.

kaiserord, kaiser_beta

Suppose we want to design a FIR filter using the Kaiser window method that will have 211 taps and a transition width of 9 Hz for a signal that is sampled at 480 Hz. Expressed as a fraction of the Nyquist frequency, the width is 9/(0.5*480) = 0.0375. The approximate attenuation (in dB) is computed as follows:

>>> from scipy.signal import kaiser_atten
>>> kaiser_atten(211, 0.0375)
64.48099630593983
kaiserord(ripple, width)

Determine the filter window parameters for the Kaiser window method.

The parameters returned by this function are generally used to create a finite impulse response filter using the window method, with either firwin or firwin2.

ripple : float

Upper bound for the deviation (in dB) of the magnitude of the filter’s frequency response from that of the desired filter (not including frequencies in any transition intervals). That is, if w is the frequency expressed as a fraction of the Nyquist frequency, A(w) is the actual frequency response of the filter and D(w) is the desired frequency response, the design requirement is that:

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

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

width : float
Width of transition region, normalized so that 1 corresponds to pi radians / sample. That is, the frequency is expressed as a fraction of the Nyquist frequency.
numtaps : int
The length of the Kaiser window.
beta : float
The beta parameter for the Kaiser window.

kaiser_beta, kaiser_atten

There are several ways to obtain the Kaiser window:

  • signal.kaiser(numtaps, beta, sym=True)
  • signal.get_window(beta, numtaps)
  • signal.get_window(('kaiser', beta), numtaps)

The empirical equations discovered by Kaiser are used.

Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.

We will use the Kaiser window method to design a lowpass FIR filter for a signal that is sampled at 1000 Hz.

We want at least 65 dB rejection in the stop band, and in the pass band the gain should vary no more than 0.5%.

We want a cutoff frequency of 175 Hz, with a transition between the pass band and the stop band of 24 Hz. That is, in the band [0, 163], the gain varies no more than 0.5%, and in the band [187, 500], the signal is attenuated by at least 65 dB.

>>> from scipy.signal import kaiserord, firwin, freqz
>>> import matplotlib.pyplot as plt
>>> fs = 1000.0
>>> cutoff = 175
>>> width = 24

The Kaiser method accepts just a single parameter to control the pass band ripple and the stop band rejection, so we use the more restrictive of the two. In this case, the pass band ripple is 0.005, or 46.02 dB, so we will use 65 dB as the design parameter.

Use kaiserord to determine the length of the filter and the parameter for the Kaiser window.

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

Use firwin to create the FIR filter.

>>> taps = firwin(numtaps, cutoff, window=('kaiser', beta),
...               scale=False, nyq=0.5*fs)

Compute the frequency response of the filter. w is the array of frequencies, and h is the corresponding complex array of frequency responses.

>>> w, h = freqz(taps, worN=8000)
>>> w *= 0.5*fs/np.pi  # Convert w to Hz.

Compute the deviation of the magnitude of the filter’s response from that of the ideal lowpass filter. Values in the transition region are set to nan, so they won’t appear in the plot.

>>> ideal = w < cutoff  # The "ideal" frequency response.
>>> deviation = np.abs(np.abs(h) - ideal)
>>> deviation[(w > cutoff - 0.5*width) & (w < cutoff + 0.5*width)] = np.nan

Plot the deviation. A close look at the left end of the stop band shows that the requirement for 65 dB attenuation is violated in the first lobe by about 0.125 dB. This is not unusual for the Kaiser window method.

>>> plt.plot(w, 20*np.log10(np.abs(deviation)))
>>> plt.xlim(0, 0.5*fs)
>>> plt.ylim(-90, -60)
>>> plt.grid(alpha=0.25)
>>> plt.axhline(-65, color='r', ls='--', alpha=0.3)
>>> plt.xlabel('Frequency (Hz)')
>>> plt.ylabel('Deviation from ideal (dB)')
>>> plt.title('Lowpass Filter Frequency Response')
>>> plt.show()
firwin(numtaps, cutoff, width=None, window="hamming", pass_zero=True, scale=True, nyq=None, fs=None)

FIR filter design using the window method.

This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if numtaps is odd and Type II if numtaps is even.

Type II filters always have zero response at the Nyquist frequency, so a ValueError exception is raised if firwin is called with numtaps even and having a passband whose right end is at the Nyquist frequency.

numtaps : int
Length of the filter (number of coefficients, i.e. the filter order + 1). numtaps must be even if a passband includes the Nyquist frequency.
cutoff : float or 1D array_like
Cutoff frequency of filter (expressed in the same units as nyq) OR an array of cutoff frequencies (that is, band edges). In the latter case, the frequencies in cutoff should be positive and monotonically increasing between 0 and nyq. The values 0 and nyq must not be included in cutoff.
width : float or None, optional
If width is not None, then assume it is the approximate width of the transition region (expressed in the same units as nyq) for use in Kaiser FIR filter design. In this case, the window argument is ignored.
window : string or tuple of string and parameter values, optional
Desired window to use. See scipy.signal.get_window for a list of windows and required parameters.
pass_zero : bool, optional
If True, the gain at the frequency 0 (i.e. the “DC gain”) is 1. Otherwise the DC gain is 0.
scale : bool, optional

Set to True to scale the coefficients so that the frequency response is exactly unity at a certain frequency. That frequency is either:

  • 0 (DC) if the first passband starts at 0 (i.e. pass_zero is True)
  • nyq (the Nyquist frequency) if the first passband ends at nyq (i.e the filter is a single band highpass filter); center of first passband otherwise
nyq : float, optional
Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in cutoff must be between 0 and nyq. Default is 1.
fs : float, optional
The sampling frequency of the signal. Each frequency in cutoff must be between 0 and fs/2. Default is 2.
h : (numtaps,) ndarray
Coefficients of length numtaps FIR filter.
ValueError
If any value in cutoff is less than or equal to 0 or greater than or equal to fs/2, if the values in cutoff are not strictly monotonically increasing, or if numtaps is even but a passband includes the Nyquist frequency.

firwin2 firls minimum_phase remez

Low-pass from 0 to f:

>>> from scipy import signal
>>> numtaps = 3
>>> f = 0.1
>>> signal.firwin(numtaps, f)
array([ 0.06799017,  0.86401967,  0.06799017])

Use a specific window function:

>>> signal.firwin(numtaps, f, window='nuttall')
array([  3.56607041e-04,   9.99286786e-01,   3.56607041e-04])

High-pass (‘stop’ from 0 to f):

>>> signal.firwin(numtaps, f, pass_zero=False)
array([-0.00859313,  0.98281375, -0.00859313])

Band-pass:

>>> f1, f2 = 0.1, 0.2
>>> signal.firwin(numtaps, [f1, f2], pass_zero=False)
array([ 0.06301614,  0.88770441,  0.06301614])

Band-stop:

>>> signal.firwin(numtaps, [f1, f2])
array([-0.00801395,  1.0160279 , -0.00801395])

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

>>> f3, f4 = 0.3, 0.4
>>> signal.firwin(numtaps, [f1, f2, f3, f4])
array([-0.01376344,  1.02752689, -0.01376344])

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

>>> signal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
array([ 0.04890915,  0.91284326,  0.04890915])
firwin2(numtaps, freq, gain, nfreqs=None, window="hamming", nyq=None, antisymmetric=False, fs=None)

FIR filter design using the window method.

From the given frequencies freq and corresponding gains gain, this function constructs an FIR filter with linear phase and (approximately) the given frequency response.

numtaps : int
The number of taps in the FIR filter. numtaps must be less than nfreqs.
freq : array_like, 1D
The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being Nyquist. The Nyquist frequency is half fs. The values in freq must be nondecreasing. A value can be repeated once to implement a discontinuity. The first value in freq must be 0, and the last value must be fs/2.
gain : array_like
The filter gains at the frequency sampling points. Certain constraints to gain values, depending on the filter type, are applied, see Notes for details.
nfreqs : int, optional
The size of the interpolation mesh used to construct the filter. For most efficient behavior, this should be a power of 2 plus 1 (e.g, 129, 257, etc). The default is one more than the smallest power of 2 that is not less than numtaps. nfreqs must be greater than numtaps.
window : string or (string, float) or float, or None, optional
Window function to use. Default is “hamming”. See scipy.signal.get_window for the complete list of possible values. If None, no window function is applied.
nyq : float, optional
Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in freq must be between 0 and nyq. Default is 1.
antisymmetric : bool, optional
Whether resulting impulse response is symmetric/antisymmetric. See Notes for more details.
fs : float, optional
The sampling frequency of the signal. Each frequency in cutoff must be between 0 and fs/2. Default is 2.
taps : ndarray
The filter coefficients of the FIR filter, as a 1-D array of length numtaps.

firls firwin minimum_phase remez

From the given set of frequencies and gains, the desired response is constructed in the frequency domain. The inverse FFT is applied to the desired response to create the associated convolution kernel, and the first numtaps coefficients of this kernel, scaled by window, are returned.

The FIR filter will have linear phase. The type of filter is determined by the value of ‘numtaps` and antisymmetric flag. There are four possible combinations:

  • odd numtaps, antisymmetric is False, type I filter is produced
  • even numtaps, antisymmetric is False, type II filter is produced
  • odd numtaps, antisymmetric is True, type III filter is produced
  • even numtaps, antisymmetric is True, type IV filter is produced

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

  • type II – zero at the Nyquist frequency
  • type III – zero at zero and Nyquist frequencies
  • type IV – zero at zero frequency

New in version 0.9.0.

[1]Oppenheim, A. V. and Schafer, R. W., “Discrete-Time Signal Processing”, Prentice-Hall, Englewood Cliffs, New Jersey (1989). (See, for example, Section 7.4.)
[2]Smith, Steven W., “The Scientist and Engineer’s Guide to Digital Signal Processing”, Ch. 17. http://www.dspguide.com/ch17/1.htm

A lowpass FIR filter with a response that is 1 on [0.0, 0.5], and that decreases linearly on [0.5, 1.0] from 1 to 0:

>>> from scipy import signal
>>> taps = signal.firwin2(150, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0])
>>> print(taps[72:78])
[-0.02286961 -0.06362756  0.57310236  0.57310236 -0.06362756 -0.02286961]
remez(numtaps, bands, desired, weight=None, Hz=None, type="bandpass", maxiter=25, grid_density=16, fs=None)

Calculate the minimax optimal filter using the Remez exchange algorithm.

Calculate the filter-coefficients for the finite impulse response (FIR) filter whose transfer function minimizes the maximum error between the desired gain and the realized gain in the specified frequency bands using the Remez exchange algorithm.

numtaps : int
The desired number of taps in the filter. The number of taps is the number of terms in the filter, or the filter order plus one.
bands : array_like
A monotonic sequence containing the band edges. All elements must be non-negative and less than half the sampling frequency as given by fs.
desired : array_like
A sequence half the size of bands containing the desired gain in each of the specified bands.
weight : array_like, optional
A relative weighting to give to each band region. The length of weight has to be half the length of bands.
Hz : scalar, optional
Deprecated. Use `fs` instead. The sampling frequency in Hz. Default is 1.
type : {‘bandpass’, ‘differentiator’, ‘hilbert’}, optional

The type of filter:

  • ‘bandpass’ : flat response in bands. This is the default.
  • ‘differentiator’ : frequency proportional response in bands.
  • ‘hilbert’ : filter with odd symmetry, that is, type III
    (for even order) or type IV (for odd order) linear phase filters.
maxiter : int, optional
Maximum number of iterations of the algorithm. Default is 25.
grid_density : int, optional
Grid density. The dense grid used in remez is of size (numtaps + 1) * grid_density. Default is 16.
fs : float, optional
The sampling frequency of the signal. Default is 1.
out : ndarray
A rank-1 array containing the coefficients of the optimal (in a minimax sense) filter.

firls firwin firwin2 minimum_phase

[1]J. H. McClellan and T. W. Parks, “A unified approach to the design of optimum FIR linear phase digital filters”, IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973.
[2]J. H. McClellan, T. W. Parks and L. R. Rabiner, “A Computer Program for Designing Optimum FIR Linear Phase Digital Filters”, IEEE Trans. Audio Electroacoust., vol. AU-21, pp. 506-525, 1973.

For a signal sampled at 100 Hz, we want to construct a filter with a passband at 20-40 Hz, and stop bands at 0-10 Hz and 45-50 Hz. Note that this means that the behavior in the frequency ranges between those bands is unspecified and may overshoot.

>>> from scipy import signal
>>> fs = 100
>>> bpass = signal.remez(72, [0, 10, 20, 40, 45, 50], [0, 1, 0], fs=fs)
>>> freq, response = signal.freqz(bpass)
>>> import matplotlib.pyplot as plt
>>> plt.semilogy(0.5*fs*freq/np.pi, np.abs(response), 'b-')
>>> plt.grid(alpha=0.25)
>>> plt.xlabel('Frequency (Hz)')
>>> plt.ylabel('Gain')
>>> plt.show()
firls(numtaps, bands, desired, weight=None, nyq=None, fs=None)

FIR filter design using least-squares error minimization.

Calculate the filter coefficients for the linear-phase finite impulse response (FIR) filter which has the best approximation to the desired frequency response described by bands and desired in the least squares sense (i.e., the integral of the weighted mean-squared error within the specified bands is minimized).

numtaps : int
The number of taps in the FIR filter. numtaps must be odd.
bands : array_like
A monotonic nondecreasing sequence containing the band edges in Hz. All elements must be non-negative and less than or equal to the Nyquist frequency given by nyq.
desired : array_like
A sequence the same size as bands containing the desired gain at the start and end point of each band.
weight : array_like, optional
A relative weighting to give to each band region when solving the least squares problem. weight has to be half the size of bands.
nyq : float, optional
Deprecated. Use `fs` instead. Nyquist frequency. Each frequency in bands must be between 0 and nyq (inclusive). Default is 1.
fs : float, optional
The sampling frequency of the signal. Each frequency in bands must be between 0 and fs/2 (inclusive). Default is 2.
coeffs : ndarray
Coefficients of the optimal (in a least squares sense) FIR filter.

firwin firwin2 minimum_phase remez

This implementation follows the algorithm given in [1]_. As noted there, least squares design has multiple advantages:

  1. Optimal in a least-squares sense.
  2. Simple, non-iterative method.
  3. The general solution can obtained by solving a linear system of equations.
  4. Allows the use of a frequency dependent weighting function.

This function constructs a Type I linear phase FIR filter, which contains an odd number of coeffs satisfying for :

The odd number of coefficients and filter symmetry avoid boundary conditions that could otherwise occur at the Nyquist and 0 frequencies (e.g., for Type II, III, or IV variants).

New in version 0.18.

[1]Ivan Selesnick, Linear-Phase Fir Filter Design By Least Squares. OpenStax CNX. Aug 9, 2005. http://cnx.org/contents/eb1ecb35-03a9-4610-ba87-41cd771c95f2@7

We want to construct a band-pass filter. Note that the behavior in the frequency ranges between our stop bands and pass bands is unspecified, and thus may overshoot depending on the parameters of our filter:

>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> fig, axs = plt.subplots(2)
>>> fs = 10.0  # Hz
>>> desired = (0, 0, 1, 1, 0, 0)
>>> for bi, bands in enumerate(((0, 1, 2, 3, 4, 5), (0, 1, 2, 4, 4.5, 5))):
...     fir_firls = signal.firls(73, bands, desired, fs=fs)
...     fir_remez = signal.remez(73, bands, desired[::2], fs=fs)
...     fir_firwin2 = signal.firwin2(73, bands, desired, fs=fs)
...     hs = list()
...     ax = axs[bi]
...     for fir in (fir_firls, fir_remez, fir_firwin2):
...         freq, response = signal.freqz(fir)
...         hs.append(ax.semilogy(0.5*fs*freq/np.pi, np.abs(response))[0])
...     for band, gains in zip(zip(bands[::2], bands[1::2]),
...                            zip(desired[::2], desired[1::2])):
...         ax.semilogy(band, np.maximum(gains, 1e-7), 'k--', linewidth=2)
...     if bi == 0:
...         ax.legend(hs, ('firls', 'remez', 'firwin2'),
...                   loc='lower center', frameon=False)
...     else:
...         ax.set_xlabel('Frequency (Hz)')
...     ax.grid(True)
...     ax.set(title='Band-pass %d-%d Hz' % bands[2:4], ylabel='Magnitude')
...
>>> fig.tight_layout()
>>> plt.show()
_dhtm(mag)

Compute the modified 1D discrete Hilbert transform

mag : ndarray
The magnitude spectrum. Should be 1D with an even length, and preferably a fast length for FFT/IFFT.
minimum_phase(h, method="homomorphic", n_fft=None)

Convert a linear-phase FIR filter to minimum phase

h : array
Linear-phase FIR filter coefficients.
method : {‘hilbert’, ‘homomorphic’}

The method to use:

‘homomorphic’ (default)
This method [4] [5] works best with filters with an odd number of taps, and the resulting minimum phase filter will have a magnitude response that approximates the square root of the the original filter’s magnitude response.
‘hilbert’
This method [1]_ is designed to be used with equiripple filters (e.g., from remez) with unity or zero gain regions.
n_fft : int
The number of points to use for the FFT. Should be at least a few times larger than the signal length (see Notes).
h_minimum : array
The minimum-phase version of the filter, with length (length(h) + 1) // 2.

firwin firwin2 remez

Both the Hilbert [1]_ or homomorphic [4] [5] methods require selection of an FFT length to estimate the complex cepstrum of the filter.

In the case of the Hilbert method, the deviation from the ideal spectrum epsilon is related to the number of stopband zeros n_stop and FFT length n_fft as:

epsilon = 2. * n_stop / n_fft

For example, with 100 stopband zeros and a FFT length of 2048, epsilon = 0.0976. If we conservatively assume that the number of stopband zeros is one less than the filter length, we can take the FFT length to be the next power of 2 that satisfies epsilon=0.01 as:

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

This gives reasonable results for both the Hilbert and homomorphic methods, and gives the value used when n_fft=None.

Alternative implementations exist for creating minimum-phase filters, including zero inversion [2]_ and spectral factorization [3] [4]. For more information, see:

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

>>> from scipy.signal import remez, minimum_phase, freqz, group_delay
>>> import matplotlib.pyplot as plt
>>> freq = [0, 0.2, 0.3, 1.0]
>>> desired = [1, 0]
>>> h_linear = remez(151, freq, desired, Hz=2.)

Convert it to minimum phase:

>>> h_min_hom = minimum_phase(h_linear, method='homomorphic')
>>> h_min_hil = minimum_phase(h_linear, method='hilbert')

Compare the three filters:

>>> fig, axs = plt.subplots(4, figsize=(4, 8))
>>> for h, style, color in zip((h_linear, h_min_hom, h_min_hil),
...                            ('-', '-', '--'), ('k', 'r', 'c')):
...     w, H = freqz(h)
...     w, gd = group_delay((h, 1))
...     w /= np.pi
...     axs[0].plot(h, color=color, linestyle=style)
...     axs[1].plot(w, np.abs(H), color=color, linestyle=style)
...     axs[2].plot(w, 20 * np.log10(np.abs(H)), color=color, linestyle=style)
...     axs[3].plot(w, gd, color=color, linestyle=style)
>>> for ax in axs:
...     ax.grid(True, color='0.5')
...     ax.fill_between(freq[1:3], *ax.get_ylim(), color='#ffeeaa', zorder=1)
>>> axs[0].set(xlim=[0, len(h_linear) - 1], ylabel='Amplitude', xlabel='Samples')
>>> axs[1].legend(['Linear', 'Min-Hom', 'Min-Hil'], title='Phase')
>>> for ax, ylim in zip(axs[1:], ([0, 1.1], [-150, 10], [-60, 60])):
...     ax.set(xlim=[0, 1], ylim=ylim, xlabel='Frequency')
>>> axs[1].set(ylabel='Magnitude')
>>> axs[2].set(ylabel='Magnitude (dB)')
>>> axs[3].set(ylabel='Group delay')
>>> plt.tight_layout()
[1]N. Damera-Venkata and B. L. Evans, “Optimal design of real and complex minimum phase digital FIR filters,” Acoustics, Speech, and Signal Processing, 1999. Proceedings., 1999 IEEE International Conference on, Phoenix, AZ, 1999, pp. 1145-1148 vol.3. doi: 10.1109/ICASSP.1999.756179
[2]X. Chen and T. W. Parks, “Design of optimal minimum phase FIR filters by direct factorization,” Signal Processing, vol. 10, no. 4, pp. 369–383, Jun. 1986.
[3]T. Saramaki, “Finite Impulse Response Filter Design,” in Handbook for Digital Signal Processing, chapter 4, New York: Wiley-Interscience, 1993.
[4](1, 2, 3) J. S. Lim, Advanced Topics in Signal Processing. Englewood Cliffs, N.J.: Prentice Hall, 1988.
[5](1, 2) A. V. Oppenheim, R. W. Schafer, and J. R. Buck, “Discrete-Time Signal Processing,” 2nd edition. Upper Saddle River, N.J.: Prentice Hall, 1999.