z-Transform in Scilab
Dominik Kehl
I would be very happy if anybody could explain to me how I get the
z-transform of a fir filter in Scilab and plot its poles and zeros in
an Real/Imag coordinate system.
Many, many thanks,
Dominik
tom_b...@agilent.com
coefficients of an equation like
y = a0*z^0 + a1*z^1 + ... +an*z^n
Then if you simply take roots(p), where p is the polynomial
representing the fir filter, you get the z-domain roots (zeros). Since
it's an fir filter, I trust it will have poles only at the origin?? :-)
You can get the polynomial representation by using the
poly(v,"z","coeff") function.
The result of roots(p) is a set of complex numbers you can plot. See
zgrid() to plot one type of z-plane grid, but you could just plot a unit
circle in the complex plane, too. Just use plot(real(r),imag(r)) for a
simple complex plane plot of complex-valued vector r (except this
is rather ugly because the roots are connected by lines)...or use plot2d
for more control, etc. A much cleaner way to do all this is to convert
the polynomial
to syslin representation and then just use plzr to plot the poles and
zeros. Form an appropriate denominator polynomial using
poly([0,0,0,...,0,1],"z","coeff") where the vector is appropriate
length.
See roots (361), poly (57), zgrid (140), plzr (114) and in general the
manual chapters on polynomial calcs (ch. 8) and plotting (ch. 2). See
also chapters 4 and 7 for some nice utilities to deal with fir and iir
functions, including frequency response plots and the like. (One of the
reasons I like SciLab is the linear system data structure and the
functions that let me deal with such structures easily.)
Example:
fircoeff=[.1,.2,0,-.4,-.6,1,-.6,-.4,0,.2,.1]
firpoly=poly(fircoeff,'z','coeff')
denomcoeff=[zeros(1,size(fircoeff,'c')-1),1]
denompoly=poly(denomcoeff,'z','coeff')
firsl=syslin('d',firpoly,denompoly)
plzr(firsl)
Hope I got that about right!
Cheers,
Tom
In article <3a5a0a73...@news.uni-potsdam.de>,
--
Tom Bruhns -- K7ITM
Sent via Deja.com
http://www.deja.com/
tom_b...@agilent.com
tom_b...@agilent.com wrote:
...
> Example:
> fircoeff=[.1,.2,0,-.4,-.6,1,-.6,-.4,0,.2,.1]
> firpoly=poly(fircoeff,'z','coeff')
> denomcoeff=[zeros(1,size(fircoeff,'c')-1),1]
> denompoly=poly(denomcoeff,'z','coeff')
> firsl=syslin('d',firpoly,denompoly)
> plzr(firsl)
I noticed after posting that, that the 12th order pole which should have
been at the origin showed up as a set of 12 poles around a
small-diameter circle, in the plot produced by plzr. The firsl linear
system representation showed z^12 in the denominator, so I don't quite
understand the plot... A defect perhaps? Comments welcome.
Dominik Kehl
Hopefully somebody can explain it.
Dominik Kehl
wow, many thanks for your comprehensive reply. I am learning a lot
from it. And I will be working on it today.
Thanks a lot,
Dominik
Dominik Kehl
>Consider that the coefficients of an fir filter simply are the
>coefficients of an equation like
>
>y = a0*z^0 + a1*z^1 + ... +an*z^n
>
I think this not correct, because positive time samples are
represented by negative powers of z. So the sequence should be:
y = a0*z^0 + a1*z^-1 + ... +an*z^-n
Please, correct me if I'm wrong.
If the above is correct, I don't know how to use the poly() function
in a way that I obtain the correct result
I hope we can solve that problem,
Dominik
tom_b...@agilent.com
course you can convert the equation I posted to the one you posted by
dividing by z^n. Please note that in the example I gave, that is
exactly what I did, and in fact, if there is excess delay, you can add
more poles at the origin. A pole at the z-plane origin represents a
one-sample delay, and of course is just a factor of "z" in the
denominator of the transfer function for the filter (system). If you
try to use the SciLab "syslin" function with a denominator with lower
order than the numerator, it complains.
I've always preferred to work with z-domain equations with positive
powers of z, so I do that with IIR filters as well. That seems to go
against some conventions, but so long as you account for the
factored-out power of z properly, the math is all correct.
Cheers,
Tom
In article <3a5b7987...@news.uni-potsdam.de>,
--
Dominik Kehl
many thanks for your reply.
tom_b...@agilent.com wrote:
>You are correct (assuming your system cannot be anticipatory ;-), but of
>course you can convert the equation I posted to the one you posted by
>dividing by z^n.
I am sure you're perfectly right, I just cannot get this conversion
done right with my bad mathematical skills.
When I divide your sequence by z^n that operation is equal to
multiplication with z^-n. When I do that (I might do that wrong though
;-) I get something like:
y = a0*z^-n + a1*z-(n-1) + ... + an*z^0
Unfortunately I don't see what denominator I need to obtain
y = a0*z^0 + a1*z^-1 + a2*z^-2 + ... + an*z^-n.
What I am trying to do is the following. I have a simple minimumphase
function:
y = 1 + 0.5*z^-1, that has a zero at z=0.5
and I want to get a plot that shows me that the zero is inside the
unit circle (at 0.5). My plan to use this procedure on more
complicated sequences to see if they are minimumphase or not.
Maybe my example is too theorectic, I don't know.
I would be really, really happy if you would continue helping me with
that and I really appreciate all you did so far,
Dominik
>Please note that in the example I gave, that is
>exactly what I did, and in fact, if there is excess delay, you can add
>more poles at the origin. A pole at the z-plane origin represents a
>one-sample delay, and of course is just a factor of "z" in the
>denominator of the transfer function for the filter (system). If you
>try to use the SciLab "syslin" function with a denominator with lower
>order than the numerator, it complains.
Okay, I think I understand that. Every pole shifts the sequence by one
sample, because it's a multiplication by z^-1.
tom_b...@agilent.com
definitions of a great many of the functions. You will find plzr
defined in the file /macros/xdess/plzr.sci. When you open it and look
at it, you realize that it can take a different set of arguments than a
syslin list, and in fact, it can handle either of two representations of
the syslin list. The one my example feeds to it is a
numerator/denominator transfer function form, and when it sees that, it
calls tr2ss to convert it to a state space form. In that conversion, I
believe that some state space matrices are being created that aren't
very well behaved. That's my guess. But in any event, there is a fix!
If, instead of the example I gave before, you do something like the
following, I believe all will work out fine:
fircoeff=[list of coefficients]
firsize=size(fircoeff,'c');
a=[zeros(1,firsize-1);eye(firsize-2,firsize-2),zeros(firsize-2,1)]
b=[1;zeros(firsize-2,1)]
c=fircoeff(2,firsize)
d=fircoeff(1)
plzr(a,b,c,d)
Note that this is an alternate form for calling plzr which is not
documented in the manual. But from the code, it is obvious that it
works. You could also, after a,b,c,d are defined, do a
firsl2=syslin('d',a,b,c,d)
plzr(firsl2)
and it would be fine. In the case of large fir filters, perhaps this
direct approach to building the state space matrices is better than
asking tf2ss to do it! It certainly gives you matrices which directly
represent how you might actually build the filter, given the obvious
unity scaling in the circular buffer. In the example I ran on my
system, tf2ss didn't give me anything like the system I expected, but
just because it made assumptions about how to represent the system in
state space doesn't mean it's "wrong." It's just a different state
space than I would have chosen. On the other hand, when we have direct
control over how it's represented, and a representation matches how we
would actually build it, there's something to be said for using that
representation. This is even more useful, I think, in working with IIR
filters.
As you may have already noted, there are functions to plot the frequency
response of the filter, too. These are easy to use once you have the
filter in a linear system (state space) form, as the example above does.
See for example repfreq (241).
Hope these ramblings help!
Cheers,
Tom
In article <3a5acd2d...@news.uni-potsdam.de>,
Dominik Kehl
express my sequence in the numerator/denominator form so that
it works in the way you described it.
y = a0*z^0 + a1*z^-1 + a2*z^-2 + ... + an*z^-n is equal to
a0*x^n + a1*x^(n-1) + ... + an
y = ---------------------------------------------
x^n
So if I put my fircoeff-vector in reverse order (btw is there a Scilab
function that does that? I have to look for it) so that I obtain:
fircoeff=[0.5 1]
firpoly=poly(fircoeff,'z','coeff')
denomcoeff=[zeros(1,size(fircoeff,'c')-1),1]
denompoly=poly(denomcoeff,'z','coeff')
firsl=syslin('d',firpoly,denompoly)
firsl =
.5 + z
-------
z
plzr(firsl)
I get the desired plot.
I guess the procedure you describe below does about the same somehow.
I will have a closer look at it tomorrow morning.
Again, many thanks for your input
Dominik
tom_b...@agilent.com
from somewhere they have been properly tested! Maybe that should have
been:
c=fircoeff(2:firsize)
My apologies if there are still other typos in it.
In article <93ihij$fb7$1...@nnrp1.deja.com>,
tom_b...@agilent.com wrote:
> fircoeff=[list of coefficients]
> firsize=size(fircoeff,'c');
> a=[zeros(1,firsize-1);eye(firsize-2,firsize-2),zeros(firsize-2,1)]
> b=[1;zeros(firsize-2,1)]
> c=fircoeff(2,firsize)
> d=fircoeff(1)
> plzr(a,b,c,d)
--
Dominik Kehl
Dominik Kehl
of our thread. Must have been pretty tired last nicht. My apologies
for that, here it is again:
Thank you sooo much, Tom. I just realized (finally) that I have to
express my sequence in the numerator/denominator form so that
it works in the way you described it.
y = a0*z^0 + a1*z^-1 + a2*z^-2 + ... + an*z^-n is equal to
a0*x^n + a1*x^(n-1) + ... + an
y = ---------------------------------------------
x^n
So if I put my fircoeff-vector in reverse order (btw is there a Scilab
function that does that? I have to look for it) so that I obtain:
fircoeff=[0.5 1]
firpoly=poly(fircoeff,'z','coeff')
denomcoeff=[zeros(1,size(fircoeff,'c')-1),1]
denompoly=poly(denomcoeff,'z','coeff')
firsl=syslin('d',firpoly,denompoly)
tom_b...@agilent.com
ke...@rz.uni-potsdam.de wrote:
...
> So if I put my fircoeff-vector in reverse order (btw is there a Scilab
> function that does that? I have to look for it)
If you do something like
j=j(size(j,'c'):-1:1)
it will reverse the order of the row-vector j. You should read the
"Introduction to Scilab" to get a whole lot of hints and insights into
things like this!
...
> I guess the procedure you describe below does about the same somehow.
> I will have a closer look at it tomorrow morning.
(I hope what I'm posting is of some use to some other "lurkers"! I know
that there are people who read this forum who are far, far beyond me
mathematically, but maybe there are others who are just trying to use
the basic tools...If this is too specific an interest to be posting
here, we could probably take it off-line.)
The way I showed (the second time) is different that what you suggest.
It puts your filter into state-space notation, where x is a vector of
the values held in the shift register (or circular buffer) of an FIR
implementation, u is the input signal, and y is the output signal. The
state space representation for a discrete time system is in general
something like:
x(k+1) = A x(k) + B u(k)
y(k) = C x(k) + D u(k)
Note that you could also say
z x = A x + B u
y = C x + D u
where u, x and y are functions of z. And for a continuous time system,
you can use the same notation, except in the s-domain, and the
time-domain equivalent of the first equation becomes
dx/dt = Ax + Bu
There are a lot of functions in SciLab that let you work easily with
state-space (or linear system) representations. Very handy! For
example, you can cascade two filters just by "multiplying" their syslin
representations. You can easily simulate the effects of feedback in
like manner, and can sum the outputs of two. See section 2.7 of the
introduction.
Here's what my (second) little code segment, posted yesterday, does.
For a simple single-input, single-output system, for an FIR filter, the
A matrix implements the shift register:
A = [ 0 0 0 0 ... 0
1 0 0 0 ... 0
0 1 0 0 ... 0
...
0 0 0 ... 1 0 ]
And B represents introducting the input into the shift register:
B = [ 1
0
...
0 ]
And C and D represent the FIR coefficients to sum the shift register
values and the input, each properly scaled:
C = [ a1 a2 a3 ... an ]
D = [ a0 ]
(Or the reverse, depending on how you defined the order!)
My little code segment just builds the A, B, C and D matrices, knowing
the vector of FIR coefficients as an input. Then plzr uses A, B, C and
D to find the poles and zeros of the system, and plot them.
You could just as well take the roots of your denominator and roots of
your numerator and plot them on an x-y plane as real,imaginary pairs,
and plot a unit circle for reference. The code for plzr shows you just
how to do that.
Cheers,
Tom
