% Skript frac_delay_2.m, in dem ein 'Fractional Delay' Filter
% untersucht wird. Arbeitet mit Modell frac_delay2.slx und Solver
% 'FixedStepDiscrete'

clear;
% -------- Parameter der Simulation
% Tiefpassfilter für den Eingangssignal
nord = 64;      fp = 0.02;    % fp <= 0,025
hTP = fir1(nord, fp*2);    % Einheitspulsantwort des FIR-Tiefpass
Ts = 0.1;
N = 20;     % (N+1 Koeffizienten des Frac-Delay-Filters)
n = 0:N;
nk = 0:0.1:N;
D =10.6;   % Sollte kleiner als N/2 sein und nur eine Dezimalstelle
% Einheitspulsantwort des Filters
alpha = 1;
%alpha = 0.85
h = alpha*sinc(alpha*(n-D));
h = h.*hanning(length(h))';
hk = alpha*sinc(alpha*(nk-D));
hk  = hk.*hanning(length(hk))';
figure(1);    clf;%--------------
if N > 20
   subplot(211),  stem(n, h);   hold on;
   plot(nk, hk);   hold off;
   title('Einheitspulsantwort des ''Fractional Delay'' Filters');
   xlabel('n');    grid on;
   subplot(212),  stem(n, h);   hold on;
   plot(nk, hk);   hold off;
   title('Einheitspulsantwort des ''Fractional Delay'' Filters (Ausschnitt)');
   xlabel('n');    grid on; 
   La = axis;    axis([N/2-10, N/2+10, La(3:4)]);
else
   stem(n, h);   hold on;
   plot(nk, hk);   hold off;
   title('Einheitspulsantwort des ''Fractional Delay'' Filters');
   xlabel('n');    grid on; 
end;
% ------- Frequenzgang des Frac-Delay-Filters
nfft = 512;
[H,w] = freqz(h,1, nfft);
[g, w1] = grpdelay(h,nfft);

figure(2);    clf;%--------------
subplot(311),   plot(w/(2*pi), 20*log10(H));
title('Amplitudengang des Fractional-Delay-Filters');
xlabel('f/fs');   grid on;
subplot(312),   plot(w/(2*pi), angle(H));
title('Phasengang des Fractional-Delay-Filters');
xlabel('f/fs');   grid on;
subplot(313),   plot(w1/(2*pi), g);
title('Gruppenlaufzeit als Anzahl Abtastwerte');
xlabel('f/fs');   grid on;
La = axis;   axis([La(1:2), 0, 1.5*D]);
% ------- Aufruf der Simulation
Ts1 = Ts/10;         % Abtastperiode für das Eingangssignal
delay1 = D*10+10;
Tsim = 20;
my_options = simset('Solver','FixedStepDiscrete');
sim('frac_delay2',[0,Tsim], my_options);
y1 = y.Data(:,1);    % Quasikontinuierlicher Eingang
y2 = y.Data(:,2);    % Ausgang des 'Fractional Delay' Filters
t = y.Time;

fehler = y1(1:10:end)-y2(1:10:end);
Q = var(fehler),
figure(3);    clf;%--------------
subplot(211), stairs(t,y1);  hold on;
stairs(t, y2);     grid on;     hold off;
title(['Verspätetes kontinuierliches und mit Fractional Delay Filter',...
    'verspätetes Signal']);
xlabel('s');     grid on;
subplot(212), stairs(t(1:10:end),fehler);
title(['Fehler des Fractional Delay Filters']);
xlabel('s');     grid on;

