Why does scipy introduce its own convention for H(z) coefficients?

Question

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?

Answer 1

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.

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.

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.

Answer 2

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, ...