% gezeigt wird die vollständige Kompensation eines reelen stabilen Pols 
% durch eine Nullstelle vs. die "Kompensation" durch eine entsprechenden 
% Nullstelle rechts -> Allpass

% 19.3.2012, G.Diehm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;
% Grafik
MarkerProperties.MarkerSize = 15;
MarkerProperties.LineWidth = 2;
%
LineProperties.LineWidth 	= 2;
LineProperties.Color		= 'k';
%
AxisProperties.LineWidth	= 2;
AxisProperties.FontSize		= 16;
%
bo1=bodeoptions;
bo1.XLim = {[0.01 100]};
bo1.YLim = {[-40 40];[-100 100]};
bo1.Grid = 'on';
bo1.PhaseVisible = 'off';
bo2=bodeoptions;
bo2.XLim = {[0.01 100]};
bo2.YLim = {[-40 40];[-100 100]};
bo2.Grid = 'on';
bo2.MagVisible = 'off';
bo2.PhaseVisible = 'on';

% Alle Systeme haben einen Pol bei -1, im letzten System kommt ein Pol bei
% -0.3 hinzu
Pol  = -1;
PolD = -0.3;
% Die Nullstellen liegen bei -1 bzw +1, zunächst aber ohne Nullstelle, im
% letzten System kommt eine Nullstelle bei +0.3 hinzu
NSteA =  [];
NSteB =  -1;
NSteC =   1;
NSteD = 0.3;
% Systeme definieren:
sysA=zpk(NSteA,Pol,1);
sysB=zpk(NSteB,Pol,1);sysB1=zpk(NSteB,[],1);
sysC=zpk(NSteC,Pol,-1);sysC1=zpk(NSteC,[],1);
sysD=zpk(NSteD,[Pol,PolD],-1);
% In den Titelzeilen der Bilder wird jeweils angezeigt, welche Nullstelle das aktuelle System hat:
Titel = ['Pole: s1 = ',num2str(Pol),'; Nullstelle s2 = '];
TitelA = [Titel,'--'];
TitelB = [Titel,num2str(NSteB)];
TitelC = [Titel,num2str(NSteC)];
TitelD = ['Pole: s1 = ',num2str(Pol),', s2 = ',num2str(PolD),'; Nullstelle s3 = ',num2str(NSteD)];
% Frequenzdaten der Systeme erzeugen
w=logspace(-2,2,100);
[mA,pA]=bode(sysA,w);
[mB,pB]=bode(sysB,w);[mB1,pB1]=bode(sysB1,w);
[mC,pC]=bode(sysC,w);[mC1,pC1]=bode(sysC1,w);
[mD,pD]=bode(sysD,w);

% plots erzeugen
fig=figure;
scr=get(0,'Screensize');
set(fig,'position',[scr(3)/4 scr(4)/3 scr(3)/2 scr(4)/3],'tag','srtdemo2');

sub1=subplot(2,2,[1 3]);
pl1=plot(real(Pol),imag(Pol),'r:x',MarkerProperties);
  hold on;
    l1=line([0 0],[-1 1]);
    set(l1,LineProperties);
    l2=line([-2 2],[0 0]);
    set(l2,LineProperties);
    axis([-2 2 -1 1]);
title(TitelA);
hold off;
sub2=subplot(2,2,2);
    m1=semilogx(w,20*log(mA(1,:)),'r');
    set(m1,'LineWidth',2);
    title('Bodediagramm');
    ylabel('Amplitude [dB]');
    axis([0.01 100 -40 40]);
    grid;    
sub3=subplot(2,2,4);
    p1=semilogx(w,pA(1,:),'r');
    set(p1,'LineWidth',2);
    xlabel('Frequenz [rad/s]');
    ylabel('Phase [°]');
    axis([0.01 100 -270 90]);
    grid;
 
pause;
set(fig,'CurrentAxes',sub1);
hold on;pl2=plot(real(NSteB),imag(NSteB),'b:o',MarkerProperties);hold off;
title(TitelB);
set(fig,'CurrentAxes',sub2);
hold on;
m2=semilogx(w,20*log(mB(1,:)),'k');
m21=semilogx(w,20*log(mB1(1,:)),'b');
hold off;
set(m2,'LineWidth',2);set(m21,'LineWidth',2);
set(fig,'CurrentAxes',sub3);
hold on;p2=semilogx(w,pB(1,:),'k');
p21=semilogx(w,pB1(1,:),'b');hold off;
set(p2,'LineWidth',2);set(p21,'LineWidth',2);

pause;
set(fig,'CurrentAxes',sub1);
set(pl2,'Visible','off');
hold on;pl3=plot(real(NSteC),imag(NSteC),'b:o',MarkerProperties);hold off;
title(TitelC);
set(fig,'CurrentAxes',sub2);
m1=semilogx(w,20*log(mA(1,:)),'r');
    set(m1,'LineWidth',2);
    title('Bodediagramm');
    ylabel('Amplitude [dB]');
    axis([0.01 100 -40 40]);
    grid;
hold on;m3=semilogx(w,20*log(mC(1,:)),'k');
m31=semilogx(w,20*log(mC1(1,:)),'b');hold off;
set(m3,'LineWidth',2);set(m31,'LineWidth',2);
set(fig,'CurrentAxes',sub3);
p1=semilogx(w,pA(1,:),'r');
    set(p1,'LineWidth',2);
    xlabel('Frequenz [rad/s]');
    ylabel('Phase [°]');
    axis([0.01 100 -270 90]);
    grid;
hold on;p3=semilogx(w,pC(1,:)-360,'k');
p31=semilogx(w,pC1(1,:)-180,'b--');hold off;
set(p3,'LineWidth',2);set(p31,'LineWidth',2);

pause;
set(fig,'CurrentAxes',sub1);
set(pl3,'Visible','off');
hold on;pl4=plot(real(NSteD),imag(NSteD),'b:o',MarkerProperties);
pl41=plot(real(PolD),imag(PolD),'r:x',MarkerProperties);hold off;
title(TitelD);
set(fig,'CurrentAxes',sub2);
m1=semilogx(w,20*log(mA(1,:)),'r');
    set(m1,'LineWidth',2);
    title('Bodediagramm');
    ylabel('Amplitude [dB]');
    axis([0.01 100 -40 40]);
    grid;
hold on;m4=semilogx(w,20*log(mD(1,:)),'k--');hold off;
set(m4,'LineWidth',2);
set(fig,'CurrentAxes',sub3);
p1=semilogx(w,pA(1,:),'r');
    set(p1,'LineWidth',2);
    xlabel('Frequenz [rad/s]');
    ylabel('Phase [°]');
    axis([0.01 100 -270 90]);
    grid;
hold on;p4=semilogx(w,pD(1,:)-360,'k');hold off;
set(p4,'LineWidth',2);
pause;
close(gcf);