% Script H_3dplot_kont_1.m in dem die Funktion H(s) über die 
% komplexe Ebene in 3D dargestellt wird
% Nach G. Doblinger, TU-Wien, 03-2001
clear
% ------- Filterwahl
nord = 7;       % Ordnung der Filter (2*nord für Bandpass)
delta_p = 1;       % Welligkeit im Durchlassbereich in dB
delta_s = 60;      % Dämpfung im Sperrbereich
% Elliptisches-Filter
%[z,p,k] = ellip(nord, delta_p, delta_s, 2*pi*100,'s');          % Tiefpass
%[z,p,k] = ellip(nord, delta_p, delta_s, 2*pi*100,'high','s');   % Hochpass   
%[z,p,k] = ellip(nord, delta_p, delta_s, 2*pi*[100,200],'s');     % Bandpass
% Butterworth-Filter
%[z,p,k] = butter(nord, 2*pi*100,'s');                  % Tiefpass
[z,p,k] = butter(nord, 2*pi*100,'high','s');           % Hochpass
%[z,p,k] = butter(nord, 2*pi*[100,200],'s');             % Bandpass
% Chebyschev-Filter
%[z,p,k] = cheby1(nord, delta_p, 2*pi*100,'s');             % Tiefpass
%[z,p,k] = cheby1(nord, delta_p, 2*pi*100,'high','s');      % Hochpass
%[z,p,k] = cheby1(nord, delta_p, 2*pi*[100,200],'s');        % Bandpass
% -------- Übertragungsfunktion
[b, a] = zp2tf(z,p,k);      % Koeffizienten des Zählers und Nenners
z_f = z;    % Zwischenspeichern der Null- und Polstellen
p_f = p;
% -------- Bereiche für die Darstellung der Übertragungsfunktion
links = min([real(z).', real(p).']);    
rechts = 0;
oben = max([imag(z).', imag(p).']);
unten = -oben;
% Komplexe Variable s, über die, die Übertragungsfunktion dargestellt
% wird
%[x, y] = meshgrid(links*2:(rechts-links)/100:rechts*2,...
%   1.2*unten:(oben-unten)/100:oben*1.2);
[x, y] = meshgrid(links*1.2:(rechts-links)/40:rechts*1.2,...
    3*unten:(oben-unten)/40:oben*3);
z = x + j * y;
% -------- Berechnung der Übertragungsfunktion im gezeigten Bereich
H = polyval(b,z)./polyval(a,z);
H = 20*log10(abs(H));
i1 = find(H > 100);
H(i1) = 100;
jl = find(H < -100);
H(jl) = -100;
close all
% -------- Animation
pos = [0.01 0.2 0.6 0.7];
hfig = figure('name','Z Domain','Units','normal','Position',pos);
%colordef(hfig,'black');
colordef(hfig,'white');
colormap('jet');
surf(x,y,H,'FaceColor','interp','EdgeColor','none','FaceLighting','phong');
hold on;

% -------- Frequenzgang (Schnittstelle mit jw)
yr = 1.2*unten:(oben-unten)/100:oben*1.2;
zr = j*yr;

Hr = (polyval(b,zr)) ./ (polyval(a,zr));
phase = unwrap(angle(Hr));
Hr = 20*log10(abs(Hr));

i1 = find(Hr > 100);    Hr(i1) = 100;
jl = find(Hr < -100);   Hr(jl) = -100;
plot3(zeros(1, length(yr)), yr, Hr);
plot3(zeros(1, length(yr)), yr, -100*ones(length(Hr),1));
hold off;     axis tight
camlight right
view([55,30]);
grid on
xlabel('Re\{z\}');      ylabel('Im\{z\}');
zlabel('|H(z)| in dB');
axis vis3d
for n = 55:15:415
   view([n,30]);
   drawnow
end
axis normal
%#####################################   
figure(2);      clf;
subplot(211), plot(yr/(2*pi), Hr - max(Hr));
title('Amplitudengang');
xlabel('Frequenz Hz');  ylabel('dB');
grid;  axis tight;

subplot(212), plot(yr/(2*pi), phase);
title('Phasengang');
xlabel('Frequenz Hz');   ylabel('Rad');
grid;  axis tight;

% --------- Koeffizienten der Übertragungsfunktion
disp('Koeffizieneten des Zaehlers '), b
disp('Koeffizieneten des Nenners '), a

figure(3);    clf;
%plot(p_f, 'x');    hold on;   plot(z_f + j*eps,'o');
zplane(z_f, p_f);
title('Pol- Nullstellen Platzierung');
xlabel('Realachse');    ylabel('Imaginärachse');
grid on;

% -------- Frequenzgang mit besserer Auflösung
figure(4);    clf;
fmin = 100/10;   fmax = 100*10;
c1 = round(log10(fmin));    c2 = round(log10(fmax));
f = logspace(c1, c2, 2000);
H = freqs(b,a,2*pi*f);

subplot(211), semilogx(f, 20*log10(abs(H)));
title('Amplitudengang des Filters');
xlabel('Frequenz Hz');   grid on;
axis tight;
subplot(212), semilogx(f, angle(H));
title('Phasengang des Filters');
xlabel('Frequenz Hz');   grid on;
axis tight;



