% Skript aktiv_tilger_3.m, in dem eine aktive Tilgung
% untersucht wird. Arbeitet mit Modell aktiv_tilger3.mdl
clear;
s = tf('s');
% -------- Parameter des Systems
m1 = 200;     k1 = 100;     c1 = 20;
f1 = sqrt(k1/m1)/(2*pi),
m2 = 1;       k2 = 10;      c2 = 10;
f2 = sqrt(k2/m2)/(2*pi),
% -------- Analyse des Systems
Ai = [m1*s^2+(c1+c2)*s + (k1+k2), -(c2*s+k2);
    -(c2*s+k2), m2*s^2+c2*s+k2];
Bi = [1,0;0,1];
H = inv(Ai)*Bi;         H = minreal(H);
[b11, a11] = tfdata(H(1,1));   [b12, a12] = tfdata(H(1,2));
[b21, a21] = tfdata(H(2,1));   [b22, a22] = tfdata(H(2,2));
b11 = b11{:};    a11 = a11{:};
b12 = b12{:};    a12 = a12{:};
b21 = b21{:};    a21 = a21{:};
b22 = b22{:};    a22 = a22{:};
% ------- PID-Regler für die aktive Tilgung
P = 0;    I = 0;     D = 2000;
Hsys = H(1,1)*feedback(1, (P+I/s+D*s/(s/100+1))*H(1,2)); 
  % Übertragungsfunktion von der Störung Fe1 bis zur Lage y1
Hsys = minreal(Hsys);           
[bsys, asys] = tfdata(Hsys);
bsys = bsys{:},    asys = asys{:},
roots(asys),
% -------- Aufruf der Simulation mit aktive Tilgung
Tfinal = 1000,         dt = 0.01;
Tnoise = 10*dt; 
sim('aktiv_tilger31', [0:dt:Tfinal]);
t = y.time;
y11 = y.signals.values(:,1);
y12 = y.signals.values(:,2);
y1 = y.signals.values(:,3);
nt = length(t);         nd = 1:nt;
fnoise = 1/Tnoise,        % Bandbreite des Rauschsignals
Tdarst = 2000/fnoise;       ndarst = fix(Tdarst/dt);
if ndarst <= nt
     nd = nt - ndarst : nt;
else nd = 1:nt;
end;    
Q = std(y1(nd)),    % Gütefaktor
figure(1);    clf;
subplot(311), plot(t(nd),y11(nd));
title(['Koordinate y11  ( Bandbreite des Rauschsignals fnoise = ',...
    num2str(fnoise),' Hz )']);
xlabel('Zeit in s');    grid on;   axis tight;
subplot(312), plot(t(nd),y11(nd),t(nd),y12(nd));
title('Koordinate y11 und y12');
xlabel('Zeit in s');    grid on;   axis tight;
subplot(313), plot(t(nd),y1(nd));
title('Koordinate y1 = y11+y12');
xlabel('Zeit in s');    grid on;   axis tight;
% -------- Spektralen Leistungsdichten
ndez = 50;          % Dezimierung des Signals y1
Fs = 1/(ndez*dt);
[Pxx,f] = pwelch(y1(1:ndez:end),hamming(512),[],512,Fs);
figure(2);     clf;
plot(f, 10*log10(Pxx));  
title(['Spektrale Leistungsdichte mit und ohne Tilgung; ',...
    'Bandbreite der Rauschkraft Fe1 =', num2str(fnoise),' Hz']);
xlabel('Hz');    ylabel('dB/Hz');     grid on;
hold on;
% -------- Aufruf der Simulation ohne aktive Tilgung
P = 0;    I = 0;    D = 0;
sim('aktiv_tilger31', [0:dt:Tfinal]);
t = y.time;
y11 = y.signals.values(:,1);
y12 = y.signals.values(:,2);
y1 = y.signals.values(:,3);
[Pxx1,f] = pwelch(y1(1:ndez:end),hamming(512),[],512,Fs);
plot(f, 10*log10(Pxx1),'r');                           % Plot the PSD.
hold off;





