%----------------------------------------------------------------------------------------------------------
% Power Line interference and base line drift simultaneously removing by Filter DxN
%----------------------------------------------------------------------------------------------------------
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 noise
ECG_PL=ECG+Power_Line.'; %Add drift to the ECG signal
ECG_PL=ECG_PL*5; % convert to uV
ECG=ECG*5;
Fdist = 0.5; % Input Signal Frequency [Hz]
Drift = sin(2*pi*Fdist/Fs*T)*Ampl; % Synthesis of 2mV (in discrets) sine wave - emulating baseline drift
ECG_PL=ECG_PL+Drift.'; %Power line noise + Drift
D=5;
N=23;
Time=0.004:0.004:RecLen*0.004; % Time vector



ECG2=zeros(size(ECG));
for i=1+D*round(N/2):RecLen-D*round(N/2)
m=0;
for j=-round(N/2)+1:round(N/2)-1
m=m+ECG_PL(i+j*D);
end
m=m/N;
ECG2(i)=ECG_PL(i)-m;
end
figure(3)
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 base line drift'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,3), plot(Time, ECG2), title('ECG filtrated'), xlabel('[Sec]'), ylabel('[uV]')
subplot(4,1,4), plot(Time, ECG-ECG2), title('ERROR'), xlabel('[Sec]'), ylabel('[uV]')

grid on

