Matlab code to detect and mark PQRST on your ECG

Homework 6

ECE/BIOM 537: Biomedical Signal Processing

Colorado State University

Student: Minh Anh Nguyen

1.      Algorithms to mark up PQRST challenge

-          the sample frequency is given (250 Hz).

-          Load ECG signal (load signals.mat;)

-          Use EXCEL (copy data from signals.mat to EXCEL data sheet) to display of ECG signal curves

-          Plot the raw ECG signal

-          Remove low frequency component.

-          Thresholding to find Peaks of Interest

-          Mark PQRST and find heart  rate

 

1.       Write algorithm to detect and mark PQRST on your ECG given to you.

 

 

3.     Determine heart rate

Heart rate is 84 beat per minute

 

4      Plot the ½ rate, ¼ rate sampling of the ECG study the impact on the algorithm perform.

             5.      Up sample the ECG to double rate and plot the same perform study.

Plot the  ¼ rate sampling of the ECG 

 

plot the up-sampling  frequency

%Homework 6

%ECE/BIOM 537: Biomedical Signal Processing

%Colorado State University

%Student: Minh Anh Nguyen

%Email: minhanhnguyen@q.com

close all;

clear all;

clc;

%%Select a filename in .mat format and load the file.

%[fname path]=uigetfile('*.mat');

%fname=strcat(path,fname);

%y1 = load(fname );

%file =load('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.mat')

load('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.mat')

disp('Contents of workspace after loading file:')

whos

fs = 250; % find the sampling rate or frequency

y1=xlsread('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.xls');

T = 1/fs;% sampling rate or frequency

% find the length of the data per second

N = length(y1);

ls = size(y1);

t = (0 : N-1) / fs;% sampling period

%t = (0 : N-1) *T;

%t = (0:1:length(y1)-1)/fs;

%subplot (2,2,2)

%plot (t,data);

figure; %subplot(1,2,1);

      plot(t,y1);

    %plot(x,y2, 'g');  

title ('plot of the original of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

y1_n=(y1-min(y1))/(max(y1)-min(y1));         % normalize between 0-1

fnyquist = fs/2;

%% find P

m1=max(y1)*.60;

P=find(y1>=m1);

y1_1500 = y1(1:1850);

t2 = 1:length(y1_1500);

figure;

plot(t2,y1_1500);

title ('plot of subset of the ECG signal')

xlabel ('time (msec)')

ylabel ('Amplitute (mv)')

grid on

%% used the snip code from this website.

%%%%http://www.mathworks.com/help/signal/examples/peak-analysis.html

%Detrending Data

%The above signal shows a baseline shift and therefore does not represent the true amplitude. In order to remove the trend, fit a low order polynomial to the signal and use the polynomial to detrend it.

[p,s,mu] = polyfit((1:numel(y1_1500))',y1_1500,6);

f_y = polyval(p,(1:numel(y1_1500))',[],mu);

ECG_data = y1_1500 - f_y;        % Detrend data

N1= length (y1_1500);

t1 = (0 : N1-1) / fs;% sampling period

figure

%plot(t1,ECG_data); grid on

plot(t2,ECG_data); grid on

ax = axis; axis([ax(1:2) -2.2 2.2])

%ax = axis; axis([ax(1:2) -3.2 3.2])

title('Detrended ECG Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

legend('Detrended ECG Signal')

%Thresholding to Find Peaks of Interest

%The QRS-complex consists of three major components: Q-wave, R-wave, S-wave. The R-waves can be detected by thresholding peaks above 0.5mV. Notice that the R-waves are separated by more than 200 samples. Use this information to remove unwanted peaks by specifying a 'MinPeakDistance'.

[~,locs_Rwave] = findpeaks(ECG_data,'MinPeakHeight',0.5,...

                                    'MinPeakDistance',120);

%Finding Local Minima in Signal

%Local minima can be detected by finding peaks on an inverted version of the original signal.

ECG_inverted = -ECG_data;

[~,locs_Swave] = findpeaks(ECG_inverted,'MinPeakHeight',0.4,...

                                        'MinPeakDistance',120);                                                           

%The following plot shows the R-waves and S-waves detected in the signal.

figure

hold on

plot(t2,ECG_data);

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

%axis([0 1850 -1.1 1.1]); grid on;

axis([0 1850 -2.2 2.2]); grid on;

legend('ECG Signal','R-waves','S-waves');

xlabel('time msec'); ylabel('Voltage(mV)')

title('R-wave and S-wave in ECG Signal')

[~,locs_Twave] = findpeaks(ECG_data,'MinPeakHeight',-0.02,...

                                      'MinPeakDistance',50);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','y');                               

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','T-wave','R-wave','S-wave');

[~,locs_Pwave] = findpeaks(ECG_data,'MinPeakHeight',-0.09,...

                                      'MinPeakDistance',25);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Pwave,ECG_data(locs_Pwave),'x','MarkerFaceColor','y');

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','g');                               

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','P-wave','T-wave','R-wave','S-wave');

[~,locs_qwave] = findpeaks(ECG_data,'MinPeakHeight',-0.2);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_qwave,ECG_data(locs_qwave),'x','MarkerFaceColor','y');

% link and zoom in to show the changes

%linkaxes(ax(1:2),'xy');

%axis(ax,[60 230 0.006 -0.04])

%Next, we try and determine the locations of the Q-waves. Thresholding the peaks to locate the Q-waves results in detection of unwanted peaks as the Q-waves are buried in noise. We filter the signal first and then find the peaks. Savitzky-Golay filtering is used to remove noise in the signal.

smoothECG = sgolayfilt(ECG_data,1,3);

figure

plot(t2,ECG_data,'b',t2,smoothECG,'r'); grid on

axis tight;

xlabel('time msec'); ylabel('Voltage(mV)');

legend('ECG Signal','Filtered Signal')

title('Filtering Noisy ECG Signal')

%We perform peak detection on the smooth signal and use logical indexing to find the locations of the Q-waves.

%[~,min_locs] = findpeaks(-smoothECG,'MinPeakDistance',29);

%[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',2);%Twave

[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',50);

% Peaks between -0.2mV and -0.5mV

%locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 &

%-smoothECG(min_locs)<-0.1); %Twave

locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 & -smoothECG(min_locs)<-0.11);

figure

hold on

plot(t2,smoothECG);

plot(locs_Qwave,smoothECG(locs_Qwave),'rs','MarkerFaceColor','g');

plot(locs_Rwave,smoothECG(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,smoothECG(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('Smooth ECG signal','T-interval','R-wave','S-wave');

%The above figure shows that the QRS-complex successfully detected in the noisy ECG signal.

%Error Between Noisy and Smooth Signal

%Notice the average difference between the QRS-complex in the raw and the detrended filtered signal.

% Values of the Extrema

[val_Qwave, val_Rwave, val_Swave] = deal(smoothECG(locs_Qwave), smoothECG(locs_Rwave), smoothECG(locs_Swave));

meanError_Qwave = mean((y1_1500(locs_Qwave) - val_Qwave))

meanError_Rwave = mean((y1_1500(locs_Rwave) - val_Rwave))

meanError_Swave = mean((y1_1500(locs_Swave) - val_Swave))

%% find PP interval

 i = 0;  %% to make the code start from 0.

 rr = 0; %% each time the code run, rr distance two peaks

 hold off % for the next graph

 rrinterval = zeros(3600,1); % create an array to strore 2 peaks

beat_count =0;

for k = 2 : length(y1)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y1(k)> y1(k-1) && y1(k) > y1(k+1) && y1(k)> 1);

     beat_count = beat_count +1;

     if beat_count ==1;

        rr =0;

     else

         rr = k-i;

         rrinterval(k)=rr;

         i=k;  

     end

    else

        rrinterval(k)= rr;

    end       

end

figure;

plot (rrinterval);

xlabel('Time in sec*10^-2'), ylabel('Distance betweeen 2 Heatbeats (R-R) in sec*10^-2'), title('R-R intervals');

%% find PP interval

%% heart rate analysis

% count the dominat peak

beat_count =0;

for k = 2 : length(y1)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y1(k)> y1(k-1) && y1(k) > y1(k+1) && y1(k)> 1)

         beat_count = beat_count +1;

    end       

end

display (k);

disp('dominant peaks');

%% divide the peak count by the duration in minute

duration_in_sec = N/fs;

duration_in_minute = duration_in_sec/60;

BPM = beat_count/duration_in_minute;

%%%  DFT to describe the signal in the frequency

NFFT = 2 ^ nextpow2(N);

Y = fft(y1, NFFT) / N;

f = (fs / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz

amp = ( 2 * abs(Y(1: NFFT / 2+1))); % Vector containing corresponding amplitudes

figure;

plot (f, amp);

title ('plot single-sided amplitude spectrume of the ECG signal')

xlabel ('frequency (Hz)')

ylabel ('|y(f)|')

grid on;

max_value=max(y1);

mean_value=mean(y1);

threshold=(max_value-mean_value)/2;

%% downsampling ½ sample frequency

close all;

clear all;

clc;

%%Select a filename in .mat format and load the file.

%[fname path]=uigetfile('*.mat');

%fname=strcat(path,fname);

%y1 = load(fname );

%file =load('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.mat')

load('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.mat')

disp('Contents of workspace after loading file:')

whos

fs = 250; % find the sampling rate or frequency

fs2 = 250*1/2;

y1=xlsread('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.xls');

T = 1/fs;% sampling rate or frequency

% find the length of the data per second

N = length(y1);

ls = size(y1);

t = (0 : N-1) / fs;% sampling period

%t = (0 : N-1) *T;

%t = (0:1:length(y1)-1)/fs;

%subplot (2,2,2)

%plot (t,data);

figure; %subplot(1,2,1);

      plot(t,y1);

    %plot(x,y2, 'g');  

title ('plot of the original of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

%%%%%%%%%%%%%

% down sampling 1/2 of frequency sample

y2 = resample(y1,fs2,fs);

N2 = length(y2);

ls2 = size(y2);

t22 = (0 : N2-1) / fs2;% sampling period

figure; %subplot(1,2,1);

      plot(t22,y2);

title ('plot of the down sampling 1/2 frequency sample of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

%y1_n=(y1-min(y1))/(max(y1)-min(y1));         % normalize between 0-1

fnyquist = fs2/2;

%% find P

m1=max(y2)*.60;

P=find(y2>=m1);

y1_1500 = y2(1:1850);

t2 = 1:length(y1_1500);

figure;

plot(t2,y1_1500);

title ('plot of subset of down sampling 1/2 frequency sample the ECG signal')

xlabel ('time (msec)')

ylabel ('Amplitute (mv)')

grid on

%% used the snip code from this website.

%%%%http://www.mathworks.com/help/signal/examples/peak-analysis.html

%Detrending Data

%The above signal shows a baseline shift and therefore does not represent the true amplitude. In order to remove the trend, fit a low order polynomial to the signal and use the polynomial to detrend it.

[p,s,mu] = polyfit((1:numel(y1_1500))',y1_1500,6);

f_y = polyval(p,(1:numel(y1_1500))',[],mu);

ECG_data = y1_1500 - f_y;        % Detrend data

N1= length (y1_1500);

t1 = (0 : N1-1) / fs2;% sampling period

figure

%plot(t1,ECG_data); grid on

plot(t2,ECG_data); grid on

ax = axis; axis([ax(1:2) -2.2 2.2])

%ax = axis; axis([ax(1:2) -3.2 3.2])

title('Detrended down sampling 1/2 frequency sample ECG Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

legend('Detrended ECG Signal')

%Thresholding to Find Peaks of Interest

%The QRS-complex consists of three major components: Q-wave, R-wave, S-wave. The R-waves can be detected by thresholding peaks above 0.5mV. Notice that the R-waves are separated by more than 200 samples. Use this information to remove unwanted peaks by specifying a 'MinPeakDistance'.

[~,locs_Rwave] = findpeaks(ECG_data,'MinPeakHeight',0.5,...

                                    'MinPeakDistance',60);

%Finding Local Minima in Signal

%Local minima can be detected by finding peaks on an inverted version of the original signal.

ECG_inverted = -ECG_data;

[~,locs_Swave] = findpeaks(ECG_inverted,'MinPeakHeight',0.4,...

                                        'MinPeakDistance',60);                                                           

%The following plot shows the R-waves and S-waves detected in the signal.

figure

hold on

plot(t2,ECG_data);

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

%axis([0 1850 -1.1 1.1]); grid on;

axis([0 1850 -2.2 2.2]); grid on;

legend('ECG Signal','R-waves','S-waves');

xlabel('time msec'); ylabel('Voltage(mV)')

title('R-wave and S-wave in down sampling 1/2 frequency sample of ECG Signal')

[~,locs_Twave] = findpeaks(ECG_data,'MinPeakHeight',-0.02,...

                                      'MinPeakDistance',25);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','y');                               

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in down sampling 1/2 frequency sample Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','T-wave','R-wave','S-wave');

[~,locs_Pwave] = findpeaks(ECG_data,'MinPeakHeight',-0.09,...

                                      'MinPeakDistance',12);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Pwave,ECG_data(locs_Pwave),'x','MarkerFaceColor','y');

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','g');                                

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in down sampling 1/2 frequency sample Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','P-wave','T-wave','R-wave','S-wave');

[~,locs_qwave] = findpeaks(ECG_data,'MinPeakHeight',-0.2);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_qwave,ECG_data(locs_qwave),'x','MarkerFaceColor','y');

% link and zoom in to show the changes

%linkaxes(ax(1:2),'xy');

%axis(ax,[60 230 0.006 -0.04])

%Next, we try and determine the locations of the Q-waves. Thresholding the peaks to locate the Q-waves results in detection of unwanted peaks as the Q-waves are buried in noise. We filter the signal first and then find the peaks. Savitzky-Golay filtering is used to remove noise in the signal.

smoothECG = sgolayfilt(ECG_data,1,3);

figure

plot(t2,ECG_data,'b',t2,smoothECG,'r'); grid on

axis tight;

xlabel('time msec'); ylabel('Voltage(mV)');

legend('ECG Signal','Filtered Signal')

title('Filtering Noisy of down sampling 1/2 frequency sample ECG Signal')

%We perform peak detection on the smooth signal and use logical indexing to find the locations of the Q-waves.

%[~,min_locs] = findpeaks(-smoothECG,'MinPeakDistance',29);

%[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',2);%Twave

[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',25);

% Peaks between -0.2mV and -0.5mV

%locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 &

%-smoothECG(min_locs)<-0.1); %Twave

locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 & -smoothECG(min_locs)<-0.11);

figure

hold on

plot(t2,smoothECG);

plot(locs_Qwave,smoothECG(locs_Qwave),'rs','MarkerFaceColor','g');

plot(locs_Rwave,smoothECG(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,smoothECG(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks down sampling 1/2 frequency sample in Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('Smooth ECG signal','T-interval','R-wave','S-wave');

%The above figure shows that the QRS-complex successfully detected in the noisy ECG signal.

%Error Between Noisy and Smooth Signal

%Notice the average difference between the QRS-complex in the raw and the detrended filtered signal.

% Values of the Extrema

[val_Qwave, val_Rwave, val_Swave] = deal(smoothECG(locs_Qwave), smoothECG(locs_Rwave), smoothECG(locs_Swave));

meanError_Qwave = mean((y1_1500(locs_Qwave) - val_Qwave))

meanError_Rwave = mean((y1_1500(locs_Rwave) - val_Rwave))

meanError_Swave = mean((y1_1500(locs_Swave) - val_Swave))

%% find PP interval

 i = 0;  %% to make the code start from 0.

 rr = 0; %% each time the code run, rr distance two peaks

 hold off % for the next graph

 rrinterval = zeros(3600,1); % create an array to strore 2 peaks

beat_count =0;

for k = 2 : length(y1)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y1(k)> y1(k-1) && y1(k) > y1(k+1) && y1(k)> 1);

     beat_count = beat_count +1;

     if beat_count ==1;

        rr =0;

     else

         rr = k-i;

         rrinterval(k)=rr;

         i=k;  

     end

    else

        rrinterval(k)= rr;

    end        

end

figure;

plot (rrinterval);

xlabel('Time in sec*10^-2'), ylabel('Distance betweeen 2 Heatbeats (R-R) in sec*10^-2'), title('R-R down sampling 1/2 frequency sample intervals');

%% find PP interval

%% heart rate analysis

% count the dominat peak

beat_count =0;

for k = 2 : length(y2)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y2(k)> y2(k-1) && y2(k) > y2(k+1) && y2(k)> 1)

         beat_count = beat_count +1;

    end        

end

display (k);

disp('dominant peaks');

%% divide the peak count by the duration in minute

duration_in_sec = N/fs2;

duration_in_minute = duration_in_sec/60;

BPM = beat_count/duration_in_minute;

%%%  DFT to describe the signal in the frequency

NFFT = 2 ^ nextpow2(N2);

Y = fft(y2, NFFT) / N2;

f = (fs2 / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz

amp = ( 2 * abs(Y(1: NFFT / 2+1))); % Vector containing corresponding amplitudes

figure;

plot (f, amp);

title ('plot single-sided amplitude spectrume of 1/2 frequency sample ECG signal')

xlabel ('frequency (Hz)')

ylabel ('|y(f)|')

grid on;

max_value=max(y1);

mean_value=mean(y1);

threshold=(max_value-mean_value)/2;

%%Downsampling ¼ sample frequency

close all;

clear all;

clc;

load('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.mat')

disp('Contents of workspace after loading file:')

whos

fs = 250; % find the sampling rate or frequency

fs1 = 250*1/2;

fs2 = 250*1/4;

y1=xlsread('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.xls');

T = 1/fs;% sampling rate or frequency

% find the length of the data per second

N = length(y1);

ls = size(y1);

t = (0 : N-1) / fs;% sampling period

figure; %subplot(1,2,1);

      plot(t,y1);

    %plot(x,y2, 'g');  

title ('plot of the original of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

% down sampling 1/2 of frequency sample

y2a = resample(y1,fs1,fs);

N1 = length(y2a);

ls1 = size(y2a);

t21 = (0 : N1-1) / fs1;% sampling period

figure; %subplot(1,2,1);

      plot(t21,y2a);

title ('plot of the down sampling 1/2 frequency sample of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

%%%%%%%%%%%%%

% down sampling 1/4 of frequency sample

y2 = resample(y1,63,250);

N2 = length(y2);

ls2 = size(y2);

t22 = (0 : N2-1) / fs2;% sampling period

figure; %subplot(1,2,1);

      plot(t22,y2);

title ('plot of the down sampling 1/4 frequency sample of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

%% find P

m1=max(y2)*.60;

P=find(y2>=m1);

y1_1500 = y2(1:1850);

t2 = 1:length(y1_1500);

figure;

plot(t2,y1_1500);

title ('plot of subset of down sampling 1/4 frequency sample the ECG signal')

xlabel ('time (msec)')

ylabel ('Amplitute (mv)')

grid on

%% used the snip code from this website.

%%%%http://www.mathworks.com/help/signal/examples/peak-analysis.html

%Detrending Data

%The above signal shows a baseline shift and therefore does not represent the true amplitude. In order to remove the trend, fit a low order polynomial to the signal and use the polynomial to detrend it.

[p,s,mu] = polyfit((1:numel(y1_1500))',y1_1500,6);

f_y = polyval(p,(1:numel(y1_1500))',[],mu);

ECG_data = y1_1500 - f_y;        % Detrend data

N1= length (y1_1500);

t1 = (0 : N1-1) / fs2;% sampling period

figure

%plot(t1,ECG_data); grid on

plot(t2,ECG_data); grid on

ax = axis; axis([ax(1:2) -2.2 2.2])

%ax = axis; axis([ax(1:2) -3.2 3.2])

title('Detrended down sampling 1/4 frequency sample ECG Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

legend('Detrended ECG Signal')

%Thresholding to Find Peaks of Interest

%The QRS-complex consists of three major components: Q-wave, R-wave, S-wave. The R-waves can be detected by thresholding peaks above 0.5mV. Notice that the R-waves are separated by more than 200 samples. Use this information to remove unwanted peaks by specifying a 'MinPeakDistance'.

[~,locs_Rwave] = findpeaks(ECG_data,'MinPeakHeight',0.5,...

                                    'MinPeakDistance',30);

%Finding Local Minima in Signal

%Local minima can be detected by finding peaks on an inverted version of the original signal.

ECG_inverted = -ECG_data;

[~,locs_Swave] = findpeaks(ECG_inverted,'MinPeakHeight',0.4,...

                                        'MinPeakDistance',30);                                                           

%The following plot shows the R-waves and S-waves detected in the signal.

figure

hold on

plot(t2,ECG_data);

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

%axis([0 1850 -1.1 1.1]); grid on;

axis([0 1850 -2.2 2.2]); grid on;

legend('ECG Signal','R-waves','S-waves');

xlabel('time msec'); ylabel('Voltage(mV)')

title('R-wave and S-wave in down sampling 1/4 frequency sample of ECG Signal')

[~,locs_Twave] = findpeaks(ECG_data,'MinPeakHeight',-0.02,...

                                      'MinPeakDistance',13);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','y');                               

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in down sampling 1/4 frequency sample Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','T-wave','R-wave','S-wave');

[~,locs_Pwave] = findpeaks(ECG_data,'MinPeakHeight',-0.09,...

                                      'MinPeakDistance',6);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Pwave,ECG_data(locs_Pwave),'x','MarkerFaceColor','y');

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','g');                               

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in down sampling 1/4 frequency sample Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','P-wave','T-wave','R-wave','S-wave');

[~,locs_qwave] = findpeaks(ECG_data,'MinPeakHeight',-0.2);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_qwave,ECG_data(locs_qwave),'x','MarkerFaceColor','y');

% link and zoom in to show the changes

%linkaxes(ax(1:2),'xy');

%axis(ax,[60 230 0.006 -0.04])

%Next, we try and determine the locations of the Q-waves. Thresholding the peaks to locate the Q-waves results in detection of unwanted peaks as the Q-waves are buried in noise. We filter the signal first and then find the peaks. Savitzky-Golay filtering is used to remove noise in the signal.

smoothECG = sgolayfilt(ECG_data,1,3);

figure

plot(t2,ECG_data,'b',t2,smoothECG,'r'); grid on

axis tight;

xlabel('time msec'); ylabel('Voltage(mV)');

legend('ECG Signal','Filtered Signal')

title('Filtering Noisy of down sampling 1/4 frequency sample ECG Signal')

%We perform peak detection on the smooth signal and use logical indexing to find the locations of the Q-waves.

%[~,min_locs] = findpeaks(-smoothECG,'MinPeakDistance',29);

%[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',2);%Twave

[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',25);

% Peaks between -0.2mV and -0.5mV

%locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 &

%-smoothECG(min_locs)<-0.1); %Twave

locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 & -smoothECG(min_locs)<-0.11);

figure

hold on

plot(t2,smoothECG);

plot(locs_Qwave,smoothECG(locs_Qwave),'rs','MarkerFaceColor','g');

plot(locs_Rwave,smoothECG(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,smoothECG(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks down sampling 1/4 frequency sample in Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('Smooth ECG signal','T-interval','R-wave','S-wave');

%The above figure shows that the QRS-complex successfully detected in the noisy ECG signal.

%Error Between Noisy and Smooth Signal

%Notice the average difference between the QRS-complex in the raw and the detrended filtered signal.

% Values of the Extrema

[val_Qwave, val_Rwave, val_Swave] = deal(smoothECG(locs_Qwave), smoothECG(locs_Rwave), smoothECG(locs_Swave));

meanError_Qwave = mean((y1_1500(locs_Qwave) - val_Qwave))

meanError_Rwave = mean((y1_1500(locs_Rwave) - val_Rwave))

meanError_Swave = mean((y1_1500(locs_Swave) - val_Swave))

%% find PP interval

 i = 0;  %% to make the code start from 0.

 rr = 0; %% each time the code run, rr distance two peaks

 hold off % for the next graph

 rrinterval = zeros(3600,1); % create an array to strore 2 peaks

beat_count =0;

for k = 2 : length(y1)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y1(k)> y1(k-1) && y1(k) > y1(k+1) && y1(k)> 1);

     beat_count = beat_count +1;

     if beat_count ==1;

        rr =0;

     else

         rr = k-i;

         rrinterval(k)=rr;

         i=k;  

     end

    else

        rrinterval(k)= rr;

    end       

end

figure;

plot (rrinterval);

xlabel('Time in sec*10^-2'), ylabel('Distance betweeen 2 Heatbeats (R-R) in sec*10^-2'), title('R-R down sampling 1/4 frequency sample intervals');

%% find PP interval

%% heart rate analysis

% count the dominat peak

beat_count =0;

for k = 2 : length(y2)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y2(k)> y2(k-1) && y2(k) > y2(k+1) && y2(k)> 1)

         beat_count = beat_count +1;

    end       

end

display (k);

disp('dominant peaks');

%% divide the peak count by the duration in minute

duration_in_sec = N/fs2;

duration_in_minute = duration_in_sec/60;

BPM = beat_count/duration_in_minute;

%%%  DFT to describe the signal in the frequency

NFFT = 2 ^ nextpow2(N2);

Y = fft(y2, NFFT) / N2;

f = (fs2 / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz

amp = ( 2 * abs(Y(1: NFFT / 2+1))); % Vector containing corresponding amplitudes

figure;

plot (f, amp);

title ('plot single-sided amplitude spectrume of 1/4 frequency sample ECG signal')

xlabel ('frequency (Hz)')

ylabel ('|y(f)|')

grid on;

max_value=max(y1);

mean_value=mean(y1);

threshold=(max_value-mean_value)/2;

%% upsampling 2 sample frequency

close all;

clear all;

clc;

load('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.mat')

disp('Contents of workspace after loading file:')

whos

fs = 250; % find the sampling rate or frequency

fs2 = 250*2;

y1=xlsread('I:\BIOM_Signal_processing\Hw5\ECGsignal_1.xls');

T = 1/fs;% sampling rate or frequency

% find the length of the data per second

N = length(y1);

ls = size(y1);

t = (0 : N-1) / fs;% sampling period

figure; %subplot(1,2,1);

      plot(t,y1);

    %plot(x,y2, 'g');  

title ('plot of the original of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

% up sampling 2 of frequency sample

y2 = resample(y1,500,250);

N2 = length(y2);

ls2 = size(y2);

t22 = (0 : N2-1) / fs2;% sampling period

figure; %subplot(1,2,1);

      plot(t22,y2);

title ('plot of the up sampling 2 frequency sample of ECG signal')

xlabel ('time (sec)')

ylabel ('Amplitute (mv)')

grid on;

%% find P

m1=max(y2)*.60;

P=find(y2>=m1);

y1_1500 = y2(1:1850);

t2 = 1:length(y1_1500);

figure;

plot(t2,y1_1500);

title ('plot of subset of upsampling 2 frequency sample the ECG signal')

xlabel ('time (msec)')

ylabel ('Amplitute (mv)')

grid on

%% used the snip code from this website.

%%%%http://www.mathworks.com/help/signal/examples/peak-analysis.html

%Detrending Data

%The above signal shows a baseline shift and therefore does not represent the true amplitude. In order to remove the trend, fit a low order polynomial to the signal and use the polynomial to detrend it.

[p,s,mu] = polyfit((1:numel(y1_1500))',y1_1500,6);

f_y = polyval(p,(1:numel(y1_1500))',[],mu);

ECG_data = y1_1500 - f_y;        % Detrend data

N1= length (y1_1500);

t1 = (0 : N1-1) / fs2;% sampling period

figure

%plot(t1,ECG_data); grid on

plot(t2,ECG_data); grid on

ax = axis; axis([ax(1:2) -2.2 2.2])

%ax = axis; axis([ax(1:2) -3.2 3.2])

title('Detrended upsampling 2 frequency sample ECG Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

legend('Detrended ECG Signal')

%Thresholding to Find Peaks of Interest

%The QRS-complex consists of three major components: Q-wave, R-wave, S-wave. The R-waves can be detected by thresholding peaks above 0.5mV. Notice that the R-waves are separated by more than 200 samples. Use this information to remove unwanted peaks by specifying a 'MinPeakDistance'.

[~,locs_Rwave] = findpeaks(ECG_data,'MinPeakHeight',0.5,...

                                    'MinPeakDistance',240);

%Finding Local Minima in Signal

%Local minima can be detected by finding peaks on an inverted version of the original signal.

ECG_inverted = -ECG_data;

[~,locs_Swave] = findpeaks(ECG_inverted,'MinPeakHeight',0.4,...

                                        'MinPeakDistance',240);                                                            

%The following plot shows the R-waves and S-waves detected in the signal.

figure

hold on

plot(t2,ECG_data);

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

%axis([0 1850 -1.1 1.1]); grid on;

axis([0 1850 -2.2 2.2]); grid on;

legend('ECG Signal','R-waves','S-waves');

xlabel('time msec'); ylabel('Voltage(mV)')

title('R-wave and S-wave in upsampling 2 frequency sample of ECG Signal')

[~,locs_Twave] = findpeaks(ECG_data,'MinPeakHeight',-0.02,...

                                      'MinPeakDistance',100);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','y');                                

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in upsampling 2 frequency sample Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','T-wave','R-wave','S-wave');

[~,locs_Pwave] = findpeaks(ECG_data,'MinPeakHeight',-0.09,...

                                      'MinPeakDistance',52);

figure;

hold on

plot(t2,ECG_data);                              

plot(locs_Pwave,ECG_data(locs_Pwave),'x','MarkerFaceColor','y');

plot(locs_Twave,ECG_data(locs_Twave),'X','MarkerFaceColor','g');                               

plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks in upsampling 2 frequency sample Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('ECG signal','P-wave','T-wave','R-wave','S-wave');

[~,locs_qwave] = findpeaks(ECG_data,'MinPeakHeight',-0.2);

figure;

hold on

plot(t2,ECG_data);                             

plot(locs_qwave,ECG_data(locs_qwave),'x','MarkerFaceColor','y');

% link and zoom in to show the changes

%linkaxes(ax(1:2),'xy');

%axis(ax,[60 230 0.006 -0.04])

%Next, we try and determine the locations of the Q-waves. Thresholding the peaks to locate the Q-waves results in detection of unwanted peaks as the Q-waves are buried in noise. We filter the signal first and then find the peaks. Savitzky-Golay filtering is used to remove noise in the signal.

smoothECG = sgolayfilt(ECG_data,1,3);

figure

plot(t2,ECG_data,'b',t2,smoothECG,'r'); grid on

axis tight;

xlabel('time msec'); ylabel('Voltage(mV)');

legend('ECG Signal','Filtered Signal')

title('Filtering Noisy of upsampling 2 frequency sample ECG Signal')

%We perform peak detection on the smooth signal and use logical indexing to find the locations of the Q-waves.

%[~,min_locs] = findpeaks(-smoothECG,'MinPeakDistance',29);

%[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',2);%Twave

[~,min_locs] = findpeaks(smoothECG,'MinPeakDistance',25);

% Peaks between -0.2mV and -0.5mV

%locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 &

%-smoothECG(min_locs)<-0.1); %Twave

locs_Qwave = min_locs(smoothECG(min_locs)>-0.3 & -smoothECG(min_locs)<-0.11);

figure

hold on

plot(t2,smoothECG);

plot(locs_Qwave,smoothECG(locs_Qwave),'rs','MarkerFaceColor','g');

plot(locs_Rwave,smoothECG(locs_Rwave),'rv','MarkerFaceColor','r');

plot(locs_Swave,smoothECG(locs_Swave),'rs','MarkerFaceColor','b');

grid on

title('Thresholding Peaks down sampling 2 frequency sample in Signal')

xlabel('time msec'); ylabel('Voltage(mV)')

ax = axis; axis([0 1850 -2.2 2.2])

legend('Smooth ECG signal','T-wave','R-wave','S-wave');

%The above figure shows that the QRS-complex successfully detected in the noisy ECG signal.

%Error Between Noisy and Smooth Signal

%Notice the average difference between the QRS-complex in the raw and the detrended filtered signal.

% Values of the Extrema

[val_Qwave, val_Rwave, val_Swave] = deal(smoothECG(locs_Qwave), smoothECG(locs_Rwave), smoothECG(locs_Swave));

meanError_Qwave = mean((y1_1500(locs_Qwave) - val_Qwave))

meanError_Rwave = mean((y1_1500(locs_Rwave) - val_Rwave))

meanError_Swave = mean((y1_1500(locs_Swave) - val_Swave))

%% find PP interval

 i = 0;  %% to make the code start from 0.

 rr = 0; %% each time the code run, rr distance two peaks

 hold off % for the next graph

 rrinterval = zeros(3600,1); % create an array to strore 2 peaks

beat_count =0;

for k = 2 : length(y1)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y1(k)> y1(k-1) && y1(k) > y1(k+1) && y1(k)> 1);

     beat_count = beat_count +1;

     if beat_count ==1;

        rr =0;

     else

         rr = k-i;

         rrinterval(k)=rr;

         i=k;  

     end

    else

        rrinterval(k)= rr;

    end       

end

figure;

plot (rrinterval);

xlabel('Time in sec*10^-2'), ylabel('Distance betweeen 2 Heatbeats (R-R) in sec*10^-2'), title('R-R down sampling 2 frequency sample intervals');

%% find PP interval

%% heart rate analysis

% count the dominat peak

beat_count =0;

for k = 2 : length(y2)-1

    %the peak has to be greater than 1 and greater than the value before it and greater then the value after it.

    if(y2(k)> y2(k-1) && y2(k) > y2(k+1) && y2(k)> 1)

         beat_count = beat_count +1;

    end       

end

display (k);

disp('dominant peaks');

%% divide the peak count by the duration in minute

duration_in_sec = N/fs2;

duration_in_minute = duration_in_sec/60;

BPM = beat_count/duration_in_minute;

%%%  DFT to describe the signal in the frequency

NFFT = 2 ^ nextpow2(N2);

Y = fft(y2, NFFT) / N2;

f = (fs2 / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz

amp = ( 2 * abs(Y(1: NFFT / 2+1))); % Vector containing corresponding amplitudes

figure;

plot (f, amp);

title ('plot single-sided amplitude spectrume of upsampling 2 frequency sample ECG signal')

xlabel ('frequency (Hz)')

ylabel ('|y(f)|')

grid on;

max_value=max(y1);

mean_value=mean(y1);

threshold=(max_value-mean_value)/2;