% Skript feder_masse_1.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;
% ------- Übertragungsfunktion
H = 1/(m1*s^2 + c1*s + k1);
b = 1;                      a = [m1, c1, k1];
w01 = sqrt(k1/m1);          f01 = w01/(2*pi), % Eigenfrequenz
fmin = f01/10;              fmax = f01*10;
e1 = round(log10(fmin));    e2 = round(log10(fmax));
f = logspace(e1,e2,5000);   w = 2*pi*f;
Hf = freqs(b,a,w);          % Frequenzgang ohne Tilger
figure(1);     clf;
subplot(211),  semilogx(f, 20*log10(abs(Hf)));
title('Amplitudengang für die Hauptmasse ohne Tilger');  
xlabel('Hz');   grid on;   ylabel('dB');
subplot(212),  semilogx(f, angle(Hf)*180/pi);
title('Phasengang für die Hauptmasse ohne Tilger');      
xlabel('Hz');   grid on;   ylabel('Grad');
% -------- 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 Xs und X1
[bt1, at1] = tfdata(Ht1);      [bt2, at2] = tfdata(Ht2);
bt1 = bt1{:};        at1 = at1{:};  % Koeffizienten der 
bt2 = bt2{:};        at2 = at2{:};  % Übertragungsfunktionen
Hf1 = freqs(bt1,at1,w);     Hf2 = freqs(bt2,at2,w); % Frequenzgänge
figure(2);    clf;
semilogx(f, 20*log10(abs(Hf)), f, 20*log10(abs(Hf1)),...
    f, 20*log10(abs(Hf2)));
title(['Amplitudengang für Y1/Fex ohne und mit Tilger ',...
      ' (nach Den Hartog) bzw. Amplitudengang für Y2/Fex']);
xlabel('Hz');     grid on;          ylabel('dB');

% -------- 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;
% -------- Homogene Lösung ohne und mit Tilger
Tfinal = 50;               dt = 0.01;
t = 0:dt:Tfinal-dt;        nt = length(t);
% Ohne Tilger
x10 = 0.1;     v10 = 0;       % Anfangsbedingungen
A0 = [0,1;-k1/m1,-c1/m1];     % Zustandsmodell ohne Tilger
B0 = [0,1/m1]';    C0 = [1,0];    D0 = 0;
my_sys_o = ss(A0,B0,C0,D0);
figure(3);      clf;
[y1,t,x1] = lsim(my_sys_o,[zeros(nt,1)],t',[x10,v10]);
subplot(311),  plot(t, y1);
title('Homogene Lösung für die Koordinate y1 ohne Tilger')
xlabel('Zeit in s');     grid on;
% Mit Tilger
my_sys_m1 = ss(A,B,C1,D1);
[yt1,t,xt1] = lsim(my_sys_m1,[zeros(nt,1)],t',[x10,v10,0,0]);
subplot(312),  plot(t, yt1);
title('Homogene Lösung für die Koordinate y1 mit Tilger')
xlabel('Zeit in s');     grid on;
my_sys_m2 = ss(A,B,C2,D2);
[yt2,t,xt2] = lsim(my_sys_m2,[zeros(nt,1)],t',[x10,v10,0,0]);
subplot(313),  plot(t, yt2);
title('Homogene Lösung für die Koordinate y2 der Tilgermasse')
xlabel('Zeit in s');     grid on;
% -------- Antwort auf Zufallsanregung
Tfinal = 1000;             dt = 0.05;
t = 0:dt:Tfinal-dt;        nt = length(t);
nst = 1;
randn('seed', 13975);
u = sqrt(nst)*randn(nt,1);             % Anregung
[ys,t,xs] =    lsim(my_sys_o,u,t',zeros(2,1));     % Ohne Tilger
[ys1,t,xts1] = lsim(my_sys_m1,u,t',zeros(4,1));    % Mit Tilger für y1
[ys2,t,xts2] = lsim(my_sys_m2,u,t',zeros(4,1));    % Mit Tilger für y2

figure(4);    clf;
nd = nt-1000+1:nt;
subplot(411), plot(t(nd), u(nd));
title('Anregungskraft Fex');    xlabel('Zeit in s');    grid on;
subplot(412), plot(t(nd), ys(nd));
title('Koordinate y1 der Hauptmasse ohne Tilger');   
xlabel('Zeit in s');    grid on;
subplot(413), plot(t(nd), ys1(nd));
title('Koordinate y1 der Hauptmasse mit Tilger');   
xlabel('Zeit in s');    grid on;
subplot(414), plot(t(nd), ys2(nd));
title('Koordinate y2 der Tilgermasse');   
xlabel('Zeit in s');    grid on;

% -------- Spektrale Leistungsdichte
nfft = 256;
[Pss,f]   =  pwelch(ys,hamming(nfft),50,nfft,1/dt);
[Pssxt1,f] = pwelch(ys1,hamming(nfft),50,nfft,1/dt);
[Pssxt2,f] = pwelch(ys2,hamming(nfft),50,nfft,1/dt);
figure(5);    clf;
plot(f, 10*log10([Pss, Pssxt1, Pssxt2]));
title(' Spektrale Leistungsdichte in dB/Hz');
xlabel('Hz');    ylabel('10*log10(Pss)');  grid on;
legend('Ohne Tilger','Mit Tilger','Tilgermasse');
% -------- Überprüfung der spektralen Leistungsdichte
% über den Satz von Parseval
Pzt = std(ys)^2,        % Leistung ohne Tilger aus Zeitbereich
Pzf = sum(Pss/(dt*nfft)), % aus der spektralen Leistungsdichte

PztT = std(ys1)^2,      % Leistung mit Tilger aus Zeitbereich
PzfT = sum(Pssxt1/(dt*nfft)), % aus der spektralen Leistungsdichte

PztmT = std(ys2)^2,     % Leistung mit Tilger aus Zeit- und 
PzfmT = sum(Pssxt2/(dt*nfft)),  % Frequenzbereich für Tilgermasse


