% Skript grp_komp_1.m, in dem mit einem Allpass-Filter 
% der nichtlineare Phasengang eines IIR-Filter 
% kompensiert wird

clear;
% -------- IIR-Tiefpassfilter das entzerrt wird
fr = 0.2*2;      % Relative Frequenz MATLAB-konv.
nord = 6;        % Ordnung des Filters
dp = 0.1;        % Welligkeit im Durchlassbereich dB
ds = 60;         % Dämpfung im Sperrbereich dB
hd = fdesign.lowpass('N,Fp,Ap,Ast',nord,fr,dp,ds);
Hd = design(hd,'ellip');
%fvtool(Hd)
% ------- Entzerrungsallpassfilter
f = 0:0.001:fr;    % Frequenzen für den Durchlassbereich
g = grpdelay(Hd,f,2); % Es wird nur der Durchlassbereich entzerrt
Gd = max(g)-g;  % Die Kompensationsgruppenlaufzeit
nord_g = 8;     % Ordnung des Allpass-Kompensationsfilter
[be, ae] = iirgrpdelay(nord_g, f, [0,fr], Gd);
                % Allpass-Kompensationsfilter
% ------- Frequenzgänge
[bd, ad] = sos2tf(Hd.sosMatrix); % Übertragungsfunktion ohne
verst = prod(Hd.ScaleValues);    % die Skalierungsfaktoren
bd = bd*verst;% Verstärkung mit den Faktoren der sos-Abschnitte
[H,w] = freqz(bd,ad);
figure(1);   clf;
subplot(221), plot(w/(2*pi), 20*log10(abs(H)));
title('Amplitudengang des Ellip-Filters');
xlabel('f/fs');    grid on;
subplot(223), plot(w/(2*pi), unwrap(angle(H)));
title('Phasengang des Ellip-Filters');
xlabel('f/fs');    grid on;
[He,w] = freqz(be,ae);
Hk = H.*He;           % Die kompensierte Übertragungsfunktion
subplot(222), plot(w/(2*pi), 20*log10(abs(Hk)));
title('Amplitudengang der Kaskade');
xlabel('f/fs');    grid on;
subplot(224), plot(w/(2*pi), unwrap(angle(Hk)));
title('Phasengang der Kaskade');
xlabel('f/fs');    grid on;
% ------- Einheitspulsantwort der Kaskade
figure(2);    clf;
bgesamt = conv(bd, be);    agesamt = conv(ad, ae);
hgesamt = filter(bgesamt, agesamt,[1, zeros(1,100)]);
stem(0:100, hgesamt);
title(['Einheitspulsantwort der Kaskade: Elliptisches-',...
    'und Kompensationsfilter']);
xlabel('Index n');   grid on;
% ------- Frequenzgang des Allpol-Kompensationsfilter
figure(3);   clf;
subplot(211), plot(w/(2*pi), 20*log10(abs(He)));
title('Amplitudengang des Allpol-Filters');
xlabel('f/fs');    grid on;
La = axis;    axis([La(1:2), -40, 10]);
subplot(212), plot(w/(2*pi), unwrap(angle(He)));
title('Phasengang des Allpol-Filters');
xlabel('f/fs');    grid on;

