% Skript analog_filter_3.m, in dem exemplarisch analoge
% Tiefpassfilter mit MATLAB-Funktionen entwickelt werden.
% Es wird das Bessel- und Chebyshev-Filter untersucht.
 
clear;
% ------- Bessel Analogfilter über Prototyptiefpassfilter
nord = 6;              % Ordnung des Filters muss < 25
% Prototypfilter
[z,p,k] = besselap(nord); % Pol- und Nullstellen bzw. Gain k
[b,a] = zp2tf(z,p,k);    % Koeffizienten der Übertragungsfunktion

figure(1),    clf;
freqs(b,a); 
title('Frequenzgang des Bessel-Tiefpassprototyps')

% ------- Umwandlung des Prototyps
fp = 600;       % Durchlassfrequenz 600 Hz
[bt1,at1] = lp2lp(b,a,2*pi*fp);
c1= round(log10(fp/100));   c2 = round(log10(fp*10));
f = logspace(c1,c2,500);

figure(2);      clf;
H1 = freqs(bt1,at1,2*pi*f);
subplot(311), semilogx(f, 20*log10(abs(H1)));
title(['Amplitudengang des Bessel-Filters für fp = ',num2str(fp), 'Hz'])
xlabel('Hz');    ylabel('dB (20 log(abs(H)))');
grid on;     hold on;     
La = axis;    axis([La(1:3), 10]);
La = axis;
plot([fp, fp], [La(3), La(4)],'r');
subplot(312), semilogx(f, unwrap(angle(H1))*180/pi);
title(['Phasengang für fp = ',num2str(fp), 'Hz'])
xlabel('Hz');    ylabel('Grad');
grid on;
subplot(313), semilogx(f, 20*log10(abs(H1)));
title(['Ausschnitt vom Amplitudengang']);
xlabel('Hz');    ylabel('dB (20 log(abs(H)))');
grid on;     hold on;     La = axis;
axis([La(1:2), -20, 5]);  hold on;
plot([fp, fp], [La(3), La(4)],'r');

% ------- Chebyshev Analogfilter über Prototyptiefpassfilter
nord = 6;              % Ordnung des Filters muss < 25
well = 3;              % Welligkeit im Durchlassbereich in dB
% Prototypfilter
[z,p,k] = cheb1ap(nord, well); % Pol- und Nullstellen bzw. Gain k
[b,a] = zp2tf(z,p,k);    % Koeffizienten der Übertragungsfunktion

figure(3),    clf;
freqs(b,a); 
title(['Frequenzgang des Chebyshev-Tiefpassprototyps mit Welligkeit = ',...
    num2str(well),' dB'])

% ------- Umwandlung des Prototyps
fp = 600;       % Durchlassfrequenz 600 Hz
[bt2,at2] = lp2lp(b,a,2*pi*fp);
c1= round(log10(fp/100));   c2 = round(log10(fp*10));
f = logspace(c1,c2,500);

figure(4);      clf;
H2 = freqs(bt2,at2,2*pi*f);
subplot(311), semilogx(f, 20*log10(abs(H2)));
title(['Amplitudengang des Chebyshevfilters für fp = ',num2str(fp),...
    ' Hz und Welligkeit = ',num2str(well),' dB']);
xlabel('Hz');    ylabel('dB (20 log(abs(H)))');
grid on;     hold on;     
La = axis;    axis([La(1:3), 10]);
La = axis;
plot([fp, fp], [La(3), La(4)],'r');
subplot(312), semilogx(f, unwrap(angle(H2))*180/pi);
title(['Phasengang für fp = ',num2str(fp), 'Hz'])
xlabel('Hz');    ylabel('Grad');
grid on;
subplot(313), semilogx(f, 20*log10(abs(H2)));
title(['Ausschnitt vom Amplitudengang']);
xlabel('Hz');    ylabel('dB (20 log(abs(H)))');
grid on;     hold on;     
La = axis;    axis([La(1:3), 10]);
La = axis;
axis([La(1:2), -5, 2]);  hold on;
plot([fp, fp], [La(3), La(4)],'r');

% ------- Antwort auf einem Puls
dt = (1/fp)/10;        % Zeitschrittweite 
t = 0:dt:0.1;      nt = length(t);
u = zeros(1,nt);       % Puls initialisieren
tau = 0.05;            % Dauer des Pulses
u = 1*(t>=0.02 & t<=tau+0.02); % Puls beginnend bei 0.02
my_system1 = tf(bt1, at1);  % System für Bessel-Filter
my_system2 = tf(bt2, at2);  % System für Chebyshev-Filter
y1 = lsim(my_system1,u,t);  % Antwort des Bessel-Filters
y2 = lsim(my_system2,u,t);  % Antwort des Chebyshev-Filters

figure(5);     clf;
subplot(211); plot(t,u);
title('Anregungspuls');
xlabel('s');     grid on;
La = axis;     axis([La(1:2),La(3), 1.2]);
subplot(212);    plot(t,y1,t,y2);
title(['Antwort des Bessel- und Chebyshev-Filters',...
    ' mit Welligkeit = ',num2str(well),' dB']);
xlabel('s');     grid on;

% ------- Impulsantworten
Timp = 0.02;     ti = 0:0.0001:Timp;
h1 = impulse(my_system1, ti);
h2 = impulse(my_system2, ti);
figure(6);    clf;
plot(ti, h1, ti, h2);
title('Impulsantwort des Bessel- und Chebyshev Filters');
xlabel('s');    grid on;

