% Filter ECG Example

figure(1)
plot(time(1:700),ecg(1:700))
title('Clean ECG')
xlabel('Time (s)')
ylabel('mV')

% 60 Hz contaminated noisy ECG
figure(2)
plot(time(1:700),ecgnse(1:700))
title('Noisy ECG')
xlabel('Time (s)')
ylabel('mV')

% Use all-zero filter
h1 = [0.5 0 0.5];
ecg1 = conv(h1,ecgnse);
figure(3)
plot(time(51:750),ecg1(51:750)) %avoid initial transient
title('Filtered ECG (all-zero)')
xlabel('Time (s)')
ylabel('mV')

% Use pole-zero filter
b = [1 0 1]*1.81/2;
a = [1 0 0.81];
ecg2 = filter(b,a,ecgnse);
figure(4)
plot(time(51:750),ecg2(51:750)) %avoid initial transient
title('Filtered ECG (pole-zero)')
xlabel('Time (s)')
ylabel('mV')

% Evaluate frequency response of each filter
[H1,w] = freqz(h1,1);
[H2,w]= freqz(b,a);
figure(5)
subplot(211)
plot(w,abs(H1))
title('All-zero filter frequency response magnitude')
xlabel('Freq (rads)')
subplot(212)
plot(w,abs(H2))
title('All-zero filter frequency response magnitude')
xlabel('Freq (rads)')

% Evaluate distortion of time domain signal
ecg3 = conv(h1,ecg);
ecg4 = filter(b,a,ecg);
figure(6)
plot(time(300:360),ecg(300:360),'k')
hold on
plot(time(300:360),ecg3(300:360),'r')
plot(time(300:360),ecg4(300:360),'g')
title('Black: Orig; Red: All-zero; Green: Pole-zero')
xlabel('Time (s)')
ylabel('mV')
