% Skript modal_analyse_1.m, in dem die Übertragungsfunktion eines
% mechanischen Systems mit einer Pulsanregung identifiziert wird
% Arbeitet mit Modell modal_analyse1.slx

clear;
% -------- Parameter des Feder-Masse-System
m2 = 5;      k2 = 10;           c2 = 0.3;
% Optimale Tilgerparameter nach Hartog
m1 = 2;      mu = m1/m2;        k1 = mu*k2/((1+mu)^2);
w1 = sqrt(k1/m1);               c1 = 2*w1*m1*sqrt(3*mu/(8*(1+mu)^3));
% Ohne Tilger
m1 = 2;      k1 = 0;            c1 = 0; 
% -------- Zustandsmodell
A = [-c1/m1, -k1/m1, c1/m1, k1/m1;
    1 0 0 0;
    c1/m2, k1/m2, -(c1+c2)/m2 -(k1+k2)/m2;
    0 0 1 0];
B = [0, 0, 1/m2, 0]';    C = eye(4,4);    D = zeros(4,1);
[eig_vekt, lambda] = eig(A);    % Eigenvektoren und Eigenwerte
lambda = diag(lambda);
f1 = imag(lambda(1))/(2*pi),
f2 = imag(lambda(3))/(2*pi),
fmax = max([f1, f2]);    Tmin = 1/fmax;

% -------- Aufruf der Simulation
dt = 0.5;           % Schrittweite der Simulation
Tsim = 200;         % Dauer der Simulation
t = 0:dt:Tsim-dt;   nt = length(t);
% Anregung mit Gausspuls
tau = 4;            % Verspätung für den Gausspuls
Tgauss = 0.1;       % Parameter für die Ausdehnung des Pulses
ug_puls = exp(-(t-tau).^2/Tgauss);
selekt1 = [0,0,0,1];      % Selektion der Lage der Masse 2
selekt2 = [0,0,1,0];      % der Geschwindigkeit der Masse 2
selekt3 = [0,1,0,0];      % der Lage der Masse 1
selekt4 = [1,0,0,0];      % der Geschwindigkeit der Masse 1

options = simset('Solver','ode4','FixedStep', dt);
sim('modal_analyse1', [0,Tsim-dt], options);
y2 = simout.data(:,1);    %  Lage der Masse 2
v2 = simout.data(:,2);    %  Geschwindigkeit der Masse 2
h2 = y2/(sum(ug_puls)*dt);% Geschätzte Impulsantwort   
nh = length(h2);    
H2 = dt*fft(h2).*exp(j*2*pi*(0:nh-1)*tau/(dt*nh)).'; 
 % Übertragungsfunktion aus der Impulsantwort mit Phasenkorrektur       
H21 = fft(y2)./(fft(ug_puls).'+eps); % Übertragungsfunktion Y(jw)/X(jw)

figure(1);   clf;
subplot(211), plot(t, ug_puls);
title('Gauss-Anregungspuls (Hammerstoss)');    xlabel('s');    grid on;
subplot(212), plot(t, y2);
title('Antwort auf Gauss-Anregungspuls (Lage Hauptmasse)');    xlabel('s'); 
grid on;

figure(2);    clf;
subplot(211), plot((0:nh-1)/(dt*nh), 20*log10(abs(H2)),...
    (0:nh-1)/(dt*nh), 20*log10(abs(H21)),'r');
title('|H(j\omega)| aus Impulsantwort und direkt aus Y(j\omega)/X(j\omega)');
xlabel('Hz');     grid on;      ylabel('dB');
subplot(212), plot((0:nh-1)/(dt*nh), angle(H2)*180/pi,...
    (0:nh-1)/(dt*nh), angle(H21)*180/pi,'r');
title('Winkel(H(j\omega)) aus Impulsantwort und direkt aus Y(j\omega)/X(j\omega)');
xlabel('Hz');     grid on;      ylabel('Grad');

figure(3);    clf;
subplot(211), plot((0:nh-1)/(dt*nh), 20*log10(abs(fft(ug_puls))));
title('|DFT| des Anregungspulses (Hammerstoss)');
xlabel('Hz');     grid on;      ylabel('dB');
subplot(212), plot((0:nh-1)/(dt*nh), 20*log10(abs(fft(y2))));
title('|DFT| des Antwort (Lage Hauptmasse)');
xlabel('Hz');     grid on;      ylabel('dB');

