% Skript all_pass_zerlegung_1.m, in dem eine Zerlegung
% einer Übertragungsfunktion in zwei parallel geschalteten
% Allpas-Filter untersucht wird

clear;
% ------- Entwurf des Filters, das zerlegt wird 
Hf = fdesign.lowpass('N,F3dB,Ap,Ast',5,0.1*2,1,60);
He = design(Hf,'ellip','FilterStructure','cascadeallpass');
[b1,a1] = tf(He.stage(1).stage(1));   % Allpass A0(z)  
[b2,a2] = tf(He.stage(1).stage(2));   % Allpass A1(z)
% ------- Frequenzgänge
[H1,w1] = freqz(b1,a1);      % Frequenzgang Allpass 1
figure(1);   clf;
subplot(321), plot(w1/(2*pi), 20*log10(abs(H1)));
title('Amplitudengang des Allpassfilters 1');
xlabel('f/fs');    ylabel('dB');    grid on;
La = axis;   axis([La(1:2), -10, 5]);
subplot(323), plot(w1/(2*pi), unwrap(angle(H1)));
title('Phasengang des Allpassfilters 1');
xlabel('f/fs');    ylabel('Rad');    grid on;

[H2,w2] = freqz(b2,a2);      % Frequenzgang Allpass 2
subplot(322), plot(w2/(2*pi), 20*log10(abs(H2)));
title('Amplitudengang des Allpassfilters 2');
xlabel('f/fs');    ylabel('dB');    grid on;
La = axis;   axis([La(1:2), -10, 5]);
subplot(324), plot(w2/(2*pi), unwrap(angle(H2)));
title('Phasengang des Allpassfilters 2');
xlabel('f/fs');    ylabel('Rad');    grid on;
subplot(313), plot(w2/(2*pi), unwrap(angle(H1))-unwrap(angle(H2)));
title('Phasendifferenz Allpassfilter 1 minus Allpassfilter 2');
xlabel('f/fs');    ylabel('Rad');    grid on;

% -------- Frequenzgang der parallelen Allpass-Struktur
H = (H1+H2)/2;                   % Komplexer Frequenzgang

figure(2);   clf;
subplot(211), plot(w2/(2*pi), 20*log10(abs(H)));
title('Amplitudengang der parallelen Struktur');
xlabel('f/fs');    ylabel('dB');    grid on;
La = axis;   axis([La(1:3), 10]);
subplot(212), plot(w2/(2*pi), unwrap(angle(H)));
title('Phasengang der parallelen Struktur');
xlabel('f/fs');    ylabel('Rad');    grid on;

[zHe, pHe, kHe] = zpk(He);   % Pol- nullstellen des IIR-Filters
z1 = roots(b1); p1 = roots(a1);% Pol- nullstellen des Allpass-Filters 1
z2 = roots(b2); p2 = roots(a2);% Pol- nullstellen des Allpass-Filters 2

figure(3);    clf;
subplot(211), plot(pHe,'*'); hold on; plot(zHe,'o');
phi = 0:0.01:2*pi;       % Kreis Radius eins
plot(cos(phi), sin(phi)); grid on;
title('Pol- Nullstellen IIR-Filter');   xlabel('Real');
ylabel('Imaginär'); axis square; axis tight; hold off;
subplot(223), plot(p1,'*'); hold on; plot(z1,'o');
phi = 0:0.01:2*pi;       % Kreis Radius eins
plot(cos(phi), sin(phi)); grid on;
title('Pol- Nullstellen Allpass-Filter 1');   xlabel('Real');
ylabel('Imaginär'); axis tight; hold off;
subplot(224), plot(p2,'*'); hold on; plot(z2,'o');
phi = 0:0.01:2*pi;       % Kreis Radius eins
plot(cos(phi), sin(phi)); grid on;
title('Pol- Nullstellen Allpass-Filter 2');   xlabel('Real');
ylabel('Imaginär'); axis tight; hold off;
% ------- Pole- und Nullstellen
betrag_pHe = abs(pHe),    betrag_zHe = abs(zHe),
betrag_p1 = abs(p1),      betrag_z1 = abs(z1),
betrag_p2 = abs(p2),      betrag_z2 = abs(z2),


