% Skript impuls_antwort_3.m, in dem die Übertragungsfunktion 
% durch einen Stoss identifiziert wird
% Arbeitet mit Modell impuls_antwort3.slx

clear;
% -------- Parameter der Schaltung
R = 50;          L = 0.01;      C = 0.1e-6;
f0 = sqrt(1/(L*C))/(2*pi);   % Resonantfrequenz
Ts = 1e-5;
% -------- Anregung mit Gausspuls
nfft = 512;
Tsim = Ts*nfft;
t = 0:Ts:Tsim-Ts;
ug_puls = exp(-(t-1e-3).^2/1e-10);

% -------- Aufruf der Simulation mit Modell impuls_antwort3.slx
options = simset('Solver','ode4','FixedStep', Ts);
sim('impuls_antwort3', [0, Tsim-Ts], options);
uc = simout.data(:,1);
i =  simout.data(:,2);
t =  simout.time;

figure(1);   clf;
subplot(211), plot(t,ug_puls)
title('Anregung als Gausspuls');
xlabel('s');     grid on;
subplot(212), plot(t, uc);
title('Antwort uc');
xlabel('s');     grid on;

% -------- Komplexer Frequenzgang (Übertragungsfunktion)
% Annäherung der FT der Impulsantwort
Uc = fft(uc);          nh = length(Uc);  % FT Annäherungen
Ug = fft(ug_puls);
Hi = Uc./(Ug+eps).';            % geschätzte Übertragungsfunktion
% Korrekte übertragungsfunktion
b = 1;     a = [L*C, R*C, 1];
[H,w] = freqs(b, a, (0:nfft-1)*2*pi/(nfft*Ts));

figure(2);     clf;
subplot(211), plot(w/(2*pi),abs(H),    (0:nh-1)/(Ts*nh),abs(Hi));
title('Amplitudengang der geschätzten Übertragungsfunktion');
xlabel('Hz');     grid on;     La = axis;
axis([0, 1/Ts, -1, La(4)]);

subplot(212), plot(w/(2*pi),angle(H),  (0:nh-1)/(Ts*nh),angle(Hi));
title('Phasengang der geschätzten Übertragungsfunktion');
xlabel('Hz');     grid on;     La = axis;
axis([0, 1/Ts, La(3:4)]);

figure(3);     clf;
plot((0:nh-1)/(Ts*nh),abs(Uc),    (0:nh-1)/(Ts*nh),abs(Ug));
title(['Betrag der FFT der Antwort ',...
    ' und der FFT des Anregungspulses']);
xlabel('Hz');     grid on;     La = axis;
axis([0, 1/Ts, -1, La(4)]);

% -------- Geschätzter Frequenzgang aus der geschätzte Impulsantwort
hges = uc/(Ts*sum(ug_puls));         % Geschätzte Impulsantwort
delta_t = 1e-3;                      % Usprung des Gausspulses
Hges = Ts*fft(hges).*exp(j*w*delta_t).'; % Geschätzte Übertragungsfunktion
% mit Kompensation der Phasenverschiebung wegen des Puls-Ursprungs 

figure(4);     clf;
subplot(211), plot(w/(2*pi),abs(H),    (0:nh-1)/(Ts*nh),abs(Hges));
title('Amplitudengang aus der geschätzten Impulsantwort ');
xlabel('Hz');     grid on;     La = axis;
axis([0, 1/Ts, -1, La(4)]);

phi_ges = angle(Hges);
subplot(212), plot(w/(2*pi),angle(H),  (0:nh-1)/(Ts*nh),angle(Hges));
title('Phasengang aus der geschätzten Impulsantwort');
xlabel('Hz');     grid on;     La = axis;
axis([0, 1/Ts, La(3:4)]);

