Infinite Impulse Reponse (IIR) Filter Design ============================================ .. author jgaeddert Joseph D. Gaeddert .. date January 1, 1975 .. location Boston, MA .. header /doc/iirdes/banner.png .. keywords infinite impulse response, IIR, filter design, iirdes, bilinear z-transform |liquid| implements infinite impulse response (IIR) filter design for the five major classes of filters (Butterworth, Chebyshev type-I, Chebyshev type-II, elliptic, and Bessel) by first computing their analog low-pass prototypes, performing a bilinear :math:`z`-transform to convert to the digital domain, then transforming to the appropriate band type (e.g. high pass) if necessary. Externally, the user may abstract the entire process by using the :api:`liquid_iirdes()` method. Furthermore, if the end result is to create a filter *object* as opposed to computing the coefficients themselves, the :api:`iirfilt_crcf_create_prototype()` method can be used to generate the object directly (see [section-filter-iirfilt]). .. _section-filter-iirdes-iirdes: liquid_iirdes(), the simplified method ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The :api:`liquid_iirdes()` method designs an IIR filter's coefficients from one of the four major types (Butterworth, Chebyshev, elliptic/Cauer, and Bessel) with as minimal an interface as possible. The user specifies the filter prototype, order, cutoff frequency, and other parameters as well as the resulting filter structure (regular or second-order sections), and the function returns the appropriate filter coefficients that meet that design. Specifically, the interface is .. code-block:: c liquid_iirdes(_ftype, _btype, _format, _n, _fc, _f0, _Ap, _As, *_B, *_A); * `_ftype` is the analog filter prototype, e.g. :api:`LIQUID_IIRDES_BUTTER` * `_btype` is the band type, e.g. :api:`LIQUID_IIRDES_BANDPASS` * `_format` is the output format of the coefficients, e.g. :api:`LIQUID_IIRDES_SOS` * `_n` is the filter order * `_fc` is the normalized cutoff frequency of the analog prototype * `_f0` is the normalized center frequency of the analog prototype (only applicable to bandpass and band-stop filter designs, ignored for low-pass and high-pass filter designs) * `_Ap` is the pass-band ripple (only applicable to Chebyshev Type-I and elliptic filter designs, ignored for Butterworth, Chebyshev Type-II, and Bessel designs) * `_As` is the stop-band ripple (only applicable to Chebyshev Type-II and elliptic filter designs, ignored for Butterworth, Chebyshev Type-I, and Bessel designs) * `_B`, `_A` are the output feed-forward (numerator) and feed-back (denominator) coefficients, respectively. The format and size of these arrays depends on the value of the `_format` and `_btype` parameters. To compute the specific lengths of the arrays, first define the effective filter order :math:`N` which is the same as the specified filter order for low- and high- pass filters, and doubled for band-pass and band-stop filters. If the the format is :api:`LIQUID_IIRDES_TF` (the regular transfer function format) then the size of `_B` and `_A` is simply :math:`N`. If, on the other hand, the format is :api:`LIQUID_IIRDES_SOS` (second-order sections format) then a few extra steps are needed: define :math:`r` as :math:`0` when :math:`N` is even and :math:`1` when :math:`N` is odd, and define :math:`L` as :math:`(N-r)/2`. The sizes of `_B` and `_A` for the second-order sections case are each :math:`3(L+r)`. As an example, the following example designs a :math:`5^{th}`-order elliptic band-pass filter with 1 dB ripple in the passband, 60 dB ripple in the stop-band, a cut-off frequency of :math:`f_c/F_s=0.2` and a center frequency :math:`f_0/F_s=0.25`; the frequency response of the resulting filter can be found in [fig-filter-iirdes_example]. .. listing doc/iirdes/iirdes.example.c .. qplot:: iirdes iirdes/iirdes.c :args: -t ellip -n 5 -b BP -r 1 -s 60 -f 0.2 -c 0.25 :width: 100% :name: fig-iirdes-example :caption: Example of the `iirdes` interface, designing a :math:`5^{th}`-order elliptic band-pass filter with 1 dB ripple in the passband, 60 dB ripple in the stop-band, a cut-off frequency of :math:`f_c/F_s=0.2` and a center frequency :math:`f_0/F_s=0.25` .. qplot:: series iirdes/iirdes_comparison.c :kwargs: { "layout": [2,1], "figsize": [12,12], "format": { "xlabel" : "Normalized Frequency [f/Fs]", "xrange" : [0, 0.5], "legend" : true }, "plots": [ { "series": [["f","H-butter"], ["f","H-cheby1"], ["f","H-cheby2"], ["f","H-ellip"], ["f","H-bessel"]], "yrange" : [0, 1], "ylabel" : "Gain" }, { "series": [["f","gd-butter"], ["f","gd-cheby1"], ["f","gd-cheby2"], ["f","gd-ellip"], ["f","gd-bessel"]], "yrange" : [0, 25], "ylabel" : "Group Delay [samples]" } ]} :width: 75% :name: fig-iirdes-comparison :caption: Comparison of the available IIR filter types in |liquid| .. _section-filter-iirdes-internal: internal description ~~~~~~~~~~~~~~~~~~~~ While the user only needs to specify the filter parameters, the internal procedure for computing the coefficients is somewhat complicated. Listed below is the step-by-step process for |liquid|'s IIR filter design procedure. * Use :api:`butterf()`, :api:`cheby1f()`, :api:`cheby2f()`, :api:`ellipf()`, :api:`besself()` to design a low-pass analog prototype :math:`H_a(s)` in terms of complex zeros, poles, and gain. The `azpkf` extension stands for "analog zeros, poles, gain (floating-point)." :api:`butter_azpkf()` Butterworth (maximally flat in the pass-band) :api:`cheby1_azpkf()` Chebyshev Type-I (equiripple in the pass-band) :api:`cheby2_azpkf()` Chebyshev Type-II (equiripple in the stop-band) :api:`ellip_azpkf()` elliptic filter (equiripple in the pass- and stop-bands) :api:`bessel_azpkf()` Bessel (maximally flat group delay) * Compute frequency pre-warping factor, :math:`m`, to set cutoff frequency (and center frequency if designing a band-pass or band-stop filter) using the :api:`iirdes_freqprewarp()` method. * Convert the low-pass analog prototype :math:`H_a(s)` to its digital equivalent :math:`H_d(z)` (also in terms of zeros, poles, and gain) using the bilinear :math:`z`-transform using the :api:`bilinear_zpkf()` method. This maps the analog zeros/poles/gain into digital zeros/poles/gain. * Transform the low-pass digital prototype to high-pass, band-pass, or band-stop using the :api:`iirdes_dzpk_lp2bp()` method. For the band-pass and band-stop cases, the number of poles and zeros will need to be doubled. * LP low-pass filter : :math:`s = m (1+z^{-1}) / (1-z^{-1})` * HP high-pass filter : :math:`s = m (1-z^{-1}) / (1+z^{-1})` * BP band-pass filter : :math:`s = m (1-c_0 z^{-1}+z^{-2}) / (1-z^{-2})` * BS band-stop filter : :math:`s = m (1-z^{-2}) / (1-c_0 z^{-1}+z^{-2})` * Transform the digital :math:`z/p/k` form of the filter to one of the two forms: * TF typical transfer function for digital iir filters of the form :math:`B(z)/A(z)`, :api:`iirdes_dzpk2tff()` * SOS second-order sections form : :math:`\prod_k{ B_k(z)/A_k(z) }`, :api:`iirdes_dzpk2sosf()`. This is the preferred method. A simplified example for this procedure is given in `examples/iirdes_example.c`. Bilinear z-transform ~~~~~~~~~~~~~~~~~~~~ The bilinear :math:`z`-transform converts an analog prototype to its digital counterpart. Given a continuous time analog transfer function in zeros/poles/gain form ("zpk") with :math:`n_z` zeros and :math:`n_p` poles .. math:: H_a(s) = k_a \frac{ (s-z_{a0})(s-z_{a1})\cdots(s-z_{an_z-1}) }{ (s-p_{a0})(s-p_{a1})\cdots(s-p_{an_p-1}) } :label: eqn-filter-iirdes-Ha the bilinear :math:`z`-transform converts :math:`H_a(s)` into the discrete transfer function :math:`H_d(z)` by mapping the :math:`s`-plane onto the :math:`z`-plane with the approximation .. math:: s \approx \frac{2}{T} \frac{1-z^{-1}}{1 + z^{-1}} :label: eqn-filter-iirdes-bilinear This maps :math:`H_a(0) \rightarrow H_d(0)` and :math:`H_a(\infty) \rightarrow H_d(\omega_s/2)`, however we are free to choose the pre-warping factor which maps the cutoff frequency :math:`\omega_c`. .. math:: s \rightarrow \omega_c \cot\left(\frac{\pi \omega_c}{\omega_s}\right) \frac{1-z^{-1}}{1+z^{-1}} :label: eqn-filter-iirdes-bilinear_prewarp Substituting this into :math:`H_a(s)` gives the discrete-time transfer function .. math:: H(z) = k_a \frac{ \left(m\frac{1-z^{-1}}{1+z^{-1}}-z_{a0}\right) \left(m\frac{1-z^{-1}}{1+z^{-1}}-z_{a1}\right) \cdots \left(m\frac{1-z^{-1}}{1+z^{-1}}-z_{an_z-1}\right) }{ \left(m\frac{1-z^{-1}}{1+z^{-1}}-p_{a0}\right) \left(m\frac{1-z^{-1}}{1+z^{-1}}-p_{a1}\right) \cdots \left(m\frac{1-z^{-1}}{1+z^{-1}}-p_{an_p-1}\right) } :label: eqn-filter-iirdes-H where :math:`m=\omega_c \cot\left(\pi \omega_c / \omega_s\right)` is the frequency pre-warping factor, computed in |liquid| via the method :api:`iirdes_freqprewarp()`. Multiplying both the numerator an denominator by :math:`(1+z^{-1})^{n_p}` and applying some algebraic manipulation results in the digital filter .. math:: H_d(s) = k_d \frac{ (1-z_{d0}z^{-1})(1-z_{d1}z^{-1})\cdots(1-z_{dn-1}z^{-1}) }{ (1-p_{d0}z^{-1})(1-p_{d1}z^{-1})\cdots(1-p_{dn-1}z^{-1}) } :label: eqn-filter-iirdes-Hd The :api:`bilinear_zpk()` method in |liquid| transforms the the analog zeros (:math:`z_{ak}`), poles (:math:`p_{ak}`), and gain (:math:`H_0`) into their digital equivalents (:math:`z_{dk}`, :math:`p_{dk}`, :math:`G_0`). For a filter with :math:`n_z` analog zeros :math:`z_{ak}` the digital zeros :math:`z_{dk}` are computed as .. math:: z_{dk} = \begin{cases} \frac{1 + m z_{ak}}{1 - m z_{ak}} & k \lt n_z \\ -1 & \text{otherwise} \end{cases} :label: eqn-filter-iirdes-bilinear_zeros where :math:`m` is the pre-warping factor. For a filter with :math:`n_p` analog poles :math:`p_{ak}` the digital poles :math:`p_{dk}` are computed as .. math:: p_{dk} = \frac{1 + m p_{ak}}{1 - m p_{ak}} :label: eqn-filter-iirdes-bilinear_poles Keeping in mind that an analog filter's order is defined by its number of poles, the digital gain can be computed as .. math:: G_0 = H_0 \prod_{k=0}^{n_p-1}{ \frac{1 - p_{dk}}{1 - z_{dk}} } :label: eqn-filter-iirdes-bilinear_gain .. _section-filter-iirdes-transformations: Filter transformations ~~~~~~~~~~~~~~~~~~~~~~ The prototype low-pass digital filter can be converted into a high-pass, band-pass, or band-stop filter using a combination of the following filter transformations in |liquid|: :api:`iirdes_dzpk_lp2hp(*_zd,*_pd,_n,*_zdt,*_pdt)` Converts a low-pass digital prototype :math:`H_d(z)` to a high-pass prototype. This is accomplished by transforming the :math:`n` zeros and poles (represented by the input arrays ``_zd`` and ``_pd``) into :math:`n` transformed zeros and poles (represented by the output arrays ``_zdt`` and ``_pdt``). :api:`iirdes_dzpk_lp2bp(*_zd,*_pd,_n,*_zdt,*_pdt)` Converts a low-pass digital prototype :math:`H_d(z)` to a band-pass prototype. This is accomplished by transforming the :math:`n` zeros and poles (represented by the input arrays ``_zd`` and ``_pd``) into :math:`2n` transformed zeros and poles (represented by the output arrays ``_zdt`` and ``_pdt``). .. _section-filter-iirdes-coefficients: Filter Coefficient Computation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The digital filter defined by [eqn-filter-iirdes-Hd] can be expanded to fit the familiar IIR transfer function as in [eqn-filter-iirfilt-Hz]. This can be accomplished using the :api:`iirdes_dzpk2tff()` method. Alternatively, the filter can be written as a set of cascaded second-order IIR filters: .. math:: H_d(z) = G_0 \left[ \frac{1 + z^{-1}} {1 - p_0 z^{-1}} \right]^r \prod_{k=1}^{L} {\left[ G_i \frac{(1-z_iz^{-1})(1-z_i^*z^{-1})} {(1-p_iz^{-1})(1-p_i^*z^{-1})} \right]} :label: eqn-filter-iirfilt-Hd-sos where :math:`r=0` when the filter order is odd, :math:`r=1` when the filter order is even, and :math:`L=(n-r)/2`. This can be accomplished using the :api:`iirdes_dzpk2sosf()` method and is preferred over the traditional transfer function design for stability reasons. .. _section-filter-iirdes-types: Available Filter Types ~~~~~~~~~~~~~~~~~~~~~~ There are currently five low-pass prototypes available for infinite impulse response filter design in |liquid|, as described below. .. toctree:: :caption: Getting Started :maxdepth: 1 groupdelay butter cheby1 cheby2 ellip bessel