3
$\begingroup$

Conventionally, the definition of the system function for a IIR digital system is:

$$H(z)=\frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}+\cdots}{1-a_{1}z^{-1}-a_{2}z^{-2}-\cdots}$$

where coefficients are the ones from the difference equation:

$$y[n]=a_{1}y[n-1]+b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]+\cdots$$

In both equations + and - are operators, and not part of the coefficient themselves. The subtraction operators come from the fact feedback terms where moved on the left side of the equation.


But Scipy definition of function tf2zpk assumes:

$$H(z)=\frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}+\cdots}{a_{0}+a_{1}z^{-1}+a_{2}z^{-2}+\cdots}$$

with only sum operators. Therefore coefficients $\small a_k$, but the first one, have to be redefined with the opposite sign to get the correct result.

The convention will be different, e.g. when calling Numpy convolve. I'm wondering the reason Scipy introduces a different convention for the coefficients, creating a source of confusion and errors. What is the benefit?

$\endgroup$
2
  • $\begingroup$ Because it's based on MATLAB $\endgroup$
    – endolith
    Commented Apr 12, 2023 at 19:19
  • 1
    $\begingroup$ Mathematical notation for a discrete-time LCCDE is $$ \sum_{k=0} a_k y[n-k] = \sum_k b_k x[n-k] $$yielding the standard/conventional transfer function $$ H(z) = \frac{ \sum_k b_k z^{-k}}{ \sum_k a_k z^{-k} } $$all algebraically positive coefficients. However, one can also write an LCCDE as $$y[n] = \sum_{k=1} \tilde{a_k} y[n-k] + \sum_k b_k x[n-k] $$ which is equivalent to $$ y[n] + \sum_{k=1}(-\tilde{a_k}) y[n-k] = \sum_k b_k x[n-k] $$yielding the algebraically negative coefficients in $\tilde{a_k}$. Note that $$ a_k = - \tilde{a_k} $$. $\endgroup$
    – Fat32
    Commented Apr 13, 2023 at 0:13

2 Answers 2

6
$\begingroup$

I agree with the OP's annoyance in that two conventions are used and believe it comes down to what is commonly used in filter design vs what is commonly used in control systems.

However the only difference I see (in the tools) is in the exponents, not the sign of the coefficients as in the OP's first link: tools targeting control system design tend to use transfer functions with polynomials having positive exponents (in s for continuous time or z for discrete time) while those targeting filter design tend to use transfer functions with negative exponents. The two forms for the transfer function commonly used for an example 2nd order IIR system would be:

$$H(z) = \frac{b_0 z^2 +b_1 z + b_2}{ z^2 +a_1 z + a_2}$$

$$H(z) = \frac{b_0 +b_1 z{-1} + b_2z^{-2}}{ 1 +a_1 z^{-1} + a_2z^{-2}}$$

In all the functions I reviewed in the MATLAB, Octave and scipy.signal tools, the transfer functions are given as above with positive coefficients, which as detailed below leads to implementations in the Direct Forms where the feedback coefficient must be negated. (I haven't found one yet, but is there a counter example in the tools where a denominator is given as $1-a_1 z^{-1} -a_2 z^{-2} - \ldots$?)

When the order of the numerator and denominator are the same, the two forms are equivalent as we can divide the numerator and denominator in the top form by $z^2$ to get the bottom form. But if the orders don't match, there will be a difference in delay between the two forms which can lead to erroneous results if the user is not aware.

The function tf2zpk is one example of a transfer function with positive exponents as typical for control system design. The functions freqz, filter in MATLAB and lfilter in scipy.signal are examples using transfer functions with negative exponents as typical for filter design. The control library in Octave and Matlab is also consistent with this convention of having positive exponents. For all functions that are both in MATLAB and Python scipy.signal the convention used for the exponents tend to match function to function. However in all cases, unlike the OP's first link, the convention is to use positive coefficients in the transfer function, leading to negating the feedback coefficients in implementation.

Note that when the difference equation is described as:

$$a_{0}y[n]= -a_{1}y[n-1] - a_{2}y[n-2] +\cdots+b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]+\cdots$$

The resulting z-transform will result in a denominator of $a_0+a_1z^{-1}+a_2z^{-2}+ \ldots$ which I would say is convention in Matlab, Octave, Python and most DSP text books. Below, I show determining the transfer function from the difference equation for the case of a simple 2nd IIR filter to help make this clearer:

$$a_{0}y[n]=-a_{1}y[n-1]-a_{2}y[n-2]+b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]$$

$$a_{0}Y(z) = -a_{1}Y(z)z^{-1} - a_{2}Y(z)z^{-2}+b_{0}X(z)+b_{1}X(z)z^{-1}+b_2X(z)z^{-2}$$

$$a_{0}Y(z) + a_{1}Y(z)z^{-1} + a_{2}Y(z)z^{-2} = b_{0}X(z)+b_{1}X(z)z^{-1}+b_2X(z)z^{-2}$$

$$Y(z)(a_{0} + a_{1}z^{-1} + a_{2}z^{-2}) = X(z)( b_{0}+b_{1}z^{-1}+b_2z^{-2})$$

$$H(z) = \frac{Y(z)}{X(z)} = \frac{b_{0}+b_{1}z^{-1}+b_2z^{-2}}{a_{0} + a_{1}z^{-1} + a_{2}z^{-2}}$$

This sign for the feedback coefficients matches MATLAB's implementation and many text box descriptions for IIR implementation using the Direct Forms as I depict in the graphic below and as shown in "Introduction to Digital Filters for Audio Applications" which is a reputable site by Dr. Julius Smith.

Filter forms

As for convention, this is how I have seen the transfer function for digital filters given (with all the coefficients in the transfer function shown as positive, leading to negative coefficients in the implementation); in my experience it has been contrary to convention to be described as in the first link provided by the OP - but this is very likely as the OP suspects due to MATLAB taking on this form. Either form can be used (as long as it is defined and used consistent to that definition, it is not "wrong"), but a simple Google image search of "Direct Form 1 Filter" confirms with little doubt which form is more common today, and therefore can be called "conventional use".

As for origination, the convention predates MATLAB (although as far as tools go we can certainly refer to it as a MATLAB convention). Matlab originated in 1984 but according to acm.org the MATLAB control toolbox appeared in 1985 and signal processing toolbox (filters) in 1987. The current help for the filter function references Oppenheim, Schafer and Buck “Discrete-Time Signal Processing” published in 1999- however I have this text and confirmed that it uses the opposite convention! Formula 6.27 on page 354 is given as:

$$H(z) = \frac{\sum_{k=0}^M b_k z^{-k}}{1-\sum_{k=1}^N a_k z^{-k}}$$

This Bell System paper from 1964 by R.M. Golden and J.F. Kaiser "Design of Wideband Sampled-Data Filters" also shows an implementation with negative coefficients: Although a is used for the numerator and b for the denominator, we note that positive coefficients are used in the transfer function while negative feedback is applied in the implementation (effectively negating the feedback coefficients). I suspect it has appeared in numerous papers and text books prior to MATLAB's implementations and they went with the viewed as the common convention.

transfer function

block diagram

As a side note, the OP mentioned numpy convolve as not following the same convention, but note that numpy convolve works with polynomials, not ratios of polynomials as these structures depict.

$\endgroup$
4
  • 2
    $\begingroup$ Thanks, you have a good point with MATLAB cloning, but I'm not yet convinced: 1/ whatever the structure, H(z) is equivalent to the conventional one, there is no need to invent a new one for each implementation, 2/ Implementation structures with conventional coefficients are also found, e.g. McClellan, another reputable expert, uses this DF-II and this TDF-II. I'm not saying you're wrong at all, just trying to understand. $\endgroup$
    – mins
    Commented Apr 12, 2023 at 10:56
  • $\begingroup$ @mins looking closer across the control and filter targeted libraries, it seems all the functions assume positive coefficients in the denominator for the transfer functions (which means changing the sign of feedback in the implementation). The difference I see is functions used in control systems tend to use positive exponents, while filter design uses negative. Have you found a counter case in the tools where the denominator is given as 1- a1 z^-1 - a2 z^-2 etc? $\endgroup$ Commented Apr 12, 2023 at 11:27
  • $\begingroup$ @mins I updated my answer to be clearer. As an additional exercise, google "Direct Form 1 Filter" select images: I suspect that negating the feedback coefficients is much more prevalent in the literature (so all positive coefficients in the transfer function form). This has been my experience as well when working in the filter design space. Thanks for the good question! $\endgroup$ Commented Apr 12, 2023 at 12:51
  • 1
    $\begingroup$ Thank you for the detailed explanation. I think you're right, I found other details (see my own answer). There is likely a link between the preferences filter design / control systems and the two different conventions (from Makhoul and from Rabiner/Shafer). I'm selecting your answer anyway. $\endgroup$
    – mins
    Commented Apr 12, 2023 at 12:56
4
$\begingroup$

Puzzled by this problem, I did some research online and found what is likely the reason for this Scipy convention. Dan Boschen provided a first element: Scipy functions are cloned from MATLAB.

So the convention is rather a MATLAB one.

I got this homework from Carnegie Mellon University:

Unfortunately, MATLAB uses the convention used in the classic review article on linear prediction and lattice filtering by John Makhoul (1975).

$$H(z)=\frac{\sum\limits _{\ell=0}^{M}b_{\ell}z^{-\ell}}{1+\sum\limits _{k=1}^{N}a_{k}z^{-k}}$$

This means that the coefficients of the denominator polynomial ak and the reflection coefficients ki produced by MATLAB are opposite in sign to what we have seen in Rabiner and Schafer and used in class discussion.

The paper mentioned is:

Spectral linear prediction: properties and application,
IEEE Trans. Acoust., Speech, Signal Process. 23, 283–296.

So MATLAB implemented this version of the equation rather than:

$$H(z)=\frac{\sum\limits _{\ell=0}^{M}b_{\ell}z^{-\ell}}{1-\sum\limits _{k=1}^{N}a_{k}z^{-k}}$$

which is stated by some authors. Since in Makhoul equation terms are added to 1 in the denominator, and in the other equation terms are subtracted from 1, coefficient for the denominator have opposites values.

I seems most authors use Makhoul convention:

$1+\sum$:

  • MATLAB, Scipy
  • Wikipedia article
  • Makhoul, Manolakis, Steven Smith, Proakis, Oppenheim, Schafer, Buck, Li Tan, Lyons, ...

$1-\sum$:

  • Rabiner, McClellan, Michael Stiber, ...
$\endgroup$
4
  • $\begingroup$ I don't actually think it is a "MATLAB" convention, to check I found one of the papers thought to first reference "Direct Form" for filter design as Malik12 provided at this post: dsp.stackexchange.com/questions/75765/…. There the coefficients are also negative and transfer function positive. I'll update my answer to include that. I suspect many of the text books that predated MATLAB (pre 1984) also had this convention which MATLAB then adopted. $\endgroup$ Commented Apr 12, 2023 at 13:07
  • $\begingroup$ @DanBoschen: I agree, MATLAB adopted the existing definition, at least this is what is stated by Carnegie Mellon. I'm currently looking at who uses which convention, and I'll add some names here. $\endgroup$
    – mins
    Commented Apr 12, 2023 at 13:12
  • $\begingroup$ yes that makes complete sense, as far as tools go-- and how subsequent tools would benefit by copying MATLAB rather than going their own path. All makes sense. $\endgroup$ Commented Apr 12, 2023 at 13:15
  • 1
    $\begingroup$ Good research! I saw MATLAB referenced Oppenheim, Schafer and Buck so had assumed that text was $1+\sum$ but just confirmed in my own copy and it is indeed $1-\sum$ so you can add that one if you like to your $1-\sum$ list. I am concurring with your opinion that MATLAB likely solidified it's current standing to be more commonly seen as $1+\sum$ today. I wish I could Google time machine the equivalent of the web in 1987 right before MATLAB released their signal processing toolbox so can't really be sure anymore what was "convention" back then. $\endgroup$ Commented Apr 12, 2023 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.