%--------------------------------------------------------------------------
fp=fopen('G:\A1001D1.DAT','r'); % Open file and create file pointer , Fs=250Hz, V=5uV/bit
RecLen=1000;
ECG=fread(fp,RecLen,'int16'); % Read 1000 discrets from the file
T = [ 1 : 1 : RecLen ]; % Time Samples
Fdist = 50; % Input Signal Frequency [Hz]
Fs = 250; % Sampling Frequency [Hz]
Ampl=2000; % Base line amplitude [uV]
Power_Line = sin(2*pi*Fdist/Fs*T)*Ampl/5; % Synthesis of 2mV (in discrets) sine wave - emulating power line interference
ECG_PL=ECG+Power_Line.'; %Add drift to the ECG signal
ECG_PL=ECG_PL*5; % convert to uV
ECG=ECG*5;
% moving average procedure
ECG1=zeros(size(ECG));
ECG1(3:RecLen-3)=(ECG_PL(1:RecLen-5)+ECG_PL(2:RecLen-4)+ECG_PL(3:RecLen-3)+ECG_PL(4:RecLen-2)+ECG_PL(5:RecLen-1))/5;
ERR1=ECG-ECG1;
Time=0.004:0.004:RecLen*0.004; % Time vector
% plot results
figure(1)
subplot(4,1,1), plot(Time, ECG), title('Original ECG'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,2), plot(Time, ECG_PL), title('ECG with power line noise'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,3), plot(Time, ECG1), title('ECG moving average'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,4), plot(Time, ERR1), title('Error'), xlabel('[Sec]'), ylabel('[uV]');
grid on
%--------------------------------------------------------------------------
% 2. Substraction procedure
%--------------------------------------------------------------------------
LinData=zeros(RecLen);
power_buf=zeros(5);
IsLin=false;
M=100; % Linearity criteria
K=0; % Linearity counter
j=3;
for i=3:RecLen-6
D1=ECG_PL(i+5)-ECG_PL(i);
D2=ECG_PL(i+6)-ECG_PL(i+1);
if (abs(D1-D2)<M)
if (K<5)
K=K+1;
end
else
K=0;
end
if (K==5)
LinData(i)=1;
avg=(ECG_PL(i-2)+ECG_PL(i-1)+ECG_PL(i)+ECG_PL(i+1)+ECG_PL(i+2))/5;
power_buf(j)=ECG_PL(i)-avg; % Update power line buffer
end
ECG1(i)=ECG_PL(i)-power_buf(j);
j=j+1;
if (j>5)
j=1;
end
end
ERR1=ECG-ECG1;
figure(2)
subplot(5,1,1), plot(Time, ECG), title('Original ECG'), xlabel('[Sec]'), ylabel('[uV]')
subplot(5,1,2), plot(Time, ECG_PL), title('ECG with base line drift'), xlabel('[Sec]'), ylabel('[uV]')
subplot(5,1,3), plot(Time, LinData), title('ECG Linearity criteria'), xlabel('[Sec]'), ylabel(''), axis([0, RecLen*0.004, 0, 2]);
subplot(5,1,4), plot(Time, ECG1), title('Cleared ECG'), xlabel('[Sec]'), ylabel('[uV]'), axis([0, RecLen*0.004, -5000, 5000]);
subplot(5,1,5), plot(Time, ERR1), title('Error'), xlabel('[Sec]'), ylabel('[uV]'), axis([0, RecLen*0.004, -100, +100]);
grid on

b=[-1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 1-1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 -1/17 0 0 0 0 ]
