z transfer function evaluates to inaccurate value near z=1

20 views (last 30 days)
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

Accepted Answer

David Goodmanson
David Goodmanson on 4 Nov 2022
Edited: David Goodmanson on 4 Nov 2022
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.
  1 Comment
Ashish Panpalia
Ashish Panpalia on 4 Nov 2022
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

Sign in to comment.

More Answers (2)

Paul
Paul on 13 Nov 2022
Edited: Paul on 13 Nov 2022
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
William Rose on 4 Nov 2022
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.
  1 Comment
Ashish Panpalia
Ashish Panpalia on 4 Nov 2022
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

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!