% Skript fourier_transf_puls_3.m, in dem die Fourier-
% Transformation einer abklingenden Exponentialfunktion
% untersucht wird. Arbeitet mit Modell fourier_transf_puls3.slx
clear;
% -------- Parameter des Systems
N = 128;           % Anzahl der Abtastwerte im Untersuchungsintervall
fs = 400;          % Abtastfrequenz
Ts = 1/fs;         % Abtastperiode
Tsig = N*Ts;       % Periode das das aperiodische Signal ergibt
% -------- Eingangssignal (abklingenden Exponentialfunktion)
a = 1/(Tsig/2.5);    % Parameter der Exponentialfunktion
%a = 1/(Tsig/1);    % Parameter der Exponentialfunktion 
dt = Ts/10;        % Schrittweite für die kontinuierliche Funktion
t = 0:dt:Tsig-dt;
yi = exp(-a*t);    % kontinuierliche Exponentialfunktion
% -------- Aufruf der Simulation
Tf = 3*Tsig;       % 3 Perioden als Dauer der Simulation
w = ones(1,N);     % Ohne Fenster-Funktion (rechteckiges Fenster)
%w = exp(-(0:N-1)/50);    % Exponential-Fenster-Funktion
wf = sqrt(N/sum(w.^2));  % Korrekturfaktor wegen der Fenster-Funktion

sim('fourier_transf_puls3', [0,Tf]);
Y = y.signals.values;  % Beträge der FFT*Ts
Y = Y(:,:,end);        % Betrag der FFT*Ts des letzten Datenpuffers
phi = y3.signals.values; % Winkel der FFT des letzten Datenpuffers
phi = phi(:,:,end);

figure(1);   clf;
xdw = y1.signals.values(:,1,end);
xd =  y2.signals.values(:,1,end);
subplot(311), plot(t, yi);
title(['Zeitkontinuierliche Exponentialfunktion mit 1/a = ',...
    num2str(1/a),'s']);
xlabel('s');    grid on;
La = axis;      axis([La(1), max(t), 0, La(4)]);

subplot(312), stem((0:N-1)*Ts, [xd, xdw]);
title(['Zeitdiskrete Sequenz ohne und mit Fenster-Funktion']);
xlabel('s');    grid on;
axis tight;

subplot(313), plot((-(N/2):(N/2-1))*fs/N, 10*log10(fftshift(Y)),'r');
%subplot(212), stem((-(N/2):(N/2-1))*fs/N, 10*log10(fftshift(Y)),'r');
title(['Betrag der FFT*Ts und korrekter Betragsspektrum']);
xlabel('Hz');     grid on;    ylabel('10*log10(Ts*|Xk|)');
axis tight;

% ------- Analytisches Fourier-Spektrum
df = fs/N;
f = -fs/2:df:fs/2-df;
Yid = 1./(sqrt(a^2 + (2*pi*f).^2));   % Korrekter Betragspektrum
phi_id = -atan2(2*pi*f,a);            % Korrekter Phasenspektrum
hold on;
plot((-(N/2):(N/2-1))*fs/N, 10*log10(Yid));
La = axis;     axis([La(1:2),-40, -5]);
hold off;   legend('FT über die DFT', 'Analytische FT');
% ------- Parseval Theorem
E_zeit = sum(xd.^2)*Ts,
E_freq = sum(Y.^2)*fs/N,

figure(2);     clf;
subplot(211), plot((-(N/2):(N/2-1))*fs/N, 10*log10(fftshift(Y)),'r'); 
title('Betragspektrum und Betrag der DFT*Ts');
xlabel('Hz');       ylabel('10*log10(| |)');       grid on;
hold on;
plot((-(N/2):(N/2-1))*fs/N, 10*log10(Yid));
hold off;
subplot(212), plot((-(N/2):(N/2-1))*fs/N, fftshift(phi),'r'); 
title('Phasenspektrum und Winkel der DFT');
xlabel('Hz');       ylabel('Rad)');       grid on;
hold on;
plot((-(N/2):(N/2-1))*fs/N, phi_id);
hold off;