% Skript zerophase_1.m in dem eine Nullphase-Filterung
% untersucht wird

clear;
% -------- Bildung eines EKG-Signals
seg = ecg(500);
x = [seg, seg, seg, seg];     nx = length(x);
noise = 0.015;
xn = x + sqrt(noise)*randn(1,nx);
figure(1);    clf;
subplot(211), plot(0:nx-1, x);
title('EKG-Sequenz ohne Rauschen');   grid on;
xlabel('n');
subplot(212), plot(0:nx-1, xn);
title('EKG-Sequenz mit Rauschen');   grid on;
xlabel('n');

% -------- Spektrum des Signals
segmentLength = 400;
nfft = 512;
[pxx1, w] = pwelch(x,segmentLength,[],nfft);
[pxx2, w] = pwelch(xn,segmentLength,[],nfft);
figure(2);    clf;
subplot(211), plot(w/(2*pi), 10*log10(pxx1));
title('Spektrale Leistungsdichte des EKG-Signals ohne Rauschen');
xlabel('f/fs');    grid on;
subplot(212), plot(w/(2*pi), 10*log10(pxx2));
title('Spektrale Leistungsdichte des EKG-Signals mit Rauschen');
xlabel('f/fs');    grid on;

% -------- Kausale Filterung des EKG-Signals
nord = 128;
h = fir1(nord, 0.1*2);
y = filter(h,1,xn);
% -------- Nullphase-Filterung mit filtfilt
yz = filtfilt(h,1,xn);

figure(3);     clf;
subplot(411), plot(0:nx-1, x);
title('EKG-Signal ohne Rauschanteil')
xlabel('n');   grid on;   axis tight;
subplot(412), plot(0:nx-1, xn); 
title('EKG-Signal mit Rauschanteil')   
xlabel('n');   grid on;   axis tight;
subplot(413), plot(0:nx-1, y); 
title('EKG-Signal mit kausaler Filterung')
xlabel('n');   grid on;   axis tight;
subplot(414), plot(0:nx-1, yz);
hold on;
plot(0:nx-1, [y(nord/2+1:end), zeros(1,nord/2)]);
hold off;
title('EKG-Signal mit Nullphase-Filterung')
xlabel('n');   grid on;   axis tight;

figure(4);     clf;
plot(0:nx-1, x);
hold on;
plot(0:nx-1, [y(nord/2+1:end), zeros(1,nord/2)]);
plot(0:nx-1, yz, 'g');
hold off;
La = axis;    axis([1, 400, La(3:4)]);
title('EKG-Signale mit Nullphase-Filterung (Ausschnitt)')
xlabel('n');   grid on;

