CPS - 1 - 2006
ZASTOSOWANIA PRZEKSZTAA
CENIA ZET
Rozwiązywanie równań różnicowych
Dyskretny system liniowy-stacjonarny można opisać równaniem różnicowym w
postaci ogólnej
NM
[ ] [ ]
"a y n - k = "b x n - k
kk
k =0 k =0
Przekształcenie ZET pozwala w efektywny sposób obliczać odpowiedzi systemu na
zadane wymuszenia.
UPrzykład:
Rozważmy przyczynowy LTI system opisany równaniem różnicowym
y n 0.9y n -1 = x n
[ ]- [ ] [ ]
Wyznaczymy odpowiedz
systemu na wymuszenie jednostkowe x[n]=1[n] dla
wartości początkowej odpowiedzi y[-1]=2.
Wykorzystamy własności przesunięcia w dziedzinie czasu dla transformaty
jednostronnej
Z Z -1
Jeżeli y n �!Ż# z to y n -1 �!Ż# z Y z + y
Y
[ ] ( ) [ ] ( ) [-1
]
Przekształcamy ( transformacja Z ) obie strony równania różnicowego
-1
Y z - 0.9 z Y z + y = X z
( ) ( ) [-1
] ( )
( )
CPS - 2 - 2006
-1
Y z 1- 0.9z = X z + 0.9y
( ) ) ( ) [-1
( ]
X z 0.9y
( ) [-1
]
Y z =+
( )
-1 -1
1- 0.9z 1- 0.9z
Ponieważ transformata skoku jednostkowego wynosi
1
X z =
( )
-1
1- z
oraz wartość początkowa
y[-1]=2
to:
11.8
Y z =+
( )
-1 -1 -1
1- 0.9z 1- z 1- 0.9z
()( )
Stosujemy rozkład na ułamki proste:
AA1.8
12
Y z =+ +
( )
-1 -1 -1
1- 0.9z 1- z 1- 0.9z
Współczynniki ułamków prostych obliczamy z zależności
11
A = = = -9
1 z=0.9
-1
10
1- z 1-
9
11
A = = =10
2 z=1
-1
9
1- 0.9z 1-
10
Stąd
-9101.8
Y z =+ + OZ z >1
( )
-1 -1 -1
1- 0.9z 1- z 1- 0.9z
CPS - 3 - 2006
Jeżeli system jest przyczynowy to oznacza, że odpowiedz
na skok jednostkowy jest
funkcją prawostronną, dla czasów dodatnich. Obszar zbieżności będzie zatem
zewnętrzem okręgu o promieniu 1.
Odwrotna transformata dla przyjętego obszaru zbieżności daje rozwiązanie w postaci
sygnału wyjściowego systemu:
nn
y n =-9 0.9 1 n +10 �"1 n +1.8 0.9 1 n
[ ] ( ) [ ] [ ] ( ) [ ]
n
Ą#10 ń#
y n = - 8.2 �" 0.9 �"1 n
[ ] ( ) [ ]
Ł#Ś#
UPrzykład
Wyznaczymy odpowiedz
impulsową systemu przyczynowego LTI opisanego
równaniem różnicowym, dla zerowych warunków poczatkowych.
y n +1 + 2y n = x n
[ ] [ ] [ ]
Wykorzystamy własności:
Z
y n +1 �!Ż# zY z - zy 0

[ ] ( ) [ ]
Z
� n �!Ż#
1
[ ]
Transformując obie strony równania różnicowego ( przy wyznaczaniu odpowiedzi
impulsowej zakłada się zerowe warunki początkowe ) otrzymamy:
zY z - zy 0 + 2Y z =1
( ) [ ] ( )
zY z + 2Y z =1
( ) ( )
Y z z + 2 =1
( )( )
CPS - 4 - 2006
Ostatecznie transformata odpowiedzi impulsowej wynosi:
1
Y z =
( )
z + 2
1 1
= �"
-1
z 1+ 2z
1
�#ś#
-1
= z
ś#1+ 2z-1
ź#
�# #
Transformata odwrotna wyrażenia w nawiasie (obszar zbieżności jest zewnętrzem
okręgu o promieniu 2) wynosi
1 n
Z
�!Ż#
-2 1 n
( ) [ ]
-1
1+ 2z
-1
P
Mnożenie przez zP powoduje w dziedzinie czasu opóz
nienie sygnału o jedną próbkę,
stąd otrzymujemy:
1 n-1
�#ś#�!Ż#
-1 Z
z -2 1 n -1
( ) [ ]
ś#1+ 2z-1
ź#
�# #
Odpowiedz
impulsowa systemu wynosi
n-1
y n =
) [ ]
[ ] (-2 1 n -1
Transmitancja systemu dyskretnego
UDefinicjaU: Transmitancję H(z) systemu LTI definiujemy jako transformatę ZET
odpowiedzi impulsowej systemu h[n].
Odpowiedz
y[n] systemu na dowolne wymuszenie x[n] oblicza się jako splot
odpowiedzi impulsowej h[n] systemu oraz wymuszenia:
y n = h n " x n
[ ] [ ] [ ]
CPS - 5 - 2006
Przekształcając obie strony wyrażenia oraz wykorzystując własność
transformaty ZET ( splotu w dziedzinie czasu ), możemy wyrazić transformatę
odpowiedzi Y[z] w postaci iloczynu transmitancji H[z] oraz transformaty wymuszenia
X[z].
Y z = H z X z
( ) ( ) ( )
Stąd transmitancja:
Y z
( )
H z =
( )
X z
( )
Transmitancja systemu w zależności od współczynników równania różnicowego
System LTI w dziedzinie czasu opisuje równanie różnicowe.
NM
[ ] [ ]
"a y n - k = "b x n - i
ki
k =0 i=0
Wykorzystując własność liniowości przekształcenia ZET oraz transformaty
przesuniętych w czasie sygnałów otrzymuje się:
NM
-k -i
( ) ( )
"a z Y Z = "b z X z
ki
k =0 i=0
Stąd transmitancja systemu może być opisana jako:
M
-i
"bz
i
Y z
( )
i=0
H z = H z =
( ) ( )
N
X z
( )
"a z-k
k
k =0
CPS - 6 - 2006
UPrzykład:
Znajdziemy równanie różnicowe opisujące system jeżeli dana jest jego transmitancja
5z + 2
H z =
( )
z2 + 3z + 2
-2
P
Pomnożymy licznik i mianownik przez zP i porównamy z zależnością ogólną
dla H(z)
b0 + b1z-1 + b2z-2
H z =
( )
a0 + a1z-1 + a2z-2
5z-1 + 2z-2
H z =
( )
1+ 3z-1 + 2z-2
Stąd wynika równanie różnicowe opisujące system ma postać:
a0 y n + a1y n -1 + a2 y n - 2 = b1x n -1 + b2x n - 2
[ ] [ ] [ ] [ ] [ ]
y n + 3y n -1 + 2y n - 2 = 5x n -1 + 2x n - 2
[ ] [ ] [ ] [ ] [ ]
Transmitancja systemu w zależności od macierzy stanu
System liniowy-stacjonarny w dziedzinie czasu opisują równania stanu
Q n +1 = AQ n + bx n
[ ] [ ] [ ]
y n = cQ n + Dx n
[ ] [ ] [ ]
CPS - 7 - 2006
Transformata Z wektora stanu wynosi:
Q1 z
Ą#ń#
( )
ó#Q z
( )Ą#
2
ó#Ą#
Q z =
( )
ó#Ą#

ó#
( )Ą#
ó#Q z Ą#
N
Ł#Ś#
Transformata pierwszego równania macierzowego wynosi
zQ z = AQ z + bX z
( ) ( ) ( )
-1
Q z = zI - A bX z
( ) ( ) ( )
Dla drugiego równania mamy
Y z = cQ z + DX z
( ) ( ) ( )
Podstawiając pierwsze do drugiego:
Ą#c -1
Y z = zI - A b + Dń# X z
( ) ( ) ( )
Ł#Ś#
Stąd transmitancja systemu z zależności od macierzy stanu A, b, c, D:
Y z
( )
H z =
( )
X z
( )
-1
H z = c zI - A b + D
( ) ( )
CPS - 8 - 2006
UPrzykład
Wyznaczymy transmitancję oraz równanie różnicowe systemu LTI, opisanego
w przestrzeni stanu:
0 1 0
Ą# ń# Ą# ń#
A = , b = c = 3 0 , D 0
[ ] [ ]
ó# ó#2Ą#,
Ł#-1 1Ą# Ł# Ś#
Ś#
Obliczamy:
1 0 0 1 z
Ą# ń# Ą# ń# Ą# -1
ń#
zI - A = z =
ó#0 1Ą# - ó# ó#1 z -1Ą#
Ł# Ś# Ł#-1 1Ą# Ł# Ś#
Ś#
z
ń#
-1 1 Ą# -1 1
zI
( - A =
)
z2 - z +1ó# -1 zĄ#
Ł# Ś#
Transmitancja systemu
-1
H z = c zI - A b + D
( ) ( )
z -1 1 0
1 ń#Ą# ń#
3 0
H z =+ 0
( ) [ ]Ą#
ó#
z2 - z +1 -1 zĄ#ó#2Ą#
Ł# Ś#Ł# Ś#
2
1
= 3 0
[ ]Ą# ń#
ó#2zĄ#
z2 - z +1
Ł# Ś#
6
=
z2 - z +1
-2
Mnożymy licznik i mianownik przez zP P
6z-2
H z =
( )
1- z-1 + z-2
Równanie różnicowe:
y n y n -1 + y n - 2 = 6x n - 2
[ ]- [ ] [ ] [ ]
CPS - 9 - 2006
UWykorzystanie funkcji Matlab a
Charakterystyki systemów LTI
MATLAB zawiera szereg specjalnych funkcji pozwalających na swobodne
przechodzenia między tymi różnymi opisami systemu.
Jeżeli wektory b i a zawierają współczynniki wielomianów funkcji
transmitancji odpowiednio licznika i mianownika, to funkcja tf2ss(b,a) określa opis
systemu w przestrzeni stanu, natomiast tf2zp(b,a) przekształca funkcję transmitancji
do postaci zero-biegunowej.
Podobnie zp2ss i zp2tf konwertują opis zero-biegunowy do opisów
odpowiednio w przestrzeni stanu oraz transmitancji.
Funkcje ss2tf i ss2zp to konwersje opisu w przestrzeni stanu do opisu w
postaci transmitancji i zero-biegunowej.
UPrzykład
Rozpatrzymy system dyskretny opisany funkcją transmitancji
0.094 z4 + 4z3 + 6z2 + 4z +1
()
H z =
( )
z4 + 0.4860z2 + 0.0177
Matlab
b=.094*[1 4 6 4 1];
a=[1 0 0.486 0 0.0177];
zplane(b,a)
CPS - 10 - 2006
1
0.8
0.6
0.4
0.2
4
0
-0.2
-0.4
-0.6
-0.8
-1
-1 -0.5 0 0.5 1
Real Part
System w postaci zero-biegunowej:
>> [r,p]=tf2zp(b,a)
r =
-1.0002
-1.0000 + 0.0002i
-1.0000 - 0.0002i
-0.9998
p =
0 + 0.6681i
0 - 0.6681i
0 + 0.1991i
0 - 0.1991i
Imaginary Part
CPS - 11 - 2006
System opisany w przestrzeni stanu:
>> [A,B,C,D ]=tf2ss(b,a)
A =
0 -0.4860 0 -0.0177
1.0000 0 0 0
0 1.0000 0 0
0 0 1.0000 0
B =
1
0
0
0
C =
0.3760 0.5183 0.3760 0.0923
D =
0.0940
Charakterystyka częstotliwościowa:
[H,w]=freqz(b,a,250);
plot(w,abs(H))
Widmo amplitudowe
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5
Czestotliwosc
System posiada zero ( o krotności 4) w punkcie z=-1 oraz cztery bieguny na osi
urojonych. Takie położenie zer implikuje, bardzo małe wartości na charakterystyce
amplitudowej dla wysokich częstotliwości.
Amplituda
CPS - 12 - 2006
Przyczynowość systemów LTI
Odpowiedz
impulsowa systemu przyczynowego jest zerowa dla n<0.
Dlatego odpowiedz
impulsowa systemu przyczynowego jest zdeterminowana
przez prawostronną transformatę odwrotną transmitancji.
Patrząc na bieguny transmitancji, to wszystkie muszą mieć moduł mniejszy niż
promień obszaru zbieżności, inaczej obszar zbieżności transmitancji jest na
zewnątrz dla wszystkich biegunów.
Stabilność systemów LTI
Wszystkie bieguny generują składniki wykładnicze (patrz: rozkład na ułamki
proste):
malejące, gdy bieguny znajdują się wewnątrz okręgu jednostkowego (sy)
b=[0 1];
a=[1 0.5];
x=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
zplane(b,a);
y=filter(b,a,x);
stem(y);
1
1
0.8
0.6
0.4
0.5
0.2
0
-0.2
0
-0.4
-0.6
-0.8
-1
-0.5
-1 -0.5 0 0.5 1 0 2 4 6 8 10 12 14 16 18
Real Part
Imaginary Part
CPS - 13 - 2006
rosnące gdy bieguny znajdują się na zewnątrz okręgu jednostkowego.
b=[0 1];
a=[1 -1.5];
x=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
zplane(b,a);
y=filter(b,a,x);
stem(y);
700
1
0.8
600
0.6
500
0.4
0.2
400
0
300
-0.2
-0.4
200
-0.6
-0.8 100
-1
0
-1 -0.5 0 0.5 1 1.5 0 2 4 6 8 10 12 14 16 18
Real Part
System jest stabilny (BIBO) jeżeli odpowiedz
impulsowa jest bezwzględnie
sumowalna, dotyczy to przypadku gdy odpowiedz
impulsowa systemu jest suma
malejących składników. Patrząc na transmitancję, oznacza to, że jej obszar
zbieżności musi zawierać okrąg jednostkowy.
Wnioski:
System jest stabilny i przyczynowy jeżeli wszystkie bieguny są położone
wewnątrz okręgu jednostkowego.
System przyczynowy, stabilny jest minimalnofazowy jeżeli wszystkie zera i
bieguny transmitancji są położone wewnątrz okręgu jednostkowego
Imaginary Part
CPS - 14 - 2006
Przykład:
Transmitancja systemu stabilnego wynosi
z-1
H z =
( )
1
1+ 2z-1 1- z-1
()()
3
Oblicz odpowiedz
impulsową systemu.
Odpowiedz
impulsowa to odwrotna transformata ZET transmitancji. Rozkład na
ułamki proste:
AB
H z =+
( )
1
1+ 2z-1 1- z-1
() ()
3
-1
1
z - 13
2
A = = =- =-
-1
11 1 1
1- z 1- �" 2 1+ 7
() (- ) ( )
33 2 6
z=-2
-1
z 3 3
B == =
-1
1+ 2z 1+ 2 �"3 7
()
1
z=
3
33
-
77
H z =+
( )
-1 -1
1
1+ 2z 1- z
() ()
3
Wybieram dla układu stabilnego OZ jako pierścień miedzy okręgami o promieniach 2
i 1/3 ( dlatego, ten obszar zawiera okrąg jednostkowy).
3
- n
7
3

(-2 1 -1
) [-n
]
7
-1
1+ 2z
()
3
n
7
3 1
1 n
( ) [ ]
7 3
-1
1
1- z
()
3
Odpowiedz
impulsowa:
nn
33 1
h n =
) [-n
[ ] (-2 1 -1 + 1 n
] ( ) [ ]
77 3
CPS - 15 - 2006
Matlab - schematy blokowe
MATLAB pozwala na konwersję modelu  zero-biegunowego lub opisanego
w przestrzeni stanu na kaskadowe połączenie układów drugiego rzędu postaci
opisanych jako:
b0 + b1z-1 + b2z-2
H z =
( )
1+ a1z-1 + a2z-2
W dziedzinie czasu taki system opisuje równanie różnicowe
Ą# ń# Ą# ń# Ą# ń# Ą# ń#
y Ą#nń# + a1y n -1Ś# + a2 y n - 2Ś# = b0x Ą#nń# + b1x n -1Ś# + b2x n - 2Ś#
Ł# Ś# Ł# Ł# Ł# Ś# Ł# Ł#
Schemat blokowy, który odpowiada powyższym równaniom:
X(z) b0 Y(z)
Ł Ł
z-1
b1
-a1
Ł Ł
z-1
-a2 b2
CPS - 16 - 2006
Współczynniki kaskadowo połączonych układów 2 rzędu otrzymuje się
stosując funkcję MATLABA zp2sos lub ss2sos. Format funkcji jest nastepujący:
>> sos=zp2sos(z,p,k)
gdzie z- wektor zawierający zera transmitancji, p- wektor zawierający bieguny
transmitancji, k- wzmocnienie. W wyniku otrzymuje się macierz 6 kolumnową o
liczbie wierszy zależnej od liczby wyznaczonych kaskad. Kolejne elementy każdego
wie
rsza zawierają współczynniki kolejno b0, b1, b2, 1, a1, a2.
UPrzykład:
(1+ jz-1)(1- jz-1)(1+ z-1)
H z =
( )
jĄ /4 jĄ /4 jĄ /8 jĄ /8
11 33
(1- e z-1)(1- e- z-1)(1- e z-1)(1- e- z-1)
22 44
1 3
System posiada zera dla z=ąj oraz z=-1, bieguny dla z = eą jĄ / 4 oraz z = eą jĄ / 8 .
2 4
Zastosujemy funkcję zp2sos:
z=[-1 -j j];
p=[.5*exp(j*pi/4), .5*exp(-j*pi/4), .75*exp(j*pi/8), .75*exp(-
j*pi/8)];
k=1;
sos=zp2sos(z,p,k);
sos =
0 1.0000 1.0000 1.0000 -0.7071 0.2500
1.0000 0 1.0000 1.0000 -1.3858 0.5625
CPS - 17 - 2006
Wyznaczone kaskady 2 rzędu mają następujące transmitancje:
z-1 + z-2
H1(z)=
1- 0.7071z-1 + 0.25z-2
1+ z-2
H2(z)=
1-1.3858z-1 + 0.5625z-2
Na tej podstawie można narysować schemat blokowy systemu.
X(z) Y(z)
Ł Ł Ł
z-1 z-1
0.0171
1.3858
Ł Ł Ł
z-1 z-1
-0.25
-0.5625
Przykład:
Zadanie polega na doborze systemu, który skutecznie będzie tłumił sygnał
wejściowy sinusoidalny. Wówczas z charakterystyki amplitudowej tego filtru będzie
można odczytać częstotliwość tłumionego sygnału.
CPS - 18 - 2006
Prostym systemem, który może pełnić funkcję filtru wycinającego
pojedynczą składową sinusoidalną jest filtr przedstawiony na rysunku:
Xn-1 Xn Xn+1
Z-1 Z-1
a a
+
en
Sygnał na wyjściu filtru opisuje równanie :
en = axn-1 + xn + axn+1
Transformata  Z równania różnicowego:
E z = X z 1+ a z-1 + zĄ#
( ) ( )Ą#
( )ń#
ó#
Ł# Ś#
Transmitancja filtru wynosi:
H z =1+ a z-1 + z
( )
( )
j�
Podstawiając z = e oblicza się charakterystykę częstotliwościową filtru:
H j� =1+ a e- j� + ej�
( ) ( )
CPS - 19 - 2006
H j� =1+ 2acos�
( )
Charakterystyka amplitudowa:
H j� = 1+ 2acos�
( )
1000
1
0.9
0.8
0.7
2
0.6
0.5
0.4
1
0.3
0.2
2/3
0.1
11/20
101/200
0
0 pi/2 pi
Filtr zerowy dla a=[1000 2 1 2/3 11/20 101/200];
Miejsce zerowe charakterystyki amplitudowej
1+ 2acos� = 0
odpowiada częstotliwości tłumionego sygnału:
1
ś#
� = arccos�# -
ś# ź#
2a
�# #
CPS - 20 - 2006
Dla częstotliwości próbkowania fp częstotliwość sygnału wyniesie:
fp fp
1
ś#
f = � = arccos�# -
ś# ź#
2Ą 2Ą 2a
�# #
Pozostaje zatem wyznaczyć współczynnik  a filtru.
Założymy, opisany filtr przetwarza sygnał wejściowy o długości M+1 próbek.
Chcemy aby filtr był tak dobrany, aby na jego wyjściu nic się nie pojawiało.
Jednak ze względu na szum w sygnale wejściowym i błędy numeryczne
wystarczy jeżeli będzie spełniony warunek minimum średniokwadratowego:
M -1
"
en 2 = 0
( )
"
"a
n=1
Podstawiając:
M -1
2
"
a xn-1 + xn+1 + xn = 0
[]
()
"
"a
n=1
M -1
" 2
2
[] []
"(a xn-1 + xn+1 + 2axn xn-1 + xn+1 + ...)= 0
"a
n=1
M -1 M -1
2
[][]
"2a xn-1 + xn+1 + "2x xn-1 + xn+1 = 0
n
n=1 n=1
M -1 M -1
2
a xn-1 + xn+1 = - xn xn-1 + xn+1
[] []
""
n=1 n=1
CPS - 21 - 2006
Współczynnik filtru:
M -1
- xn xn-1 + xn+1
[]
"
n=1
a =
M -1
2
xn-1 + xn+1
[]
"
n=1
Stąd ostatecznie częstotliwość sygnału:
M -1
�# 2 ś#
xn-1 + xn+1 ź#
[]
"
ś#
fp
n=1
ś#ź#
f = arccos
M -1
2Ą
ś#
2 xn xn-1 + xn+1 ź#
[]ź#
"
ś#
�# n=1 #
Wyniki symulacji:
M=40
51.5
Estim.
True
50
49.5
0
0.3 0.6
time (s)
U*)Przygotowano na podstawie: S. Haykin, B.Van Veen  Signals and System New
York 1999
frequency (Hz)