% Skript rad_aufhaeng_4.m, in dem die Radaufhängung
% eines Fahrzeugs simuliert wird, die mit Zufallssignale 
% angeregt wird.
% Arbeitet mit Simulink-Modell rad_aufhaeng4.slx
clear;
s = tf('s');
% ------- Parameter des Systems
ms = 300;   c = 12000;     k = 50000;
ma = 30;    ct = 10000;    kt = 100000;

f10 = sqrt(k/ms)/(2*pi),   % Eigenfrequenzen ohne Dämpfung
f20 = sqrt(kt/ma)/(2*pi),
% ------- Matrizen des Systems
M = [ms 0;0 ma];   C = [c -c;-c ct+c];   K = [k -k; -k kt+k];
% ------- Poly-Eigenwerte
[X,e] = polyeig(K, C, M);
[e, k] = sort(e),
X = X(:,k);
f1 = abs(imag(e(1)))/(2*pi),    f2 = abs(imag(e(3)))/(2*pi),

% ------- Partikuläre Lösung über Übertragungsfunktion
D = [0;ct];     B = [0;kt];
H = inv(s^2*M + s*C + K)*[s*D + B],
H1 = H(1),      H2 = H(2), % Transfer-Function Daten

[b, a] = tfdata(H);
b1 = b{1};    a1 = a{1};   % Koeffizienten der Übertragungs-
                           % funktion von y zu x1 
b2 = b{2};    a2 = a{2};   % Koeffizienten der Übertragungs-
                           % funktion von y zu x2
S0 = 2e-8;                 % Spektrale Leistungsdichte für die 
             % Anregung in m^2/Hz bei einer Geschwindigkeit von 30 m/s
% -------- Aufruf der Simulation                           
dt = 0.01;    Tfinal = 100;
varianz = S0/dt;   % Varianz der Anregung für der Random-Number Block
nfft = 256;                % Buffer-Größe
wind = sum(hann(nfft));    % Hanning-Fenster Werte
sim('rad_aufhaeng4',[0:dt:Tfinal]);                           

t = x.time;
x1 = x.signals.values(:,1);
x2 = x.signals.values(:,2);

figure(1);    clf;
nd = 1:500;
subplot(211), plot(t(nd), x1(nd));
title('Lage x1 der Masse ms');
xlabel('Zeit in s');     grid on;
subplot(212), plot(t(nd), x2(nd));
title('Lage x2 der Masse ma');
xlabel('Zeit in s');     grid on;

% ------- Spektrale Leistungsdichten
X1fft = Xfft(:,1,end);   % Die letzten gemittelten spektralen
X2fft = Xfft(:,2,end);   % Leistungsdichten

figure(2);    clf;
subplot(211), plot((-nfft/2:nfft/2-1)/(dt*nfft),...
    10*log10(fftshift(X1fft)));
title('Spektrale Leistungsdichte für x1');    grid on;
xlabel('Hz');    ylabel('dBW/Hz')
subplot(212), plot((-nfft/2:nfft/2-1)/(dt*nfft),...
    10*log10(fftshift(X2fft)));
title('Spektrale Leistungsdichte für x2');    grid on;
xlabel('Hz');    ylabel('dBW/Hz')

% ------- Parseval Theorem
P1t = std(x1)^2,     % Leistungen über Zeit
P2t = std(x2)^2,
P1f = sum(X1fft)/(dt*nfft),   % Leistungen über spektrale
P2f = sum(X2fft)/(dt*nfft),   % Leistungsdichte
                           
                           
                           
                           