tf2sos

Convert digital filter transfer function data to second-order sections form

Syntax

[sos,g] = tf2sos(b,a)

[sos,g] = tf2sos(b,a,order)

[sos,g] = tf2sos(b,a,order,scale)

sos = tf2sos(___)

Description

[sos,g] = tf2sos(b,a) finds a matrix sos in second-order section form with gain g that is equivalent to the digital filter represented by transfer function coefficient vectors b and a.

[sos,g] = tf2sos(b,a,order) specifies the order of the rows in sos.

[sos,g] = tf2sos(b,a,order,scale) specifies the scaling of the gain and numerator coefficients of all second-order sections.

sos = tf2sos(___) embeds the overall system gain in the first section.

Examples

Second-Order Section Implementation of a Butterworth Filter

Design a Butterworth 4th-order lowpass filter using the function butter. Specify the cutoff frequency as half the Nyquist frequency. Implement the filter as second-order sections. Verify that the two representations are identical by comparing their numerators and denominators.

[nm,dn] = butter(4,0.5);
[ss,gn] = tf2sos(nm,dn);
numers = [conv(ss(1,1:3),ss(2,1:3))*gn;nm]

numers = 2×5

    0.0940    0.3759    0.5639    0.3759    0.0940
    0.0940    0.3759    0.5639    0.3759    0.0940

denoms = [conv(ss(1,4:6),ss(2,4:6));dn]

denoms = 2×5

    1.0000    0.0000    0.4860   -0.0000    0.0177
    1.0000    0.0000    0.4860   -0.0000    0.0177

Mass-Spring System

A one-dimensional discrete-time oscillating system consists of a unit mass, $$m$$, attached to a wall by a spring of unit elastic constant. A sensor samples the acceleration, $$a$$, of the mass at $$F_s=5$$ Hz.

Generate 50 time samples. Define the sampling interval $$\Delta t=1/F_s$$.

Fs = 5;
dt = 1/Fs;
N = 50;
t = dt*(0:N-1);
u = [1 zeros(1,N-1)];

The transfer function of the system has an analytic expression:

$$H(z)=\frac{1-z^{-1}(1+\cos \Delta t)+z^{-2}\cos \Delta t}{1-2z^{-1}\cos \Delta t+z^{-2}}.$$

The system is excited with a unit impulse in the positive direction. Compute the time evolution of the system using the transfer function. Plot the response.

bf = [1 -(1+cos(dt)) cos(dt)];
af = [1 -2*cos(dt) 1];
yf = filter(bf,af,u);
stem(t,yf,"o")
xlabel("t")

Compute the time-dependent acceleration using the second-order sections representation of the transfer function to filter the input. Plot the result. The result is the same in both cases.

sos = tf2sos(bf,af);
yt = sosfilt(sos,u);
stem(t,yt,"filled")

Input Arguments

b, a — Transfer function coefficients

vectors

Transfer function coefficients, specified as vectors. Express the transfer function in terms of b and a as

$$H(z)=\frac{B(z)}{A(z)}=\frac{b_1+b_2{z}^{-1}+\cdots +b_{n+1}{z}^{-n}}{a_1+a_2{z}^{-1}+\cdots +a_{m+1}{z}^{-m}}.$$

Example: b = [1 3 3 1]/6 and a = [3 0 1 0]/3 specify a third-order Butterworth filter with normalized 3 dB frequency 0.5π rad/sample.

Data Types: double

order — Row order

'up' (default) | 'down'

Row order, specified as one of the following:

• 'up' — Order the sections so the first row of sos contains the poles farthest from the unit circle.

• 'down' — Order the sections so the first row of sos contains the poles closest to the unit circle.

Data Types: char

scale — Scaling of gain and numerator coefficients

'none' (default) | 'inf' | 'two'

Scaling of gain and numerator coefficients, specified as one of the following:

• 'none' — Apply no scaling.

• 'inf' — Apply infinity-norm scaling.

• 'two' — Apply 2-norm scaling.

Using infinity-norm scaling with 'up'-ordering minimizes the probability of overflow in the realization. Using 2-norm scaling with 'down'-ordering minimizes the peak round-off noise.

Note

Infinity-norm and 2-norm scaling are appropriate only for direct-form II implementations.

Data Types: char

Output Arguments

sos — Second-order section representation

matrix

Second-order section representation, returned as a matrix. sos is an L-by-6 matrix

$$\text{sos}=\left[\begin{array}{cccccc}{b}_{01}& b_{11}& b_{21}& 1& a_{11}& a_{21}\\ b_{02}& b_{12}& b_{22}& 1& a_{12}& a_{22}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ b_{0L}& b_{1L}& b_{2L}& 1& a_{1L}& a_{2L}\end{array}\right]$$

whose rows contain the numerator and denominator coefficients b_ik and a_ik of the second-order sections of H(z):

$$H(z)=g{\displaystyle \prod _{k=1}^L{H}_k}(z)=g{\displaystyle \prod _{k=1}^L\frac{b_{0k}+b_{1k}{z}^{-1}+b_{2k}{z}^{-2}}{1+a_{1k}{z}^{-1}+a_{2k}{z}^{-2}}}.$$

g — Overall system gain

real scalar

Overall system gain, returned as a real scalar.

If you call tf2sos with one output argument, the function embeds the overall system gain in the first section, H₁(z), so that

$$H(z)={\displaystyle \prod _{k=1}^L{H}_k(z)}.$$

Note

Embedding the gain in the first section when scaling a direct-form II structure is not recommended and can result in erratic scaling. To avoid embedding the gain, use tf2sos with two outputs.

Algorithms

tf2sos uses a four-step algorithm to determine the second-order section representation for an input transfer function system:

1. It finds the poles and zeros of the system given by b and a.

2. It uses the function zp2sos, which first groups the zeros and poles into complex conjugate pairs using the cplxpair function. zp2sos then forms the second-order sections by matching the pole and zero pairs according to the following rules:

   1. Match the poles closest to the unit circle with the zeros closest to those poles.

   2. Match the poles next closest to the unit circle with the zeros closest to those poles.

   3. Continue until all of the poles and zeros are matched.

   tf2sos groups real poles into sections with the real poles closest to them in absolute value. The same rule holds for real zeros.

3. It orders the sections according to the proximity of the pole pairs to the unit circle. tf2sos normally orders the sections with poles closest to the unit circle last in the cascade. You can tell tf2sos to order the sections in the reverse order by specifying order as 'down'.

4. tf2sos scales the sections by the norm specified in scale. For arbitrary H(ω), the scaling is defined by

   $${\Vert H\Vert }_p={\left[\frac{1}{2\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\left|H(\omega )\right|}^{p}d\omega }\right]}^{1/p}$$

   where p can be either ∞ or 2. See the references for details on the scaling. The algorithm follows this scaling in an attempt to minimize overflow or peak round-off noise in fixed-point filter implementations.

References

[1] Jackson, L. B. Digital Filters and Signal Processing. 3rd ed. Boston: Kluwer Academic Publishers, 1996.

[2] Mitra, S. K. Digital Signal Processing: A Computer-Based Approach. New York: McGraw-Hill, 1998.

[3] Vaidyanathan, P. P. “Robust Digital Filter Structures.” Handbook for Digital Signal Processing (S. K. Mitra and J. F. Kaiser, eds.). New York: John Wiley & Sons, 1993.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• Any character or string input must be a constant at compile time.

Version History

Introduced before R2006a