%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% EE 675 - Fall 2015
%%%% Simple Code to show How to Open the ECG data files
%%%%
%%%% Your project should try to filter the ECG with artifact
%%%% Target is to make it look like the original
%%%%
%%%% Look at "readme.txt" file for content description of data files
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
load('s118.mat');
load('s119.mat');

range= 1:10000;
fs=  360;       % sampling frequency = 500 Sa/s
%range = range*1/fs;         % convert time axis into seconds


% figure(1)
% sig= s118em12;
% subplot (3,1,1), plot(range*1/fs, s118(1,range)), title('Original ECG Signal - 118')
% subplot (3,1,2), plot(range*1/fs, sig(1,range)), title('ECG Signal with Electode Motion')
% subplot (3,1,3), plot(range*1/fs, -s118(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 

figure(1)
sig= s118bw12;
subplot (3,1,1), plot(range*1/fs, s118(1,range)), title('Original ECG Signal - 118')
subplot (3,1,2), plot(range*1/fs, sig(1,range)), title('ECG Signal with Baseline Wander')
subplot (3,1,3), plot(range*1/fs, -s118(1,range)+sig(1,range)), title('Difference')
xlabel('Time Samples')


% figure(5)
% sig= s118ma12;
% subplot (3,1,1), plot(range*1/fs, s118(1,range)), title('Original ECG Signal - 118')
% subplot (3,1,2), plot(range*1/fs, sig(1,range)), title('ECG Signal with EMG Noise')
% subplot (3,1,3), plot(range*1/fs, -s118(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
% figure(2)
% sig= s119em12;
% subplot (3,1,1), plot(range*1/fs, s119(1,range)), title('Original ECG Signal - 119')
% subplot (3,1,2), plot(range*1/fs, sig(1,range)), title('ECG Signal with Electode Motion')
% subplot (3,1,3), plot(range*1/fs, -s119(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
% figure(4)
% sig= s119bw12;
% subplot (3,1,1), plot(range*1/fs, s119(2,range)), title('Original ECG Signal - 119')
% subplot (3,1,2), plot(range*1/fs, sig(2,range)), title('ECG Signal with Baseline Wander')
% subplot (3,1,3), plot(range*1/fs, -s119(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
% figure(6)
% sig= s119ma12;
% subplot (3,1,1), plot(s119(1,range)), title('Original ECG Signal - 119')
% subplot (3,1,2), plot(sig(1,range)), title('ECG Signal with EMG Noise')
% subplot (3,1,3), plot(-s119(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: you can read Testdata1 and testdata2 the same way
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function Hd = ecg_bw1
%ECG_BW1 Returns a discrete-time filter object.

%
% MATLAB Code
% Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
%
% Generated on: 11-Oct-2016 19:08:10
%

% Equiripple Highpass filter designed using the FIRPM function.

% All frequency values are in Hz.
Fs = 360;  % Sampling Frequency

Fstop = 0.1;             % Stopband Frequency
Fpass = 1;               % Passband Frequency
Dstop = 0.1;           % Stopband Attenuation
Dpass = 0.057501127785;  % Passband Ripple
dens  = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop, Fpass]/(Fs/2), [0 1], [Dstop, Dpass]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

N
% [EOF]
N2= mean(grpdelay(Hd));  
sig= s118bw12(1,:);
sigf= filter(Hd, sig);
sigf= sigf(N2+1:end);

figure
subplot (3,1,1), plot(sig(1,1:10000-N2)), title('ECG Signal with Baseline Wander')
subplot (3,1,2), plot(sigf), title('ECG Signal Filtered Using FIR Filter ')
subplot (3,1,3), plot( -s118bw12(1,1:10000-N2)+sigf), title('Difference')
xlabel('Time Samples')


%function Hd = ecg_bw2
%ECG_BW2 Returns a discrete-time filter object.

%
% MATLAB Code
% Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
%
% Generated on: 11-Oct-2016 20:14:31
%

% Butterworth Highpass filter designed using FDESIGN.HIGHPASS.

% All frequency values are in Hz.
Fs = 360;  % Sampling Frequency

Fstop = 0.5;         % Stopband Frequency
Fpass = 1;           % Passband Frequency
Astop = 60;          % Stopband Attenuation (dB)
Apass = 1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);

% [EOF]

N2= round(mean(grpdelay(Hd)));  
sig= s118bw12(1,:);
sigf= filter(Hd, sig);
sigf= -sigf(N2+1:end);

figure
subplot (3,1,1), plot(sig(1,1:10000-N2)), title('ECG Signal with Baseline Wander')
subplot (3,1,2), plot(sigf), title('ECG Signal Filtered Using IIR Filter')
subplot (3,1,3), plot( -s118bw12(1,1:10000-N2)+sigf), title('Difference')
xlabel('Time Samples')


