%--------------------------------------------------------------------------
% Baseline drift removing
%--------------------------------------------------------------------------
fp=fopen('G:\A1001D1.DAT','r'); % Open file and create file pointer , Fs=250Hz, V=5uV/bit
RecLen=3000;
ECG=fread(fp,RecLen,'int16'); % Read 3000 discrets from file
T = [ 1 : 1 : RecLen ]; % Time Samples
Fdist = 0.25; % Input Signal Frequency [Hz]
Fs = 250; % Sampling Frequency [Hz]
Ampl=2000; % Base line amplitude [uV]
BLD = sin(2*pi*Fdist/Fs*T)*Ampl/5; % Synthesis of 2mV sine wave - emulating base line drift
ECG_BLD=ECG+BLD.'; %Add drift to the ECG signal
ECG_BLD=ECG_BLD*5; % convert to uV
ECG=ECG*5;
%Calculating filter coefficents for Fc=0.64 Hz
Fc=0.64; % [Hz]
c1 = 1/[1+tan(Fc*pi/Fs)];
c2 = [1-tan(Fc*pi/Fs)]/[1+tan(Fc*pi/Fs)];
b=[c1 -c1]; a=[1 -c2]; % Set filter matrix coeff
ECG1=filter(b,a,ECG_BLD); % Filtering in forward direction
ECGr=fliplr(ECG1); % Reverse ECG
P=filter(b,a,ECGr); % Filtering in opposite direction
ECG2=fliplr(P); % Reverse ECG
ERR1=ECG-ECG1;
ERR2=ECG-ECG2;
Time=0.004:0.004:RecLen*0.004; % Time vector
% plot results
figure(1)
subplot(4,1,1), plot(Time, ECG), xlabel('[Sec]'), ylabel('[uV]')
title('Original ECG')
subplot(4,1,2), plot(Time, ECG_BLD), xlabel('[Sec]'), ylabel('[uV]')
title('ECG with base line drift')
subplot(4,1,3), plot(Time, ECG1) , xlabel('[Sec]'), ylabel('[uV]')
title('ECG first filtration')
subplot(4,1,4), plot(Time, ECG2), xlabel('[Sec]'), ylabel('[uV]')
title('ECG second filtration')
figure(2)
subplot(4,1,1), plot(Time, ECG), title('Original ECG'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,2), plot(Time, ECG_BLD), title('ECG with base line drift'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,3), plot(Time, ECG1), title('ECG first filtration'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,4), plot(Time, ERR1), title('Error first filtration'), xlabel('[Sec]'), ylabel('[uV]'), axis([0, 12, -2000, 2000]);
grid on
figure(3)
subplot(4,1,1), plot(Time, ECG), title('Original ECG'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,2), plot(Time, ECG_BLD), title('ECG with base line drift'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,3), plot(Time, ECG2), title('ECG second filtration'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,4), plot(Time, ERR2), title('Error second filtration'), xlabel('[Sec]'), ylabel('[uV]'), axis([0, 12, -2000, 2000]);
grid on
figure(4)
subplot(4,1,1), plot(Time, ECG), title('Original ECG'), xlabel('[Sec]'), ylabel('[uV]'), axis([4, 8, -5000, 5000]);
subplot(4,1,2), plot(Time, ECG_BLD), title('ECG with base line drift'), xlabel('[Sec]'), ylabel('[uV]'), axis([4, 8, -5000, 5000]);
subplot(4,1,3), plot(Time, ECG1), title('ECG first filtration'), xlabel('[Sec]'), ylabel('[uV]'), axis([4, 8, -5000, 5000]);
subplot(4,1,4), plot(Time, ECG2), title('ECG second filtration'), xlabel('[Sec]'), ylabel('[uV]'), axis([4, 8, -5000, 5000]);
fp=fopen('G:\A1001D1.DAT','r'); % Open file and create file pointer , Fs=250Hz, V=5uV/bit
RecLen=10000;
ECG=fread(fp,RecLen,'int16'); % Read 10000 discrets from the file
ECG1=filter(b,a,ECG); % Filtering in forward direction
ECGr=fliplr(ECG1); % Reverse ECG
P=filter(b,a,ECGr); % Filtering in opposite direction
ECG2=fliplr(P); % Reverse ECG
ECG=ECG*0.005; % Convert to mV
ECG1=ECG1*0.005;
ECG2=ECG2*0.005;
Time=0.004:0.004:RecLen*0.004; % Time vector
% plot results
figure(5)
subplot(3,1,1), plot(Time, ECG), title('Original ECG from AHA DB file A2001D1.DAT'), xlabel('[Sec]'), ylabel('[mV]');
subplot(3,1,2), plot(Time, ECG1), title('ECG first filtration'), xlabel('[Sec]'), ylabel('[mV]');
subplot(3,1,3), plot(Time, ECG2), title('ECG second filtration'), xlabel('[Sec]'), ylabel('[mV]');