#
# Test_komplexe_Stellen_02.m
#
clc;
clear all;
close all;
format long;
# pkg load ltfat;
pkg load control;
fig=0;
#
global Vini Vonn Tdel Tris Tfal Tonn Tper;
#
i=sqrt(-1);
immax=6;
#
feck=1.0e6;
ima=0.5*sqrt(2);
rea=sqrt(1-ima^2);
#
nz(1)=1;
np(1)=2;
pol(1,1)=-1.0e4;
zer(1,1)=-feck;
pol(1,2)=-1.0e8;
#
nz(2)=2;
np(2)=3;
pol(2,1)=-1.0e4;
zer(2,1)=-feck;
zer(2,2)=-7.0e6*(rea+i*ima);
pol(2,2)=-1.0e7*(rea+i*ima);
pol(2,3)=-1.0e8;
#
nz(3)=2;
np(3)=3;
pol(3,1)=-1.0e4;
zer(3,1)=-feck;
zer(3,2)=-7.0e6*(rea+i*ima);
pol(3,2)=-1.0e7;
pol(3,3)=-1.0e8;
#
nz(4)=2;
np(4)=3;
pol(4,1)=-1.0e4;
zer(4,1)=-feck;
zer(4,2)=-7.0e6;
pol(4,2)=-1.0e7*(rea+i*ima);
pol(4,3)=-1.0e8;
#
nz(5)=2;
np(5)=3;
pol(5,1)=-1.0e4;
zer(5,1)=-feck;
zer(5,2)=-7.0e6;
pol(5,2)=-1.0e7;
pol(5,3)=-1.0e8;
#
nz(6)=1;
np(6)=2;
pol(6,1)=-1.0e4;
zer(6,1)=-0.9*feck*(rea+i*ima);
pol(6,2)=-1.0e8;
#
itmax=8001;
freq=logspace(2,10,itmax);
tmax=600e-6;
t2=linspace(0,tmax,itmax)';
#
# PULSE aus LTspice
#
Vini=0;
Vonn=-1.35;
Tdel=3e-6;
Tris=4e-9;
Tfal=85e-6;
Tonn=0;
Tper=tmax;
#
# Vektor des Input-Impuls
#
fin=ffin(t2);
#
# Zusammenstellung der Übertragungsfunktionen aus Polen und Nullstellen
#
xsys1 = zpk([zer(1,1)],[pol(1,1) pol(1,2)],1);
xsys2 = zpk([zer(2,1) zer(2,2)],[pol(2,1) pol(2,2) pol(2,3)],1);
xsys3 = zpk([zer(3,1) zer(3,2)],[pol(3,1) pol(3,2) pol(3,3)],1);
xsys4 = zpk([zer(4,1) zer(4,2)],[pol(4,1) pol(4,2) pol(4,3)],1);
xsys5 = zpk([zer(5,1) zer(5,2)],[pol(5,1) pol(5,2) pol(5,3)],1);
xsys6 = zpk([zer(6,1)],[pol(6,1) pol(6,2)],1);
[mag(:,1), pha(:,1)] = bode (xsys1,freq);
[mag(:,2), pha(:,2)] = bode (xsys2,freq);
[mag(:,3), pha(:,3)] = bode (xsys3,freq);
[mag(:,4), pha(:,4)] = bode (xsys4,freq);
[mag(:,5), pha(:,5)] = bode (xsys5,freq);
[mag(:,6), pha(:,6)] = bode (xsys6,freq);
#
for m=1:immax
  nmax(m)=1;
  absmag(m)=max(mag(:,m));
  for n=1:itmax
    if absmag(m)==mag(n,m);
      nmax(m)=n;
    endif
  endfor
  fakt(m)=1/mag(nmax(m),m);
endfor
#
xsys1 = zpk([zer(1,1)],[pol(1,1) pol(1,2)],fakt(1));
xsys2 = zpk([zer(2,1) zer(2,2)],[pol(2,1) pol(2,2) pol(2,3)],fakt(2));
xsys3 = zpk([zer(3,1) zer(3,2)],[pol(3,1) pol(3,2) pol(3,3)],fakt(3));
xsys4 = zpk([zer(4,1) zer(4,2)],[pol(4,1) pol(4,2) pol(4,3)],fakt(4));
xsys5 = zpk([zer(5,1) zer(5,2)],[pol(5,1) pol(5,2) pol(5,3)],fakt(5));
xsys6 = zpk([zer(6,1)],[pol(6,1) pol(6,2)],fakt(6));
[mag(:,1), pha(:,1)] = bode (xsys1,freq);
[mag(:,2), pha(:,2)] = bode (xsys2,freq);
[mag(:,3), pha(:,3)] = bode (xsys3,freq);
[mag(:,4), pha(:,4)] = bode (xsys4,freq);
[mag(:,5), pha(:,5)] = bode (xsys5,freq);
[mag(:,6), pha(:,6)] = bode (xsys6,freq);
#
for m=1:immax
  difphi=0.0;
  for n=1:itmax-1
    if (pha(n+1,m)-pha(n,m)) > 160
      difphi=180.0;
    endif
    if (pha(n+1,m)-pha(n,m)) > 330
      difphi=360.0;
    endif
    if (pha(n+1,m)-pha(n,m)) > 690
      difphi=720.0;
    endif
    if difphi > 0.0
      pha(n+1,m)=pha(n+1,m)-difphi;
    endif
  endfor
endfor
#
printf("--- Pole und Nullstellen ermitteln.\n");
for m=1:immax
  phi=pha(:,m)/90.0-0.5;
  nz2(m)=0;
  np2(m)=0;
  for o=fix(min(phi)):fix(max(phi))
    printf("m=%d  o=%d  max(phi)=%g  min(phi)=%g\n",m,o, ...
               max(pha(:,m)),min(pha(:,m)));
    for n=1:itmax-1
      if phi(n) < o && o < phi(n+1)
        nz2(m)=nz2(m)+1;
        zer2(m,nz2(m))=freq(n)+(o-phi(n))*(freq(n+1)-freq(n))/(phi(n+1)-phi(n));
       printf("Nullstelle %d  n=%d  freq(n)=%g interpol. Frequenzwert=%g\n", ...
                                               nz2(m),n,freq(n),zer2(m,nz2(m)));
      endif
      if phi(n) > o && o > phi(n+1)
        np2(m)=np2(m)+1;
        pol2(m,np2(m))=freq(n)+(o-phi(n))*(freq(n+1)-freq(n))/(phi(n+1)-phi(n));
       printf("Polstelle  %d  n=%d  freq(n)=%g interpol. Frequenzwert=%g\n", ...
                                               np2(m),n,freq(n),pol2(m,np2(m)));
      endif
    endfor
  endfor
endfor
printf("--- Pole und Nullstellen sind ermittelt.\n");
#
# Kurvendisskussion
#
for m=1:immax
ma1(:,m)=diffqt(mag(:,m),freq)';
ma2(:,m)=diffqt(freq.*ma1(:,m),freq)';
ph1(:,m)=diffqt(pha(:,m),freq)';
ph2(:,m)=diffqt(freq.*ph1(:,m),freq)';
#
#fig=figure(fig+1);
#f1=1/max(abs(mag(:,m)));
#f2=1/max(abs(freq.*ma1(:,m)'));
#f3=1/max(abs(freq.*ma2(:,m)'));
#semilogx(freq,f1*mag(:,m),"linewidth",1,"r", ...
#         freq,f2*freq.*ma1(:,m)',"linewidth",1,"linestyle","-.","b", ...
#         freq,f3*freq.*ma2(:,m)',"linewidth",1,"linestyle",":","m")
#title('Magnitude');
#xlabel('Frequenz [Hz]');
#ylabel('auf max. Amplitude normiert');
#legend("Magnitude normiert","1. Ableitung normiert","2. Ableitung normiert");
#grid
#
#fig=figure(fig+1);
#f2=1/max(abs(freq.*ph1(:,m)'));
#f3=1/max(abs(freq.*ph2(:,m)'));
#semilogx(freq,pha(:,m)/45,"linewidth",1,"r", ...
#         freq,f2*freq.*ph1(:,m)',"linewidth",1,"linestyle","-.","b", ...
#         freq,f3*freq.*ph2(:,m)',"linewidth",1,"linestyle",":","m")
#title('Phase');
#xlabel('Frequenz [Hz]');
#ylabel('Phase [Grad/45]');
#legend("Phase","1. Ableitung normiert","2. Ableitung normiert");
#grid
#
# Nullstellensuche in 1. und 2. Ableitung
#
printf("--- Nullstellen der 1. und 2. Ableitungen ermittelt\n");
#
na1(m)=0;
na2(m)=0;
np1(m)=0;
np2(m)=0;
for n=1:itmax-1
  if ma1(n,m)*ma1(n+1,m) <= 0.0
    na1(m)=na1(m)+1;
    sta1(na1(m),m)=freq(n)-freq(n).*ma1(n,m)*(freq(n+1)-freq(n))/ ...
                                     (freq(n+1).*ma1(n+1,m)-freq(n).*ma1(n,m));
   printf("Nullstelle ma1 %d  n=%d  freq(n)=%g interpol. Frequenzwert=%g\n", ...
                                          na1(m),n,freq(n),sta1(na1(m),m));
  endif
  if ma2(n,m)*ma2(n+1,m) <= 0.0
    na2(m)=na2(m)+1;
    sta2(na2(m),m)=freq(n)-freq(n).*ma2(n,m)*(freq(n+1)-freq(n))/ ...
                                     (freq(n+1).*ma2(n+1,m)-freq(n).*ma2(n,m));
   printf("Nullstelle ma2 %d  n=%d  freq(n)=%g interpol. Frequenzwert=%g\n", ...
                                          na2(m),n,freq(n),sta2(na2(m),m));
  endif
  if ph1(n,m)*ph1(n+1,m) <= 0.0
    np1(m)=np1(m)+1;
    stp1(np1(m),m)=freq(n)-freq(n).*ph1(n,m)*(freq(n+1)-freq(n))/ ...
                                     (freq(n+1).*ph1(n+1,m)-freq(n).*ph1(n,m));
   printf("Nullstelle ph1 %d  n=%d  freq(n)=%g interpol. Frequenzwert=%g\n", ...
                                          np1(m),n,freq(n),stp1(np1(m),m));
  endif
  if ph2(n,m)*ph2(n+1,m) <= 0.0
    np2(m)=np2(m)+1;
    stp2(np2(m),m)=freq(n)-freq(n).*ph2(n,m)*(freq(n+1)-freq(n))/ ...
                                     (freq(n+1).*ph2(n+1,m)-freq(n).*ph2(n,m));
   printf("Nullstelle ph2 %d  n=%d  freq(n)=%g interpol. Frequenzwert=%g\n", ...
                                          np2(m),n,freq(n),stp2(np2(m),m));
  endif
endfor
endfor
printf("--- Nullstellen der 1. und 2. Ableitungen wurden ermittelt\n");
#
fig=figure(fig+1);
semilogx(freq,freq.*ma1',"linewidth",1)
title('1. Ableitung der Magnituden');
xlabel('Frequenz [Hz]');
ylabel('Steigung der Amplitude');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#
fig=figure(fig+1);
semilogx(freq,freq.*ma2',"linewidth",1)
title('2. Ableitung der Magnituden');
xlabel('Frequenz [Hz]');
ylabel('Krümmung der Amplitude');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#
fig=figure(fig+1);
semilogx(freq,freq.*ph1',"linewidth",1)
title('1. Ableitung der Phase');
xlabel('Frequenz [Hz]');
ylabel('Steigung der Phase');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#
fig=figure(fig+1);
semilogx(freq,freq.*ph2',"linewidth",1)
title('2. Ableitung der Phase');
xlabel('Frequenz [Hz]');
ylabel('Krümmung der Phase');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#
# Plot der Frequenzbereichsdaten
#
fig=figure(fig+1);
semilogx(freq,20*log10(mag),"linewidth",1)
title('Magnitude');
xlabel('Frequenz [Hz]');
ylabel('Verstärkung [dB]');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#
fig=figure(fig+1);
semilogx(freq,pha/45,"linewidth",1)
title('Phase');
xlabel('Frequenz [Hz]');
ylabel('Phase [Deg/45°]');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#
fig=figure(fig+1);
bode(xsys1,xsys2,xsys3,xsys4,xsys5,xsys6);
#
# Berechnung der Zeitbereichsdaten aus den Übertragungsfunktionen
#
fn(:,1) = lsim(xsys1,fin,t2(itmax));
fn(:,2) = lsim(xsys2,fin,t2(itmax));
fn(:,3) = lsim(xsys3,fin,t2(itmax));
fn(:,4) = lsim(xsys4,fin,t2(itmax));
fn(:,5) = lsim(xsys5,fin,t2(itmax));
fn(:,6) = lsim(xsys6,fin,t2(itmax));
#
fig=figure(fig+1);
plot(10^6*t2,fn',"linewidth",1)
title('Systemoutput');
xlabel('Zeit [μs]');
ylabel('Amplitude [V]');
legend('xsys1','xsys2','xsys3','xsys4','xsys5','xsys6');
grid
#