tf2sos input argument

I have a question about the input argument to the method tf2sos.

The documentation says it takes the coefficients as vectors of the transfer function

.

I suspect that what is actually meant is that it takes as vectors of the transfer function

.

I have tried the following code, but I can not make sense of it.

clear all;

close all;

ts = 1/10000;

z = tf('z', ts);

b0 = 1;

b1 = 2;

b2 = 3;

a0 = 1;

a1 = -2;

a2 = 5;

H = tf([b2, b1, b0], [a2, a1, a0], ts);

Hi = (b0 + b1*z^-1 + b2*z^-2)/(a0 + a1*z^-1 + a2*z^-2);

[sos_H,g_H] = tf2sos(H.Numerator{1}, H.Denominator{1});

[sos_Hi,g_Hi] = tf2sos(Hi.Numerator{1}, Hi.Denominator{1});

K_H = g_H*(sos_H(1,1) + sos_H(1,2)*z^-1 + sos_H(1,3)*z^-2)/(sos_H(1,4) + sos_H(1,5)*z^-1 + sos_H(1,6)*z^-2);

K_Hi = g_Hi*(sos_Hi(1,1) + sos_Hi(1,2)*z^-1 + sos_Hi(1,3)*z^-2)/(sos_Hi(1,4) + sos_Hi(1,5)*z^-1 + sos_Hi(1,6)*z^-2);

figure()

hold on

step(H);

step(K_H, '--');

legend();

grid();

figure()

hold on;

step(Hi);

step(K_Hi, '--');

legend();

grid();

What am I doing wrong?

Answer 1

[Paul]

Hi Nicolai,

As far as I know, discrete time transfer functions in the Signal Processing Toolbox (SPT), whether b/a or sos form, are always represented with coefficients in ascending powers of z^-1 (I'd be very curious to see any doc page that shows otherwise), exactly as shown on the doc page for tf2sos as cited in the question:

.

[Quote from OP in a comment to a since-deleted answer]

"... it doesn't matter which notation (z or z^-1) is used to create the 'b' and 'a' vector. They will have the same entries and create the same transfer function."

That statement is correct ONLY if the order of B(z) is the same as the order of A(z), i.e., m = n, which is clear in that case by multiplying both B(z) and A(z) by z^n (or, more precisely if numel(a) == numel(b)). But it's incorrect otherwise.

Consider an example where that is not the case: H(z) = (0 + 4*z^-1)/(1 + 2*z^-1 + 3*z^-2). For this H(z), we enter b and a for the SPT as

b = [0,4]; % note the leading zero
a = [1,2,3];

Now we convert to sos form

sos = tf2sos(b,a)

sos = 1×6

0 4.0000 0 1.0000 2.0000 3.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

and we see a single section with numerator 0 + 4*z^-1 + 0*z^-2 and denominator 1 + 2*z^-1 + 3*z^-2.

The frequency response of H(z) can be computed with freqz in the SPT, which also assumes b/a (and sos sections) with ascending powers of z^-1

[htf,w] = freqz(b,a);
hsos = freqz([sos;[0,0,1,0,0,1]],w);  % need to append a second section with unity gain when using freqz with a single sos section

The Control Systems Toolbox (CST) allows the tf form to be specified in terms of either z^-1 or z.

Because our original b and a were specified using SPT conventions, we use

Hzi = tf(b,a,-1,'Variable','z^-1')

Hzi =
 
        4 z^-1
  -------------------
  1 + 2 z^-1 + 3 z^-2
 
Sample time: unspecified
Discrete-time transfer function.

We can switch the variable to z if we want

Hz = Hzi; Hz.Variable = 'z'

Hz =
 
       4 z
  -------------
  z^2 + 2 z + 3
 
Sample time: unspecified
Discrete-time transfer function.

and we get the expected result as if we had multiplied the top and bottom of Hzi by z^2.

If we want to get Hz directly, we'd specify the num and den coefficients in descending powers of z, which requires a trailing zero in the numerator to represent 4*z + 0

tf([4,0],[1 2 3],-1)  % default Variable = 'z'

ans =
 
       4 z
  -------------
  z^2 + 2 z + 3
 
Sample time: unspecified
Discrete-time transfer function.

All of these forms have the same frequency response, as they should.

figure
opts = bodeoptions;
opts.PhaseWrapping = 'on';
bode(Hz,Hzi,frd(htf,w,-1),frd(hsos,w,-1),w,opts);

Warning: The sampling time "Ts" was set to +1 (Ts=-1 is ambiguous for FRD models).

Warning: The sampling time "Ts" was set to +1 (Ts=-1 is ambiguous for FRD models).

Summary: The SPT always uses b/a (and sos sections) in ascending powers of z^-1 for discrete-time filters. The CST allows the tf form to be specified in either ascending powers of z^-1 or descending powers of z, with the latter being the default. Make sure to use leading or trailing zeros or both in b and a to specify zero coefficients as needed for any form chosen.

[Paul]

To be safe, one can always specify b and a of the same length using leading and/or trailing zeros as needed, in which case b and a are the same for either z or z^-1 forms of the same transfer function.

Using the example above

b = [0,4,0];
a = [1,2,3];
Hz = tf(b,a,-1)

Hz =
 
       4 z
  -------------
  z^2 + 2 z + 3
 
Sample time: unspecified
Discrete-time transfer function.

Hzi = tf(b,a,-1,'Variable','z^-1')

Hzi =
 
        4 z^-1
  -------------------
  1 + 2 z^-1 + 3 z^-2
 
Sample time: unspecified
Discrete-time transfer function.

Hz - Hzi

ans =
 
  0
 
Static gain.