il filtro interpolatore

Tecniche avanzate di controllo 12/13
Applicazione della decomposizione serie parallelo:
il filtro interpolatore
Ruzza Marco - Vidali Elia
May 11, 2013
In questo documento trattiamo un’applicazione della decomposizione serie parallelo di una
funzione di trasferimento. Per trattare tale applicazione ci rifacciamo al testo di un’ esercitazione
di Stima e Filtraggio svolta lo scorso anno, i cui dati preliminari sono qui riportati:
Il processo x(·), stazionario del secondo ordine, é descritto dal modello ARMA
(1 − 0.2z −1 + 0.9z −2 )x(t) = (1 − 1.4z −1 + 1.2z −2 )n(t) t ∈ Z
con n(·) un rumore bianco di varianza σ 2 . Del processo si considera la versione rumorosa
z(t) = x(t) + w(t)
t∈Z
dove w(·) é un rumore colorato, scorrelato da x(·), di densitá spettrale
Sw (ejθ ) =
λ2
2 − cos θ
θ ∈ [−π, π].
Infine, si ottiene il processo delle osservazioni y(·) tramite un trasduttore di trasferenza
F (z) =
1 − 2z −1
1 − 0.8z −1
|z| > 0.8
pilotato da z(·).
In particolare consideriamo la funzione di trasferimento di un filtro interpolatore.
La stima lineare a minima varianza d’errore x̂(t|Z) rappresenta la stima di x(t) considerando
le misure { y(s), s ∈ Z}. Utilizzando la teoria di Wiener-Kolmogorov si ha che la funzione di
trasferimento del filtro che produce tale stima é data dalla seguente espressione:
H∞ (z) =
Sxy (z)
Sy (z)
con Sxy (z) spettro incrociato dei processi x(t) e y(t) e con Sy (z) spettro di y(t).
Con opportuno codice Matlab (che qui non viene riportato, dato che va oltre l’obiettivo di
questo elaborato) si ottiene
z 7 − 7.367z 6 + 20.19z 5 − 31.75z 4 + 30.77z 3 − 18.51z 2 + 6.253z − 0.8
z 7 − 10.07z 6 + 31.7z 5 − 54.02z 4 + 61.33z 3 − 39.2z 2 + 17.13z − 2
A questo punto per implementare il filtro interpolatore é possibile seguire due metodi:
H∞ (z) =
1
• Metodo delle “due passate”
• Mayne-Fraser
La realizzazione secondo il metodo delle “due-passate” si basa sulla possibilitá di rappresentare
la funzione di trasferimento H∞ (z) del filtro interpolatore mediante il prodotto:
H∞ (z) = H− (z)H+ (z)
con H+ (z) che rappresenta la funzione di trasferimento causale e H− (z) che rappresenta la
parte anticausale.
Il metodo introdotto da Mayne e Fraser consiste nell’implementazione del filtro interpolatore di trasferenza razionale H∞ (z) mediante l’uso in parallelo di due filtri, entrambi stabili,
l’uno causale di trasferenza K+ (z) e l’altro anticausale di trasferenza K− (z), corrispondenti alla
scomposizione additiva:
H∞ (z) = K− (z) + K+ (z)
Per la separazione di H∞ (z) in H− (z) e in H+ (z) secondo le specifiche dettate nel documento
é stata generata una funzione chiamata two steps.m (si veda sezione 3). Eseguendo tale funzione
si ottiene:
H+ (z) =
z 4 − 2.235z 3 + 2.294z 2 − 1.14z + 0.1786
z 4 − 9.219z 3 + 23.19z 2 − 27.82z + 20.65
H− (z) =
1
z 3 − 5.132z 2 + 6.425z − 4.478
z 3 − 0.8475z 2 + 0.6993z − 0.09687
Approccio Stima e Filtraggio
Per calcolare K+ (z) e K− (z) é stata quindi generato il codice riportato in MF.m (sezione 3),
ottenendo:
K+ (z) =
K− (z) =
2
0.1654z 3 − 0.2109z 2 + 0.1116z − 0.02273
z 3 − 0.8475z 2 + 0.6993z − 0.09687
0.8346z 4 − 4.924z 3 + 9.539z 2 − 9.595z + 3.415
z 4 − 9.219z 3 + 23.19z 2 − 27.82z + 20.65
Approccio Tecniche Avanzate di Controllo
Il nostro scopo é quello di verificare, secondo quanto visto a lezione, che si possono ottenere le
funzioni di trasferimento K+ (z) e K− (z) anche attraverso le realizzazioni di stato. Per avere una
notazione simile a quella delle dispense, poniamo:
W1 (z) = H− (z)
W2 (z) = H+ (z)
A questo punto, é stato generato il file TAC.m (sezione 3), ottenendo:
Z1 (z) =
0.8355z 4 − 4.932z 3 + 9.564z 2 − 9.621z + 3.426
z 4 − 9.222z 3 + 23.18z 2 − 27.82z + 20.64
2
0.1645z 3 − 0.2099z 2 + 0.1111z − 0.0227
z 3 − 0.8478z 2 + 0.6994z 1 − 0.0969
dove Z1 (z) e Z2 (z) sono state cosı́ nominate in analogia alle dispense, e chiaramente corrispondono rispettivamente alle K− (z) e K+ (z) .
Z2 (z) =
Le figure 1 e 2 dimostrano come i due tipi di realizzazione serie (i.e. “due passate”) e
parallelo (i.e. Mayne-Fraser) siano entrambe in grado di ricostruire, con precisione paragonabile,
il processo x(·).
15
x(t)
x(t|Z) 2S
10
5
0
−5
−10
−15
100
110
120
130
140
150
Figure 1: Confronto tra x(t) e interpolatore two step
15
x(t)
x(t|Z) MF
10
5
0
−5
−10
−15
100
105
110
115
120
125
130
135
140
145
150
Figure 2: Confronto tra x(t) e interpolatore Mayne Fraser
3
3
Codice sorgente
two steps.m
function [H_pos H_neg]=two_steps(H)
[H_num H_den]=tfdata(H,’v’);
zeri=roots(H_num);
poli=roots(H_den);
poli_stab=[];
poli_in=[];
%separo poli stabili e poli instabili
for i=1:length(poli)
if(abs(poli(i))<1)
poli_stab=[poli_stab poli(i)];
else
poli_in=[poli_in poli(i)];
end
end
%costruisco i numeratori delle funzioni H+ e HH_pos_num=[];
H_neg_num=[];
%se grado d+ maggione grado n allora n+=n
if(length(poli_stab)>=length(zeri))
H_pos_num=zeri;
% altrimenti individuo zeri nell’origine e li metto in n+ finche grado
% relativo H+ =0 se grado relativo >0 completo n+ con altri zeri di n
else
ind=find(zeri);
n_0=length(zeri)-length(ind);
if(length(ind)==length(zeri))
for i=1:length(poli_stab)
H_pos_num=[H_pos_num zeri(i)];
end
for i=length(poli_stab)+1:length(zeri)
H_neg_num=[H_neg_num zeri(i)];
end
else
for i=1:(length(zeri)-length(poli_stab)) %sistemo nun neg
if(length(ind)>=1)
H_neg_num=[H_neg_num zeri(ind(1))];
ind=ind(2:length(ind));
end
end
for i=1:n_0
H_pos_num=[H_pos_num 0];
end
for i=1:length(ind)
H_pos_num=[H_pos_num zeri(ind(1))];
end
end
end
H_pos=tf(poly(H_pos_num),poly(poli_stab),-1);
H_neg=tf(poly(H_neg_num),poly(poli_in),-1);
end
4
MF.m
%funzione che calcola la parte causale a anticausale di una funzione
function [ca cab]=MF(H)
H=minreal(H);
[num den]=tfdata(H,’v’);
[r p k]=residue(num,den);
rc=[];
pc=[];
%Calcolo parte causale
for i=1:length(r)
if abs(p(i))<=1
rc=[rc r(i)];
pc=[pc p(i)];
end
end
rac=[];
pac=[];
%Calcolo parte anticausale
for i=1:length(r)
if abs(p(i))>1
rac=[rac r(i)];
pac=[pac p(i)];
end
end
[nac dac]=residue(rac,pac,[]);
ac=tf(nac,dac,-1);
%Calcolo f(0)
f_0=evalfr(ac,0);
if length(k)==0
f_0=f_0;
else
f_0=f_0+k(length(k));
end
[n d]=residue(rc,pc,f_0);
ca=tf(n,d,-1);
[n d]=residue(rc,pc,f_0/2);
cab=tf(n,d,-1);
end
TAC.m
H_inf_num=[1,-7.367,20.19,-31.75,30.77,-18.51,6.253,-0.8];
H_inf_den=[1,-10.07,31.7,-54.02,61.33,-39.2,17.13,-2];
H_inf=tf(H_inf_num,H_inf_den,-1);
[p,r]=two_steps(H_inf);
W_2=ss(p) %corrisponde a H+
W_1=ss(r) %corrisponde a H%%costruisco la realizzazione di stato di W
A=[W_1.a W_1.b*W_2.c
zeros(3,4) W_2.a]
5
B=[W_1.b+W_2.d
W_2.b]
C=[W_1.c W_1.d*W_2.c]
D=[W_1.d*W_2.d]
%trovo la matrice P adatta che esiste dato che sigma(A_1) non ha elementi
%in comuni con sigma(A_2)
P=lyap(W_1.a,-W_2.a,W_1.b*W_2.c)
T=[eye(4)
P
zeros(3,4)
eye(3)]
%mostro che effettivamente la matrice T corretta
T^(-1)*A*T
%costruisco le realizzazioni di stato di Z_1 e Z_2
G_1=W_1.b*W_2.d-P*W_2.b
G_2=W_2.b
H_1=W_1.c
H_2=W_1.c*P+W_1.d*W_2.c
%gradi di libert
J_1=0.8355
J_2=W_1.d*W_2.d-J_1
Z_1=ss(W_1.a,G_1,H_1,J_1)
%corrisponde a KZ_2=ss(W_2.a,G_2,H_2,J_2)
%corrisponde a K+
%calcolo fdt di Z_1
[num1,den1]=ss2tf(W_1.a,G_1,H_1,J_1)
tf1=tf(num1,den1)
%calcolo fdt di Z_2
[num2,den2]=ss2tf(W_2.a,G_2,H_2,J_2)
tf2=tf(num2,den2)
%calcolo fdt del parallelo
tft=tf1+tf2
6