Remez

openseize.filtering.fir.Remez

Bases: FIR

A Parks-McClellan optimal Chebyshev FIR filter.

This FIR is designed by minimizing:

E(f) = W(f) |A(f) - D(f)|

where E(f) is a the weighted difference of the actual frequency response A(F) from the desired frequency response D(f) for f in [0, nyquist] and W(f) is a set of weights.

Attributes:

Name	Type	Description
bands	Sequence	A monitonically increasing sequence of pass and stop band edge frequencies that must include 0 and the nyquist frequencies.
desired	Sequence	A sequence of desired gains in each band of bands.
fs	int	The sampling rate of the digital system.
gpass	float	The maximum ripple in the pass band(s) (dB).
gstop	float	float The minimum attenuation in the stop band(s) (dB).
kwargs	float	Any valid kwarg for scipy's signal.remez function. numtaps (int): The number of taps used to construct this filter. This overrides Remez's internal numtaps property. weight (Sequence): A sequence the same length as desired with weights to relatively weight each band. If no value is passed the weights are the inverse of the percentage loss and attenuation in the pass and stop bands respectively. maxiter (int): The number of iterations to test for convergence. grid_density (int): Resolution of the grid remez uses to minimize E(f). If algorithm converges but has unexpected ripple, increasing this value can help. The default is 16.

Examples:

>>> # Design a bandpass filter that passes 400 to 800 Hz
>>> bands = [0, 300, 400, 800, 900, 2500]
>>> desired = [0, 1, 0]
>>> filt = Remez(bands, desired, fs=5000, gpass=.5, gstop=40)
>>> filt.plot()

Notes

The Remez algorithm may fail to converge. It is recommended that any filter designed with Remez be visually inspected before use.

References

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. Remez algorithm: https://en.wikipedia.org/wiki/Remez_algorithm

Source code in openseize/filtering/fir.py

class Remez(FIR):
    """A Parks-McClellan optimal Chebyshev FIR filter.

    This FIR is designed by minimizing:

    ```math
    E(f) = W(f) |A(f) - D(f)|
    ```

    where E(f) is a the weighted difference of the actual frequency
    response A(F) from the desired frequency response D(f) for f in
    [0, nyquist] and W(f) is a set of weights.

    Attributes:
        bands (Sequence):
            A monitonically increasing sequence of pass and stop band edge
            frequencies that must include 0 and the nyquist frequencies.
        desired (Sequence):
            A sequence of desired gains in each band of bands.
        fs (int):
            The sampling rate of the digital system.
        gpass (float):
           The maximum ripple in the pass band(s) (dB).
        gstop: float
            The minimum attenuation in the stop band(s) (dB).
        kwargs:
            Any valid kwarg for scipy's signal.remez function.

            - numtaps (int):
                The number of taps used to construct this filter. This
                overrides Remez's internal numtaps property.
            - weight (Sequence):
                A sequence the same length as desired with weights to
                relatively weight each band. If no value is passed
                the weights are the inverse of the percentage loss and
                attenuation in the pass and stop bands respectively.
            - maxiter (int):
                The number of iterations to test for convergence.
            - grid_density (int):
                Resolution of the grid remez uses to minimize E(f). If
                algorithm converges but has unexpected ripple, increasing
                this value can help. The default is 16.

    Examples:
        >>> # Design a bandpass filter that passes 400 to 800 Hz
        >>> bands = [0, 300, 400, 800, 900, 2500]
        >>> desired = [0, 1, 0]
        >>> filt = Remez(bands, desired, fs=5000, gpass=.5, gstop=40)
        >>> filt.plot()

    Notes:
        The Remez algorithm may fail to converge. It is recommended that any
        filter designed with Remez be visually inspected before use.

    References:
        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. Remez algorithm: https://en.wikipedia.org/wiki/Remez_algorithm

    """

    def __init__(self,
                 bands: Sequence[float],
                 desired: Sequence[float],
                 fs: int,
                 gpass: float = 1,
                 gstop: float = 40,
                 **kwargs):
        """Initialize this Remez FIR.

        Args:
            bands:
                A monotonically increasing sequence of pass & stop band edge
                frequencies than includes the 0 and nyquist frequencies.
            desired:
                A sequence containing the desired gains (1 or 0) for each
                band in bands.
            fs:
                The sampling rate of the digital system
            gpass:
                The maximum ripple in the pass band (dB). Default is 1 dB
                which is an amplitude loss of ~ 11%. If more than 1 pass
                band is supplied in bands, the same maximum loss will be
                applied to all bands.
            gstop:
                The minimum attenuation in the stop band (dB). The default
                is 40 dB which is an amplitude loss of 99%. If more than
                1 stop band is supplied in bands, the same minimum
                attenuation is applied to all stop bands.
            kwargs:
                Any valid kwarg for scipy's signal.remez function.

                - numtaps (int):
                    The number of taps used to construct this filter. This
                    overrides Remez's internal numtaps property.
                - weight (Sequence):
                    A sequence the same length as desired with weights to
                    relatively weight each band. If no value is passed
                    the weights are the inverse of the percentage loss and
                    attenuation in the pass and stop bands respectively.
                - maxiter (int):
                    The number of iterations to test for convergence.
                - grid_density (int):
                    Resolution of the grid remez uses to minimize E(f). If
                    algorithm converges but has unexpected ripple,
                    increasing this value can help. The default is 16.
        """

        self.bands = np.array(bands).reshape(-1, 2)
        self.desired = np.array(desired, dtype=bool)

        # construct fpass and fstop from bands for plotting
        fp = self.bands[self.desired].flatten()
        fpass = fp[np.logical_and(fp > 0, fp < fs / 2)]
        fst = self.bands[~self.desired].flatten()
        fstop = fst[np.logical_and(fst > 0, fst < fs / 2)]

        # transform gpass and gstop to amplitudes
        self.delta_pass = 1 - 10 ** (-gpass / 20)
        self.delta_stop = 10 ** (-gstop / 20)
        self.delta = (self.delta_pass * self.desired +
                      self.delta_stop * (1 - self.desired))

        super().__init__(fpass, fstop, gpass, gstop, fs, **kwargs)

    @property
    def btype(self):
        """Return the string band type of this filter."""

        fp, fs = self.fpass, self.fstop
        if len(fp) < 2:
            btype = "lowpass" if fp < fs else "highpass"
        elif len(fp) == 2:
            btype = "bandstop" if fp[0] < fs[0] else "bandpass"
        else:
            btype = "multiband"
        return btype

    @property
    def numtaps(self):
        """Estimates the number of taps needed to meet this filters pass and
        stop band specifications.

        This is the Bellanger estimate for the number of taps. Strictly it
        does not apply to multiband filters as it applies a single pass and
        a single stop attenuation to each band. As such the frequency
        response of the filter should be checked to ensure that the pass and
        stop band criteria are being met.

        References:
            M. Bellanger, Digital Processing of Signals: Theory and
            Practice (3rd Edition), Wiley, Hoboken, NJ, 2000.
        """

        dp, ds = self.delta_pass, self.delta_stop
        n = -2 / 3 * np.log10(10 * dp * ds) * self.fs / self.width
        ntaps = int(np.ceil(n))
        return ntaps + 1 if ntaps % 2 == 0 else ntaps

    def _build(self, **kwargs):
        """Returns the coefficients of this Remez."""

        # Get kwargs or use defaults to pass to scipy remez
        ntaps = kwargs.pop("numtaps", self.numtaps)
        weight = kwargs.pop("weight", 1 / self.delta)
        maxiter = kwargs.pop("maxiter", 25)
        grid_density = kwargs.pop("grid_density", 16)

        return sps.remez(ntaps,
                         self.bands.flatten(),
                         self.desired,
                         weight=weight,
                         maxiter=maxiter,
                         grid_density=grid_density,
                         fs=self.fs)

__init__(bands, desired, fs, gpass=1, gstop=40, **kwargs)

Initialize this Remez FIR.

Parameters:

Name	Type	Description	Default
bands	Sequence[float]	A monotonically increasing sequence of pass & stop band edge frequencies than includes the 0 and nyquist frequencies.	required
desired	Sequence[float]	A sequence containing the desired gains (1 or 0) for each band in bands.	required
fs	int	The sampling rate of the digital system	required
gpass	float	The maximum ripple in the pass band (dB). Default is 1 dB which is an amplitude loss of ~ 11%. If more than 1 pass band is supplied in bands, the same maximum loss will be applied to all bands.	1
gstop	float	The minimum attenuation in the stop band (dB). The default is 40 dB which is an amplitude loss of 99%. If more than 1 stop band is supplied in bands, the same minimum attenuation is applied to all stop bands.	40
kwargs		Any valid kwarg for scipy's signal.remez function. numtaps (int): The number of taps used to construct this filter. This overrides Remez's internal numtaps property. weight (Sequence): A sequence the same length as desired with weights to relatively weight each band. If no value is passed the weights are the inverse of the percentage loss and attenuation in the pass and stop bands respectively. maxiter (int): The number of iterations to test for convergence. grid_density (int): Resolution of the grid remez uses to minimize E(f). If algorithm converges but has unexpected ripple, increasing this value can help. The default is 16.	{}

Source code in openseize/filtering/fir.py

def __init__(self,
             bands: Sequence[float],
             desired: Sequence[float],
             fs: int,
             gpass: float = 1,
             gstop: float = 40,
             **kwargs):
    """Initialize this Remez FIR.

    Args:
        bands:
            A monotonically increasing sequence of pass & stop band edge
            frequencies than includes the 0 and nyquist frequencies.
        desired:
            A sequence containing the desired gains (1 or 0) for each
            band in bands.
        fs:
            The sampling rate of the digital system
        gpass:
            The maximum ripple in the pass band (dB). Default is 1 dB
            which is an amplitude loss of ~ 11%. If more than 1 pass
            band is supplied in bands, the same maximum loss will be
            applied to all bands.
        gstop:
            The minimum attenuation in the stop band (dB). The default
            is 40 dB which is an amplitude loss of 99%. If more than
            1 stop band is supplied in bands, the same minimum
            attenuation is applied to all stop bands.
        kwargs:
            Any valid kwarg for scipy's signal.remez function.

            - numtaps (int):
                The number of taps used to construct this filter. This
                overrides Remez's internal numtaps property.
            - weight (Sequence):
                A sequence the same length as desired with weights to
                relatively weight each band. If no value is passed
                the weights are the inverse of the percentage loss and
                attenuation in the pass and stop bands respectively.
            - maxiter (int):
                The number of iterations to test for convergence.
            - grid_density (int):
                Resolution of the grid remez uses to minimize E(f). If
                algorithm converges but has unexpected ripple,
                increasing this value can help. The default is 16.
    """

    self.bands = np.array(bands).reshape(-1, 2)
    self.desired = np.array(desired, dtype=bool)

    # construct fpass and fstop from bands for plotting
    fp = self.bands[self.desired].flatten()
    fpass = fp[np.logical_and(fp > 0, fp < fs / 2)]
    fst = self.bands[~self.desired].flatten()
    fstop = fst[np.logical_and(fst > 0, fst < fs / 2)]

    # transform gpass and gstop to amplitudes
    self.delta_pass = 1 - 10 ** (-gpass / 20)
    self.delta_stop = 10 ** (-gstop / 20)
    self.delta = (self.delta_pass * self.desired +
                  self.delta_stop * (1 - self.desired))

    super().__init__(fpass, fstop, gpass, gstop, fs, **kwargs)

btype() property

Return the string band type of this filter.

Source code in openseize/filtering/fir.py

@property
def btype(self):
    """Return the string band type of this filter."""

    fp, fs = self.fpass, self.fstop
    if len(fp) < 2:
        btype = "lowpass" if fp < fs else "highpass"
    elif len(fp) == 2:
        btype = "bandstop" if fp[0] < fs[0] else "bandpass"
    else:
        btype = "multiband"
    return btype

numtaps() property

Estimates the number of taps needed to meet this filters pass and stop band specifications.

This is the Bellanger estimate for the number of taps. Strictly it does not apply to multiband filters as it applies a single pass and a single stop attenuation to each band. As such the frequency response of the filter should be checked to ensure that the pass and stop band criteria are being met.

References

M. Bellanger, Digital Processing of Signals: Theory and Practice (3rd Edition), Wiley, Hoboken, NJ, 2000.

Source code in openseize/filtering/fir.py

@property
def numtaps(self):
    """Estimates the number of taps needed to meet this filters pass and
    stop band specifications.

    This is the Bellanger estimate for the number of taps. Strictly it
    does not apply to multiband filters as it applies a single pass and
    a single stop attenuation to each band. As such the frequency
    response of the filter should be checked to ensure that the pass and
    stop band criteria are being met.

    References:
        M. Bellanger, Digital Processing of Signals: Theory and
        Practice (3rd Edition), Wiley, Hoboken, NJ, 2000.
    """

    dp, ds = self.delta_pass, self.delta_stop
    n = -2 / 3 * np.log10(10 * dp * ds) * self.fs / self.width
    ntaps = int(np.ceil(n))
    return ntaps + 1 if ntaps % 2 == 0 else ntaps

Bases and Mixins

FIR Base

Bases: abc.ABC, mixins.ViewInstance, FIRViewer

Base class for finite impulse response filters.

Attributes:

Name	Type	Description
fpass	np.ndarray	1-D numpy array of start and stop edge frequencies of this filter's passband(s).
fstop	np.ndarray	1-D numpy array of start and stop edge frequencies of this filter's stopband(s).
gpass	float	Maximum ripple in the passband(s) in dB.
gstop	float	Minimum attenuation in the stopbands in dB.
fs	int	The sampling rate of the digital system.
nyq	float	The nyquist rate of the digital system, fs/2.
width	float	The minimum transition width between the pass and stopbands.
coeffs	np.ndarray	A 1-D numpy array of filter coeffecients.

Notes

This FIR ABC defines the common and expected methods of all concrete FIR filters in the openseize.filtering.fir module. Inheritors must override abstract methods & properties of this base to be instantiable.

ftype() property

Returns the string name of this FIR filter.

pass_attenuation() property

Converts the max passband ripple, gpass, into a pass band attenuation in dB.

cutoff() property

Returns an ndarray of the -6 dB points of each transition band.

__call__(data, chunksize, axis=-1, mode='same', **kwargs)

Apply this filter to an ndarray or producer of ndarrays.

Parameters:

Name	Type	Description	Default
data	Union[Producer, np.ndarray]	The data to be filtered.	required
chunksize	int	The number of samples to hold in memory during filtering.	required
axis	int	The axis of data along which to apply the filter. If data is multidimensional, the filter will be independently applied along all slices of axis.	-1
mode	str	A numpy convolve mode; one of 'full', 'same', 'valid'. Full: This mode includes all points of the convolution of the filter window and data. This mode does not compensate for the delay introduced by the filter. Same: This mode (Default) returns data of the same size as the input. This mode adjust for the delay introduced by the filter. Valid: This mode returns values only when the filter and data completely overlap. The result using this mode is to shift the data (num_taps - 1) / 2 samples to the left of the input data.	'same'
kwargs		Any valid keyword argument for the producer constructor.	{}

Returns:

Type	Description
Union[Producer, np.ndarray]	Filtered result with type matching input 'data' parameter.

Viewer Mixin

A collection of common plotting methods for both IIR, FIR and Parks-McClellan filters.

All filters in openseize have the ability to plot their impulse response and frequency response to a matplotlib figure called the Viewer. This mixin is inherited by specific IIR, FIR and ParksMcClellan Viewers in this file. Each of these specific viewers is inherited by the corresponding filter type (i.e. IIR, FIR, ParksMcClellan) in the openseize filtering module.

plot(size=(8, 6), gridalpha=0.3, worN=2048, rope=-100, axarr=None, show=True)

Plots the impulse and frequency response of this filter.

Parameters:

Name	Type	Description	Default
size	Tuple[int, int]	tuple The figure size to display for the plots. Default is 8 x 6.	(8, 6)
gridalpha	float	float in [0, 1] The alpha transparency of each subplots grid. Default is 0.3	0.3
worN	int	int The number of frequencies to compute the gain and phase responses over. Default is 2048 frequencies.	2048
rope	float	float For plotting, all values below this region of practical equivalence will be set to this value. Default is -100 dB. Any filter response smaller than this will be set to -100 for plotting.	-100
axarr	Optional[Sequence[plt.Axes]]	A Matplotlib axis array. An optional axis array to plot the impulse and frequency responses to. Default None means a new axis is created.	None