MATLAB Program for ECG analysis You are to design a series of flters to enhance

Question

MATLAB Program for ECG analysis

You are to design a series of flters to enhance an ECG signal that is corrupted

with noise from the typical 60 Hz. electrical AC sources. This can be done by

applying a series of notch filters, three should work, to remove the fundamental 60 Hz.

noise and its first two harmonics, 120 Hz and 180 Hz.

Next we wish to perform a calculation to determine the heart rate of the sample

data, by using detection of zero crossings. Since the original signal also has

"noise" from sources internal to the body, such as muscle activiations, it will be

necessary to apply a bandpas filter. It will suffice to use a 4th order Chebyshev

filter with a bandpass region of 0.25 - 40 Hz. Use a prototype 2nd order Chebyshev

low pass filter ( supplied ), pre-warp the frequencies in code ( do not use the option

in the bilinear function to do this automatically,), use the lp2bp function followed

by the bilinear function to arrive at a numberator and denominator polynomial.

Then apply each of the individual notch filters followed by the bandpass filter

to the supplied data file ecgbn.dat. Plot the original signal, the signal after

removal of the 60 Hz. and harmonic noise, and then the signal following the bandpass

filter ( use subplot for a 3x1 plot). Ttile each plot. Plot the zeros and poles for

each filter and title each plot. Scale axes if needed. Plot the magnitude and phase

for all filters.

For the notch filters use the pole placement design method, code a function to return

the numerator and denominator coefficients given a sampling rate, 3 dB bandwidth,

and passband center ( explicitly coding this for each section will result in loss

of points. Can use and external .m file or simply place the function at the end

of your code. )

Use good coding style, which means identify all variables used at the beginning of

your code, and use those variables in function calls. Directly coding information

such as sampling rate and frequency band edges in function calls is bad style, leading

to higher chances of making errors if (when) code is re-used.

Count the number of zero-crossings using the function provided.

1. The second order notch filter equations will be provided. These were derived using a simple design method called “pole placement” which is suitable for filters

that have a narrow pass or stop band.

2. The notch filters should be designed to remove 60 Hz., 120 Hz, and 180 Hz frequencies. Use a 4 Hz. 3-dB bandwidth and a 600 Hz. Sampling rate.

3. Design the bandpass filter using the supplied analog 2nd order Chebyshev prototype, apply pre-warping ( do not use the bilinear functions option here ),

the use MATLAB’s lp2bp function, followed by the bilinear function to obtain

numerator and denominator coeficients.

4. Apply each filter to the supplied ecgbn.dat file. Plot the original data, the data after removal of the 60 Hz. Frequency and harmonics, and the data after all filters

have been applied. Use the zplane function to plot poles and zeroes of each filter, along with their magnitude and phase respone.

5. Run the supplied zero crossing function to determine the heart rate of this

     sample data.

6. Answer these questions:

Question one. There is a programming principle called “DRY”, or Do not Repeat Yourself”. This means it is poor programming practice to have multiple sections of code perform the same function. Look up on the Internet about DRY and write a SHORT paragraph on how it applies to this program. This of why I want you to write a function to perform the second order notch calculations and why I don’t like directly placing frequency data into function calls.

Question two. The original data has a bandwidth of 0.1 – 250 Hz, and the resulting data has been bandpass filtered. Comment on the suitability of using the “cleaned” data for medical analysis. What steps in this processing chain can interfere with proper analysis and how could you modify this sequence of filters to provide data more suitable for analysis.

Explanation / Answer

x=0.01:0.01:2;
default=input('Press 1 if u want default ecg signal else press 2: ');
if(default==1)
li=30/72;  
  
a_pwav=0.25;
d_pwav=0.09;
t_pwav=0.16;  

a_qwav=0.025;
d_qwav=0.066;
t_qwav=0.166;
  
a_qrswav=1.6;
d_qrswav=0.11;
  
a_swav=0.25;
d_swav=0.066;
t_swav=0.09;
  
a_twav=0.35;
d_twav=0.142;
t_twav=0.2;
  
a_uwav=0.035;
d_uwav=0.0476;
t_uwav=0.433;
else
rate=input(' enter the heart beat rate :');
li=30/rate;
  
%p wave specifications
fprintf(' p wave specifications ');
d=input('Enter 1 for default specification else press 2: ');
if(d==1)
a_pwav=0.25;
d_pwav=0.09;
t_pwav=0.16;
else
a_pwav=input('amplitude = ');
d_pwav=input('duration = ');
t_pwav=input('p-r interval = ');
d=0;
end   
  
  
%q wave specifications
fprintf(' q wave specifications ');
d=input('Enter 1 for default specification else press 2: ');
if(d==1)
a_qwav=0.025;
d_qwav=0.066;
t_qwav=0.166;
else
a_qwav=input('amplitude = ');
d_qwav=input('duration = ');
t_qwav=0.166;
d=0;
end   
  
  
  
%qrs wave specifications
fprintf(' qrs wave specifications ');
d=input('Enter 1 for default specification else press 2: ');
if(d==1)
a_qrswav=1.6;
d_qrswav=0.11;
else
a_qrswav=input('amplitude = ');
d_qrswav=input('duration = ');
d=0;
end   
  
  
  
%s wave specifications
fprintf(' s wave specifications ');
d=input('Enter 1 for default specification else press 2: ');
if(d==1)
a_swav=0.25;
d_swav=0.066;
t_swav=0.09;
else
a_swav=input('amplitude = ');
d_swav=input('duration = ');
t_swav=0.09;
d=0;
end   
  
  
%t wave specifications
fprintf(' t wave specifications ');
d=input('Enter 1 for default specification else press 2: ');
if(d==1)
a_twav=0.35;
d_twav=0.142;
t_twav=0.2;
else
a_twav=input('amplitude = ');
d_twav=input('duration = ');
t_twav=input('s-t interval = ');
d=0;
end   
  
  
%u wave specifications
fprintf(' u wave specifications ');
d=input('Enter 1 for default specification else press 2: ');
if(d==1)
a_uwav=0.035;
d_uwav=0.0476;
t_uwav=0.433;
else
a_uwav=input('amplitude = ');
d_uwav=input('duration = ');
t_uwav=0.433;
d=0;
end   
  

  
end

pwav=p_wav(x,a_pwav,d_pwav,t_pwav,li);

%qwav output
qwav=q_wav(x,a_qwav,d_qwav,t_qwav,li);

  
%qrswav output
qrswav=qrs_wav(x,a_qrswav,d_qrswav,li);

%swav output
swav=s_wav(x,a_swav,d_swav,t_swav,li);

%twav output
twav=t_wav(x,a_twav,d_twav,t_twav,li);

%uwav output
uwav=u_wav(x,a_uwav,d_uwav,t_uwav,li);

%ecg output
ecg=pwav+qrswav+twav+swav+qwav+uwav;
figure(1)
plot(x,ecg);

Related Questions

Navigate

• Browse (All)
• Browse M
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.