% Skript anti_alias_2.m, in dem ein Antialiasingfilter
% Typ Butterworth untersucht wird. 
% Arbeitet mit Modell anti_alias3.slx

clear;
% -------- Parameter des Modells
fs = 1000;      % Abtastfrequenz
f1 = 50;        % Signale die nicht zu Aliasing führen
ampl1 = 1;      % Amplitude
f2 = 150;  
ampl2 = 1/3;
f3 = 250;       
ampl3 = 1/5;
f4 = 1200;        % Signal das zu Aliasing führt
%f4 = 1850;       % Signal das zu Aliasing führt
%f4 = 1150;       % Signal das zu Aliasing führt
%f4 = 2150;       % Signal das zu Aliasing führt
ampl4 = 1;
% ------- Parameter des Butterworth_Filters
fp = fs/2;     % Beginn des Sperrbereichs 
nord = 6;       % Ordnung des Filters
% ------- Frequenzgang des Chebyshev Typ 2 Filters
[z,p,k] = buttap(nord);
[b,a] = zp2tf(z,p,k);
% Umwandlung in einem gewünschten Tiefpassfilter
[bt, at] = lp2lp(b,a, 2*pi*fp);
c1= round(log10(fs/100));   c2 = round(log10(fs*10));
f = logspace(c1,c2,500);
H = freqs(bt,at,2*pi*f);

figure(1);      clf;
subplot(311), semilogx(f, 20*log10(abs(H)));
title(['Amplitudengang des Butterworth_Filters für fp = ',...
    num2str(fp), 'Hz'])
xlabel('Hz');    ylabel('dB (20 log(abs(H)))');
grid on;     hold on;     
La = axis;    axis([La(1:2),-100, 10]);

% ------- Aufruf der Simulation
Tsim = 5;
nfft = 2048;     
w = sum(hann(nfft).^2)/nfft;  % Mittelere Leistung Fensterfunktion
option = simset('Solver','ode45','MaxStep',1e-4);
sim('anti_alias3',[0,Tsim]);
t  = y.time;
yd = y.data(:,1);    % Zeitdiskretisierter Ausgang
yf = y.data(:,2);    % Ausgang des antialiasing Filters
yt = y.data(:,3);    % Summe aller Eingänge
yi = y.data(:,4);    % Summe der Nutz-Eingänge

subplot(312), plot(t, yt);
La = axis;    axis([t(1000), t(2000), La(3:4)]);
title('Summe aller Eingangssignale');
xlabel('s');  grid on; 
subplot(313), plot(t,yf, t, yd, 'r');
 axis tight;  La = axis;    axis([t(1000), t(2000), La(3:4)]);
title('Ausgang des Antialiasingfilters und des A/D-Wandlers');
xlabel('s');  grid on;   

figure(2);    clf;
plot(t,yi,t,yf);
La = axis;    axis([t(1000), t(2000), La(3:4)]);
title('Das Nutzsignal und das Signal am Ausgang des Filters');
xlabel('Zeit in s');     grid on;


% -------- Leistungsspektrum
Y = Y.signals.values(:,1,end);
figure(3);    clf;
plot((0:nfft-1)*fs/nfft, 10*log10(Y));
title('Leistungsspektrum in dBW');
xlabel('Hz');    grid on;

% -------- Parseval Theorem
Pgesamt = sum(Y),      % Gesamte Leistung aus dem Spektrum
nt = length(t);
Pzeit = sum(yd.^2)/nt, % Gesamte Leistung aus dem Zeitsignal
