% Skript feder_masse_13.m, in dem ein einfaches Feder-Masse-System
% mit Tilger untersucht wird
clear;
s = tf('s');
% ------- Parameter des Systems
m1 = 100;     k1 = 5000;     c1 = 20;
% -------- Klassische Tilgung nach Den Hartog
m2 = 5;                    mu = m2/m1;
k2 = mu*k1/((1+mu)^2);
c2 = 2*m2*sqrt(k2/m2)*sqrt(3*mu/(8*(1+mu)));  % optimale Dämpfung
%c2 = 1;                               % nicht optimale Dämpfung
Ai = [m1*s^2+(c1+c2)*s+(k1+k2), -(c2*s+k2);   % Laplace-
    -(c2*s+k2), m2*s^2+c2*s+k2];              % Transformationsmatrix
Bi = [1;0];
Ht = inv(Ai)*Bi;     Ht = minreal(Ht);  % Übertragungsfunktionen
Ht1 = Ht(1);         Ht2 = Ht(2);   % von der Anregung bis Y1 und Y2
[bt1, at1] = tfdata(Ht1);      [bt2, at2] = tfdata(Ht2);
bt1 = bt1{:};        at1 = at1{:};  % Koeffizienten der 
bt2 = bt2{:};        at2 = at2{:};  % Übertragungsfunktionen
% -------- Zustandsmodell für y1 und y2
A = [0 1 0 0;-(k1+k2)/m1,-(c1+c2)/m1, k2/m1,c2/m1;
    0 0 0 1;k2/m2,c2/m2,-k2/m2,-c2/m2];
B = [0 1/m1 0 0]';
C1 = [1 0 0 0];   D1 = 0;
C2 = [0 0 1 0];   D2 = 0;
% -------- Antwort auf Stoss
Tfinal = 50;               dt = 0.05;
t = 0:dt:Tfinal-dt;        nt = length(t);
%ua = [0:3,3:-1:0];         % Stoss (Trapezförmig) 
%ua = ones(1,5);            % Stoss Rechteckpuls
ua = exp(-((0:6)-3).^2);   % Stoss Gausspuls
u = [ua, zeros(1,nt-length(ua))];  % Anregung
x10 = 0;     v10 = 0;       % Anfangsbedingungen
my_sys_m1 = ss(A,B,C1,D1);
[yt1,t,xt1] = lsim(my_sys_m1,u',t',[x10,v10,0,0]); % Antwort
figure(1);    clf;
nud = 1:fix(nt/20);
subplot(211),  plot(t(nud),u(nud));
title('Anregungspuls');    xlabel('Zeit in s');    grid on;
subplot(212),  plot(t,yt1)
title('Antwort y1');    xlabel('Zeit in s');    grid on;
% ------- Identifikation mit FFT-Spektren
nfft = length(u);
U = fft(u);        Y1 = fft(yt1);
Huy = Y1./(U.');            % Identifizierter Frequenzgang
nd = fix(nfft/4);
figure(2);      clf;
subplot(211), plot((0:nfft-1)/(dt*nfft), 20*log10(abs(U)));
title(' Betrag der FFT des Anregungspulses');
xlabel('Hz');    ylabel('dB');      grid on;
La = axis;       axis([La(1:2), -150, 1.2*max(20*log10(abs(U)))]);
subplot(212), plot((0:nfft-1)/(dt*nfft), 20*log10(abs(Y1)));
title(' Betrag der FFT der Antwort');
xlabel('Hz');    ylabel('dB');      grid on;
La = axis;       axis([La(1:2), -150, 0.8*max(20*log10(abs(Y1)))]);
% ------- Identifizierter und idealer Amplitudengang
f = (0:nd)/(dt*nfft);      w = j*2*pi*f;
Hid = polyval(bt1,w)./polyval(at1,w);
figure(3);      clf;
plot((0:nd)/(dt*nfft), 20*log10(abs(Huy(1:nd+1))));
hold on;
plot((0:nd)/(dt*nfft), 20*log10(abs(Hid(1:nd+1))),'r');
hold off
  title(['Identifizierter und idealer Amplitudengang der Antwort',...
   ' mit Tilger']);
xlabel('Hz');     ylabel('dB');     grid on;



