% rif. "Progetto in frequenza: esercizio riassuntivo"
% specifiche a), b) e c'')

close all;clear all;

% processo modificato Fhat(s) ottenuto dal processo P(s)
% con l'aggiunta di un polo nell'origine e opportuno guadagno

% la pulsazione di attraversamento desiderata e' w*_t=30
% il margine di fase desiderato e' m*_fi=55

numFhat=[100];
denFhat=conv([1,0],conv([0.1, 1],[0.01, 1]));

% diagrammi di Bode Fhat(s)

figure(1)
bode(numFhat,denFhat)
pause

% calcolo valori attuali pulsazione di attraversamento e margine di fase 

figure(2)
margin(numFhat,denFhat)
pause

% calcolo modulo e fase di Fhat in corrispondenza alla 
% pulsazione di attraversamento desiderata 

w=30; 
[magw,phasew]=bode(numFhat,denFhat,w);
phasew
magwdB=20*log10(magw)
pause

% in corrispondenza a w*_t=30 e' necessario ottenere simultaneamente   
% un anticipo di fase di circa 55-2=53 gradi e un'attenuazione di 18 dB
%
% un'anticipatrice adeguata introdurra' certamente un'amplificazione
% notevole, il che richiedera' l'uso di un'attenuatrice

% primo passo: scelta rete anticipatrice
% in previsione della successiva introduzione della rete 
% attenuatrice, conviene ottenere un anticipo di fase lievemente maggiore
% pari a 53+3=56 gradi

% per ottenere un anticipo di fase di circa 56 gradi aumentando
% il meno possibile la pulsazione di attraversamento si sceglie m_a=11
% con w tau_a =3 (con w=w*_t=30) dai diagrammi universali

tau_a=3/30;
m_a=11;
numR_a=[tau_a, 1];
denR_a=[tau_a/m_a, 1];

% calcolo modulo e fase del processo intermedio Fint(s)= R_a(s) Fhat(s) 
% in corrispondenza alla pulsazione di attraversamento desiderata 

numFint=conv(numFhat,numR_a);
denFint=conv(denFhat,denR_a);

w=30;
[magw,phasew]=bode(numFint,denFint,w);
magwdB=20*log10(magw)
phasew
pause

% secondo passo: scelta rete attenuatrice
% per ottenere una attenuazione di circa magwdB=9.77 dB contenendo entro
% 3 gradi il ritardo di fase si sceglie m_i=3
% e w tau_i=60 (con w=w*_t=30) dai diagrammi universali

tau_i=60/30;
m_i=3;
numR_i=[tau_i/m_i, 1];
denR_i=[tau_i, 1];

% funzione compensatrice complessiva

figure(3)
bode(conv(numR_a,numR_i),conv(denR_a,denR_i))
pause

% verifica soddisfacimento specifiche sui 
% diagrammi di Bode di F(s)=R_i(s)F_int(s)=R_a(s)R_i(s)Fhat(s)=G(s)P(s)

numF=conv(numFint,numR_i);
denF=conv(denFint,denR_i);

w1=logspace(-1,3,1000);
[magFhat,phaseFhat]=bode(numFhat,denFhat,w1);
[magFint,phaseFint]=bode(numFint,denFint,w1);
[magF,phaseF]=bode(numF,denF,w1);
magFhatdB=20*log10(magFhat);
magFintdB=20*log10(magFint);
magFdB=20*log10(magF);
figure(4)
subplot(211)
semilogx(w1,magFhatdB,'r',w1,magFintdB,'g',w1,magFdB,'b','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
title('processo modificato (rosso), intermedio (verde) e compensato (blu)')
subplot(212)
semilogx(w1,phaseFhat,'r',w1,phaseFint,'g',w1,phaseF,'b',w1,-180*(1+0*w1),'w-','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
pause

% zoom diagrammi precedenti nella zona di interesse

w2=logspace(0,2,100);
[magFhat,phaseFhat]=bode(numFhat,denFhat,w2);
[magFint,phaseFint]=bode(numFint,denFint,w2);
[magF,phaseF]=bode(numF,denF,w2);
magFhatdB=20*log10(magFhat);
magFintdB=20*log10(magFint);
magFdB=20*log10(magF);
figure(5)
subplot(211)
semilogx(w2,magFhatdB,'r',w2,magFintdB,'g',w2,magFdB,'b','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Magnitude (dB)')
axis([1 100 -50 40])
title('processo modificato (rosso), intermedio (verde) e compensato (blu)')
subplot(212)
semilogx(w2,phaseFhat,'r',w2,phaseFint,'g',w2,phaseF,'b',w2,-180*(1+0*w2),'w-','LineWidth',2),grid
xlabel('Frequency (rad/sec)')
ylabel('Phase (deg)')
axis([1 100 -270 -90])
pause

% calcolo valori pulsazione di attraversamento e margine di fase 
% dopo la compensazione

figure(6)
margin(numF,denF)
pause

% se le specifiche su margine di fase e pulsazione di attraversamento
% provenivano da specifiche su modulo alla risonanza e banda passante
% ad anello chiuso si possono osservare i diagrammi di Bode
% della funzione di trasferimento ad anello chiuso W(s)=F(s)/1+F(s)...

[numW,denW] = cloop(numF,denF,-1);

w3=logspace(0,2,100);
[magW,phaseW]=bode(numW,denW,w3);
magWdB=20*log10(magW);

figure(7)
semilogx(w3,magWdB,'-',w3,-3+0*w3),grid
xlabel('pulsazione (rad/sec)')
ylabel('modulo (dB)')
title('risposta armonica sistema retroazionato con compensazione')
pause

clsys=tf(numW,denW);

bandwidth(clsys)                    % banda passante finale
20*log10(getPeakGain(clsys))        % modulo alla risonanza finale

% risposta indiciale finale 
t=0:0.01:2;
figure(10)
step(clsys,t)
title('risposta indiciale sistema retroazionato con compensazione')
