Matlab code to study the effects of noise in ECG signals

The goal of this assignment is to examine the effects of noise in signals.

1)      Simulate and plot 2 minutes of white noise, at the sampling rate in your ECG signals. 

2)      Add the noise to the ECG signal that was assigned to you. Do this in such a way that the signal-to-noise ratio (SNR) is 9 dB. It may be helpful to remember that the SNR can be calculated as10log10(σs2σn2), where σs is the standard deviation of the signal and σn is the standard deviation of the noise

 

3) Using the algorithm you designed to detect PQRST, examine the results of the system with the additive noise.  

Matlab code:

ECE/BIOM 537: Biomedical Signal Processing

Colorado State University                                                           

Student: Minh Anh Nguyen

Email: minhanhnguyen@q.com

close all;

clear all;

clc;

%%Simulate and plot 2 minutes of white noise, at the sampling rate in your ECG signals.

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

fs = 250 % find the sampling rate or frequency

T = 1/fs;% sampling rate or frequency or Time infos (for white noise)

N = length(y1); ls = size(y1);% find the length of the data per second

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

% White noise excitation

tspan = 0:T:120; % I saw from the plots in the paper that the simulation last 120 seconds

tstart = tspan(1); tend = tspan (end);

% generate Gaussian (normally-distributed) white noise

noise = randn(length(tspan), 1);

figure; subplot(211), plot (tspan,noise);

title ('plot 2 minutes of white noise at sample rate of 250Hz   ');

xlabel ('time (sec)');ylabel ('Amplitute'); grid on;

S = std(noise)

nbins = 150;

subplot (212), hist(noise,nbins),title ('histogram of 2 minutes of white noise at sample rate of 250Hz');

max_time = 120;    % Duration of your signal in seconds.

tmax = T:T:max_time;    % This is our time vector.

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

fignum = 0;

load('\BIOM_Signal_processing\Hw5\ECGsignal_1.mat') % contains hr_sig and fs

% Make time axis for ECG signal

tx = [0:length(hr_sig)-1]/fs;

fignum = fignum + 1; figure(fignum)

plot(tx,hr_sig), xlabel('Time (s)'), ylabel('Amplitude (mV)'), title('Zoom into original ECG signal'), xlim([30.3,31]) % Used to zoom in on single ECG waveform

%% add white noise

figure; SNRdb = 9; %Given SNRdb or Target SNR = 9 dB

variance = 1/(T*(10^(SNRdb/10)));

msgbox(strcat('variance of noise= ',mat2str(variance),''));

y = awgn(y1,SNRdb,'measured','linear');

subplot(211), plot(t,y1,t,y); % Plot both signals.

legend('Original signal','Signal with AWGN'); title('ECG signal with white noise at SNR = 9dB'), xlabel ('time (sec)'), ylabel ('Amplitute (mv)'),

subplot(212), plot(t,y); xlabel('Time (s)'), ylabel('Amplitude (mV)'), title('Zoom into ECG signal with white noise'), xlim([30.3,31]) % Used to zoom in on single ECG waveform

W = sqrt(variance).*randn(1,size(y1,2)); %Gaussian white noise W

y1_noise = y1 + W; %Add the noise

figure; subplot (3,1,1), plot(t,y1_noise,t,y); % Plot both signals.

legend('Original signal','Signal with AWGN'); title('ECG signal with white noise at SNR = 9dB'),xlabel ('time (sec)'), ylabel ('Amplitute (mv)'), %% add white noise

%% noise and DFT

s_noisy = awgn(y1,SNRdb);

NFFT = 2 ^ nextpow2(N); Y = fft(s_noisy, 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

subplot(3,1,2), plot (f, amp);title ('plot single-sided amplitude spectrume of the ECG signal with white noise'); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;

% Add noise

ECGwithNoise = awgn(y1,SNRdb,'measured');

figure; plot(t,ECGwithNoise); xlabel('Time (s)'), ylabel('Amplitude (mV)'), title('Zoom into ECG signal with white noise'), xlim([30.3,31]) % Used to zoom in on single ECG waveform

% Seperate noise

noise = ECGwithNoise - y1;

% Compute SNR input

snri = 10*log10 ( sum(ECGwithNoise.^2)./sum(noise.^2) );

S2 = sum(y1.^2); % signal power

N2 = S2*10^(SNRdb/10); % noise power

 std_ECG= std2(y1);

 msgbox(strcat('the standard deviation of original ECG= ',mat2str(std_ECG),''));

 std_ECG_noise= std2(ECGwithNoise);

 msgbox(strcat('the standard deviation of ECG with noise= ',mat2str(std_ECG_noise),''));

%% create a subset to zoom into the signal make easy to verify mark position

y1_1500 = ECGwithNoise(1:1850);

t2 = 1:length(y1_1500);

figure; plot(t2,y1_1500);

title ('plot of subset of the ECG signal with white noise')

xlabel ('time (msec)'), ylabel ('Amplitute (mv)'); grid on

%c) Write code to automatically detect the various features of the ECG (PQRST) and use that to mark the ECG waveform features

%% 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(t2,ECG_data); grid on

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

title('Detrended ECG Signal with white noise')

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

legend('Detrended ECG Signal with white noise')

%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 -2.2 2.2]); grid on;

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

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

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

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

                                      'MinPeakDistance',40);                         

%% The following code detect and mark T                                

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 ECG Signal with white noise')

xlabel('time msec'); ylabel('Voltage(mV)'); ax = axis; axis([0 1850 -2.2 2.2]); legend('ECG signal with white noise','T-wave','R-wave','S-wave');

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

                                      'MinPeakDistance',20);                           

%% The following code detect and mark P                               

figure; hold on; plot(t2,ECG_data);

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

plot(locs_Pwave,ECG_data(locs_Pwave),'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 ECG Signal with white noise')

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

ax = axis; axis([0 1850 -2.2 2.2]), legend('ECG signal with white noise','P-wave', 'T-Wave','R-wave','S-wave');                             

%% The following code detect and mark Q

ECG_inverted = -ECG_data;

[~,locs_Qwave] = findpeaks(ECG_inverted,'MinPeakHeight',0.06,...

                                     'MinPeakDistance',22);                                                           

figure; hold on, plot(t2,ECG_data);

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

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

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

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

plot(locs_Qwave,ECG_data(locs_Qwave),'o','MarkerFaceColor','b');

grid on, title('Thresholding Peaks in ECG Signal with white noise')

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

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

legend('ECG signal with white noise','P-wave', 'T-Wave','R-wave','S-wave', 'Q-Wave');