%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 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;


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


% figure(3)
% sig= s118bw12;
% subplot (3,1,1), plot(s118(1,range)), title('Original ECG Signal - 118')
% subplot (3,1,2), plot(sig(1,range)), title('ECG Signal with Baseline Wander')
% subplot (3,1,3), plot(-s118(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
% figure(5)
% sig= s118ma12;
% subplot (3,1,1), plot(s118(1,range)), title('Original ECG Signal - 118')
% subplot (3,1,2), plot(sig(1,range)), title('ECG Signal with EMG Noise')
% subplot (3,1,3), plot(-s118(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
% figure(2)
% sig= s119em12;
% 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 Electode Motion')
% subplot (3,1,3), plot(-s119(1,range)+sig(1,range)), title('Difference')
% xlabel('Time Samples')
% 
% 
% figure(4)
% sig= s119bw12;
% subplot (3,1,1), plot(s119(2,range)), title('Original ECG Signal - 119')
% subplot (3,1,2), plot(sig(2,range)), title('ECG Signal with Baseline Wander')
% subplot (3,1,3), plot(-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 = bpf3
%BPF3 Returns a discrete-time filter object.

%
% MATLAB Code
% Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
%
% Generated on: 18-Oct-2016 19:47:27
%

% FIR Window Bandpass filter designed using the FIR1 function.

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

Fstop1 = 0.5;             % First Stopband Frequency
Fpass1 = 1;               % First Passband Frequency
Fpass2 = 30;              % Second Passband Frequency
Fstop2 = 50;              % Second Stopband Frequency
Dstop1 = 0.1;             % First Stopband Attenuation
Dpass  = 0.057501127785;  % Passband Ripple
Dstop2 = 0.001;           % Second Stopband Attenuation
flag   = 'scale';         % Sampling Flag

% Calculate the order from the parameters using KAISERORD.
[N,Wn,BETA,TYPE] = kaiserord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 ...
                             1 0], [Dstop1 Dpass Dstop2]);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
Hd = dfilt.dffir(b);

% [EOF]

sig= s118ma12(1,:);
sigf= filter(Hd, sig);
%sigout= filtfilt(Hd.sosMatrix,Hd.ScaleValues,sig);
sigout= filtfilt(Hd.Numerator,1.0,sig);

figure
subplot (4,1,1), plot(sig), title('Original ECG Signal - 118')
subplot (4,1,2), plot(sigf), title('ECG Signal filtered')
subplot (4,1,3), plot(sigout), title('ECG Signal Zero-Phase filtered')
subplot (4,1,4), plot( -sig+sigout), title('Difference')
xlabel('Time Samples')


%-------------------------------------------------------------
%function Hd = bpf4
%BPF4 Returns a discrete-time filter object.

%
% MATLAB Code
% Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
%
% Generated on: 18-Oct-2016 19:54:37
%

% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.

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

Fstop1 = 0.5;         % First Stopband Frequency
Fpass1 = 1;           % First Passband Frequency
Fpass2 = 30;          % Second Passband Frequency
Fstop2 = 50;          % Second Stopband Frequency
Astop1 = 20;          % First Stopband Attenuation (dB)
Apass  = 1;           % Passband Ripple (dB)
Astop2 = 60;          % Second Stopband Attenuation (dB)
match  = 'both';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, ...
                      Astop2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);

% [EOF]

sig= s118ma12(1,:);
sigf= filter(Hd, sig);
sigout= filtfilt(Hd.sosMatrix,Hd.ScaleValues,sig);
%sigout= filtfilt(Hd.Numerator,1.0,sig);

figure
subplot (4,1,1), plot(sig), title('Original ECG Signal - 118')
subplot (4,1,2), plot(sigf), title('ECG Signal filtered')
subplot (4,1,3), plot(sigout), title('ECG Signal Zero-Phase filtered')
subplot (4,1,4), plot( -sig+sigout), title('Difference')
xlabel('Time Samples')



