z transfer function evaluates to inaccurate value near z=1

I have a simple 6th order CIC (sinc2) filter made by follwoing commands:

z = tf('z',1);

tf = ((1-z^(-2))/(1-z^(-1)))^6;

This Transfer Function is known to have a magnitude = 64 at z=1.

I am trying to evaluate the magnitude of this transfer function using evalfr command as follows:

abs(evalfr(tf,exp(-2i*pi*1/10000)))

It evaluates to 8.0001 which is completely wrong.

However if I evaluate this function at a value slightly away from z=1, then it gives accurate answer:

abs(evalfr(tf,exp(-2i*pi*1/100))) >>> 63.8108

I also tried to plot the magnitude response of this filter using Filter Visualization Tool (fvtool). Even here I see the magnitude response very inaccurate when f is close to 0 (which corresponds to the case when z is close to 1).

regards,

Ashish

Answer 1

[David Goodmanson]

Hi Ashish,

This function,

((1-1/z^2)/(1-1/z))^6

is identically equal to

((z+1)/z)^6

so you can just use that. And it does equal 64 at z = 1.

[Ashish Panpalia]

Thanks David,

That solves the problem.

So basically we need to cancel the zeros and poles before going for it's evaluation. Otherwise I guess there would be accuracy problems.

I also saw something mentioned about accuracy for system representation in the following link:

https://in.mathworks.com/help/control/ug/using-the-right-model-representation.html

For this expression it was easy to solve on paper or simplify it. But is there any command in matlab which simplifies the transfer function by cancelling poles and zeros?

I tried minreal, but somehow it's not working

regards,

Ashish

Answer 2

[Paul]

Hi Ashish,

This problem illustrates a couple of things to keep in mind when using the Control System Toolbox

Here is the original code

format long e
z = tf('z',1);

However, we shouldn't use tf as a variable name because tf is a very important CST function that we don't want to shadow with a variable name

h = ((1-z^(-2))/(1-z^(-1)))^6

h =
 
   z^18 - 6 z^16 + 15 z^14 - 20 z^12 + 15 z^10 - 6 z^8 + z^6
  -----------------------------------------------------------
  z^18 - 6 z^17 + 15 z^16 - 20 z^15 + 15 z^14 - 6 z^13 + z^12
 
Sample time: 1 seconds
Discrete-time transfer function.

When using tf as above, we get very high order polynomials. The denominator polynmial has at least one root exactly at z = 1

polyval(h.den{:},1)

ans = 
     0

So bad things will happen when evaluating h at z = 1

evalfr(h,1)

ans = 
   Inf

However, we also know that the poles at z = 1 should cancel with zeros at z = 1. So we can try minreal

h = minreal(h)

h =
 
  z^12 + 1.688e-14 z^11 - 6 z^10 - 7.816e-14 z^9 + 15 z^8 + 1.563e-13 z^7 - 20 z^6 - 1.457e-13 z^5 + 15 z^4 + 6.395e-14 z^3 - 6 z^2 - 1.354e-14 z + 1
  ---------------------------------------------------------------------------------------------------------------------------------------------------
                                                z^12 - 6 z^11 + 15 z^10 - 20 z^9 + 15 z^8 - 6 z^7 + z^6
 
Sample time: 1 seconds
Discrete-time transfer function.

However, minreal doesn't work very well because root finding of high-order polynomials with multiple, repeated roots is hard (or so I've been led to believe).

evalfr(h,1)

ans = 
   Inf

Alternatively, we could have used minreal from the start

h = (minreal((1-z^(-2))/(1-z^(-1))))^6

h =
 
  z^6 + 6 z^5 + 15 z^4 + 20 z^3 + 15 z^2 + 6 z + 1
  ------------------------------------------------
                        z^6
 
Sample time: 1 seconds
Discrete-time transfer function.

And then we get the expected result

evalfr(h,1)

ans = 
    64

Genereally speaking, the zpk form is preferred over the tf form (becuase of numerlcal robustess concerns when dealing with polynomials)

z = zpk('z',1);
h = ((1-z^(-2))/(1-z^(-1)))^6

h =
 
  z^6 (z+1)^6 (z-1)^6
  -------------------
     z^12 (z-1)^6
 
Sample time: 1 seconds
Discrete-time zero/pole/gain model.

This actually looks pretty good. Unfortunately, the display is misleading, as the poles (and zeros) are not exactly at z = 1.

h.p{:}.'

ans = 1×18

0 0 9.999999999999997e-01 0 0 9.999999999999997e-01 0 0 9.999999999999997e-01 0 0 9.999999999999997e-01 0 0 9.999999999999997e-01 0 0 9.999999999999997e-01

I'm a little bit surprised by this result, but I'm sure it's a result of how the CST has to do the algebraic manipulations in z

Consequently, evalfr still fails (and maybe this result is worse than the tf case because it's not obviously wrong)

evalfr((h),1)

ans = 
     8.779149519890415e-02

However, the poles and zeros of h are close enough to unity that minreal works as desired

h = minreal(h)

h =
 
  (z+1)^6
  -------
    z^6
 
Sample time: 1 seconds
Discrete-time zero/pole/gain model.

evalfr(h,1)

ans = 
    64

Of course, we could have used minreal from the outset

h = (minreal((1-z^(-2))/(1-z^(-1))))^6

h =
 
  (z+1)^6
  -------
    z^6
 
Sample time: 1 seconds
Discrete-time zero/pole/gain model.

evalfr(h,1)

ans = 
    64

Also, it's worth pointing out that even this approach is not without issue. To wit:

h = minreal((1-z^-2)/(1-z^-1))

h =
 
  (z+1)
  -----
    z
 
Sample time: 1 seconds
Discrete-time zero/pole/gain model.

h.P{:}

ans = 
     0

h.Z{:}

ans = 
    -9.999999999999999e-01

So, even with this simple case, the numerical zero isn't exactly at -1. Continuing

h = h^6;
h.P{:}.'

ans = 1×6

0 0 0 0 0 0

h.Z{:}.'

ans = 1×6

-9.999999999999999e-01 -9.999999999999999e-01 -9.999999999999999e-01 -9.999999999999999e-01 -9.999999999999999e-01 -9.999999999999999e-01

isequal(evalfr(h,1),64)

ans = logical
   1

I guess the rounding works in our favor and we get the exact, expected result.

The best approach would be to start with

h = zpk(-1,0,1,1)

h =
 
  (z+1)
  -----
    z
 
Sample time: 1 seconds
Discrete-time zero/pole/gain model.

h.Z{:}

ans = 
    -1

Then

h = h^6

h =
 
  (z+1)^6
  -------
    z^6
 
Sample time: 1 seconds
Discrete-time zero/pole/gain model.

h.Z{:} % zeros exactly at z = -1

ans = 6×1

-1 -1 -1 -1 -1 -1

isequal(evalfr(h,1),64)

ans = logical
   1

Summary: avoid the tf form if possible (the doc talks about this), simplify earlier rather than later, avoid or minimize "transfer function algebra" if possible.

[William Rose]

@Ashish Panpalia, I learned a lot from @Paul's answer.

[Ashish Panpalia]

I agree with you @William Rose. The explanation from @Paul has uncluttered lots of things. Thanks a lot @Paul

Answer 3

[William Rose]

@Ashish Panpalia,

The function is tricky, as z approaches 1, since the numerator and denominator both tend toward zero.

Why do you think G(z=1)=64?

which goes toward as z tends toward +1.

L'Hopital's rule predicts something similar: d(numerator)/dz is finite and positive, d(denominator)/dz goes to zero, so G(z) goes to infinity as z goes toward +1.

You tried approaching the point z=1 from below on the unit circle, with values such as z=exp(-i*(2*pi/100)), z=exp(-i*(2*pi/1e3)), etc. The strange results you observed are due to round off error as the numerator and denominator both approach zero. If you approach z along the real axis, from more positive values, or along the real axis from less positive values, you also get wrong answers, but they are quite different wrong answers. It just shows that the function G(z) is not nice near z=1.

[Ashish Panpalia]

Hello William,

the power of 6 is applied to both Numerator and Denominator. Sorry, the expression in my question is not very readable. So the transfer function looks like :

This one resolves to nice 64 when you apply L'Hopital's rule till the 6th derivative.

regards,

Ashish