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

clear;
% -------- Parameter der Schaltung
R = 50;          L = 0.01;      C = 0.1e-6;
f0 = sqrt(1/(L*C))/(2*pi);   % Resonantfrequenz
% -------- Anregung mit Gausspuls
nfft = 1024;
Tsim = 0.005;        Ts = Tsim/nfft;
t = 0:Ts:Tsim-Ts;    nt = length(t);
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);
tsim =  simout.time;

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

% -------- Komplexer Frequenzgang (Übertragungsfunktion)
% Annäherung der FT der Impulsantwort
%uc_disk = interp1(tsim,uc,(0:nfft-1)*Ts); % Diskretisierung der Antwort
%ug_disk = interp1(tsim,ug_puls,(0:nfft-1)*Ts,'cubic'); % Diskretisierung der Anregung
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, 15000, -1, 10]);
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, 15000, La(3:4)]);

figure(3);     clf;
plot((0:nh-1)/(Ts*nh),abs(Uc*Ts),    (0:nh-1)/(Ts*nh),abs(Ug*Ts));
title('Amplitudengang der geschätzten Übertragungsfunktion und |FFT| des Anregungspulses');
xlabel('Hz');     grid on;     La = axis;



