% stidemo
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

clear, clf
randn('seed',0)
B = [1 1]; A = [1 1 1]; n = 101;
w = linspace(0,pi,n)';
H = polyval(B,j*w)./polyval(A,j*w);
HH = H + randn(n,1)/7 + randn(n,1)*j/7;
for k = 1:5
 [A(k,:),B(k,:),RAB,Hfit] = sti(w,HH,2,2,k,eps);
 RAB = 2*sqrt(diag(RAB));
 Asd(k,:) = RAB(1:2)'; Bsd(k,:) = RAB(3:4)';
end

subplot(2,2,1)
plot(w,[real(HH) imag(HH)],'-')
hold on
plot(w,[real(H) imag(H)],':')
plot(w,[real(Hfit) imag(Hfit)],'--')
xlabel('f [rad/s]'),title('Transfer function')
set(gca,'Xtick',[0 pi/2 pi],'Xticklabels',['  0 ';'pi/2';' pi '])
axis([0 pi -2 2])
hold off

subplot(2,2,2)
pos = get(gca,'pos'); h = gca;
set(gca,'pos',[pos(1:3) pos(4)/4])
errorbar(1:5,A(:,2),Asd(:,1),Asd(:,1),'o')
text(5.6,.9,'a1')
xlabel('Iteration #')
h(2) = axes('pos',[pos(1) pos(2)+pos(4)/4 pos(3) pos(4)/4])
errorbar(1:5,A(:,3),Asd(:,2),Asd(:,2),'o')
text(5.6,.9,'a0')
h(3) = axes('pos',[pos(1) pos(2)+pos(4)/2 pos(3) pos(4)/4])
errorbar(1:5,B(:,1),Bsd(:,1),Bsd(:,1),'o')
text(5.6,.9,'b1')
h(4) = axes('pos',[pos(1) pos(2)+pos(4)/4*3 pos(3) pos(4)/4])
errorbar(1:5,B(:,2),Bsd(:,2),Bsd(:,2),'o')
text(5.6,.9,'b0')
title('Parameters')
set(h,'Xtick',1:5,'Ytick',[.6 .8 1 1.2],'xlim',[.5 5.5],'ylim',[.5 1.3])
set(h(2:4),'Xticklabels','')

[AA,BB] = meshgrid(0:5,1:6);
J = zeros(size(AA));
for k = 1:prod(size(AA))
 [A,B,RAB,Hfit] = sti(w,HH,AA(k),BB(k),100,1e-6);
 J(k) = sqrt(mean(abs(HH-Hfit).^2));
end
subplot(2,2,3)
set(mesh(AA,BB-1,J),'EdgeColor','r'), view(150,30)
set(gca,'xtick',0:5,'ytick',0:5), grid on
axis([0 5 0 5 0 1]), axis square
xlabel('deg(A)'), ylabel('deg(B)')
title('Loss function')

HH = [H+randn(n,1)/14+randn(n,1)*j/14 H/10+randn(n,1)/14+randn(n,1)*j/14];
[A,B,RAB,Hfit] = sti(w,HH,2,2);
subplot(2,2,4)
pos = get(gca,'pos');
plot(w,[real(HH(:,2)) imag(HH(:,2))],'-r',w,[real(Hfit(:,2)) imag(Hfit(:,2))],'--y')
xlabel('f [rad/s]')
axis([0 pi -.5 .5])
set(gca,'pos',[pos(1) pos(2) pos(3) pos(4)/2],'Xtick',[0 pi/2 pi],'Xticklabels',['  0 ';'pi/2';' pi '],'Ytick',[-.2 0 .2])
axes('pos',[pos(1) pos(2)+pos(4)/2 pos(3) pos(4)/2])
plot(w,[real(HH(:,1)) imag(HH(:,1))],'-r',w,[real(Hfit(:,1)) imag(Hfit(:,1))],'--y')
set(gca,'Xtick',[0 pi/2 pi],'Xticklabels','','ytick',[-1 0 1])
axis([0 pi -2 2])
title('Transfer function')
