Detecting QRS complex in ECG signal

Hi guys,

I have a few ECG mat files of USUSRP-AF patients with me. I was told that the signals have been sampled at 512Hz. I'm working on a project to prototype a PaceMaker. Therefore I need to filter these ECG signals to find the QRS complex by applying the necessary filters as appropriate. I'm very new to this MATLAB environment. What filters I need to use (lowpass, hihgpass, bandpass)? Please advice and appreciate if someone could include the commands to accomplish this task.

I have tried few examples in search results from Google. However, there wasn't a complete solution or clear instructions or guidance on how to do it. Therefore, I decided to post here so that an expert in this area would reply.

Thanks for your precious time.

Cheers, Rumi

Answer 1

[Star Strider]

When you say ‘EKG’ do you mean full surface 12-lead, Lead-II rhythm strips, or intracardiac recordings (taken with an electrophysiology catheter)? Atrial fibrillation (AF) and most other rhythm disturbances are rarely isolated pathology. The ventricular (QRS) rhythm is by definition irregularly-irregular, with a quite uneven baseline.

I certainly suggest a bandpass filter, but I’ve never analysed an AF EKG from a signal processing perspective, so I’m not certain what the best passband would be. The low end must be high enough to exclude baseline drift, and the high end low enough to filter out the atrial signal and of course exclude 50-60 Hz power frequencies. Ventricular rhythms can vary from around 50 bpm (0.8 Hz) to 300 bpm (5 Hz). I would suggest a 0.20 Hz low-frequency passband limit but your sampling frequency will only allow you at most 250 Hz at the high end.

You need to do FFTs of at least Lead II of each EKG record to see what the spectra have in common. (Posting two or more of your randomly-selected EKG mat-files, preferably of full 12-lead EKGs here would help.)

[Mohamed]

Hi, thanks for the swift response. I can see 12 columns in the data matrix in the .mat files. So I think it is a 12-lead recording. I have limited knowledge in the ECG domain. I have attached the mat file which contains recordings of four patients. It would be great if you could take a look at these and suggest the best option and help me with the MATLAB commands. Thanks for your assistance.

[Star Strider]

I’m interested in what you’re doing.

No mat files appeared yet. (Use the ‘paperclip’ or ‘staple’ icon to upload the files.)

I'll keep this open and keep checking back.

[Mohamed]

I'm trying to attach a 3MB zip file. But the Attach button doesn't respond when I select the file and press the button. Any idea why?

[Mohamed]

Ok, I think the file have now been attached :)

[Star Strider]

The files arrived!

We’ve been inundated by spam of late, so that could be part of the problem — a sort of DOS attack. I’ve noticed that it’s been slow for the past few days.

I’ll take a look at the files and reply here. Could take some time, depending on what I need to do.

[Mohamed]

No worries. I understand. Take your time and thank you very much for your kind help.

[Mohamed]

Hi,

Appreciate you're busy. Did you get a chance to take a look at those files?

Thanks.

[Star Strider]

Yes.

I have the routine written, the filtering was straightforward, but when I plotted the 12-lead EKGs discovered that all the precordial ,leads ( V_1 - V_6 ) are the same. This is not physiological, so whoever compiled these had no knowledge of electrocardiology or electrocardiography or they’d have caught the error. I believe the repeated precordial lead tracings are lead V_1 ( they look like V_1 ), but I don’t know for sure. The problem with this is that the tip of an transvenous pacemaker electrode essentially ‘sees’ lead V_3 .

I generated vectorcardiograms to make the EKGs easier to read. The problem comes in ‘thresholding’ the R-waves in order to identify the QRS complexes.

I attached the filtered 12-lead EKGs in standard configuration. They are difficult to read in terms of identifying other abnormalities, specifically ventricular conduction defects. (These are very sick hearts.)

Everything else is working. I need to be satisfied with the robustness of my QRS detection criteria, difficult because most of these recordings display electrical alternans, then I’ll clean up my code, put it in a function statement, and post it. I’ll leave the argument list and outputs of the function undefined so you can define them to return what you want. I also have to finish comment-documenting it so that its function and logic are transparent.

I’ve been working on it most of the week, and intend to post it as an attachment tomorrow. It’s been quite an adventure!

With normal EKGs this would have been a breeze. But normal EKGs don’t need therapeutic intervention.

[Star Strider]

——————————————————— COMPLETED! ———————————————————

I attached the function file ‘EKG_Explore.m’ to my original Answer. I’m satisfied that it does everything I can make it do, considering that the hearts in the EKG files have significant ventricular conduction defects and other abnormalities. The presence of the atrial fibrillation ‘noise’ prevents me from characterising the conduction defects and cardiac disease in any detail.

The function itself first filters out baseline drift and power-line interference with a bandpass filter, then detects the R-wave onset, and generates a vector of indices for the R-wave onsets. It plots and prints-to-PDF the full 12-lead scalar EKG tracings in the conventional format. (MATLAB does not plot them in conventional EKG order, so indexing them correctly was an interesting exercise.) Since EKG reading is something of an art, I also plotted them as frontal (coronal) and 3-D spatial vectorcardiograms to make them easier to interpret. I broke down the segments into the PR interval, QRS complex, and ST segment by colour.

I left the argument list and output undefined so you can fill them as you wish. I’m not certain what you want it to take as arguments or return as output. (Right now, it takes no arguments and returns no output.) If you want to suppress any of the plots, you can loop around them with an ‘if’ block.

I did my best to ‘comment-document’ it to make its logic and operation as transparent as possible.

Although it took more hours than I thought it would to program and get running as I wanted it to, I actually had fun doing it!

[Mohamed]

Hi,

You've done an excellent job and thanks for taking time out of your personal busy routine and helping me on this. I'm not very sure I understand all what you have said in the reply as my knowledge in this area is very very limited. However, I will first download the files you have attached and have a read and try to understand. Finally, if there is anything I still don't understand or need clarification, I will come back to you if that's OK with you.

Once again thanks for the time and effort.

Appreciate a million times!

Cheers, Rumi

[Star Strider]

My pleasure!

That’s fine. EKG analysis and interpretation is definitely not a trivial exercise. I did my best to make my function as transparent at possible.

Start by running the file to see what it produces, then read through it. You may want to add input arguments and outputs to it, but I need to know what those are if you want me to make any changes to the code.

EDIT — If it does what you want it to, the highest expression of appreciation here on MATLAB Answers is to Accept the Answer that most closely solves your problem.

[Mohamed]

Hi,

Definitely, I will accept the answer you have submitted. There is no doubt in that. But, I need couple of things to clarify before I close it.

As I mentioned the very first time I posted this question, the main requirement here is to produce an output (the heart beat signal) from matlab which can then be read by my micro-controller. Correct me if I'm wrong, the heart beat signal will be 'R' from 'QRS' segment right? Then the idea idea is to detect any missing beat sequence let's say 4 second delay and then my FPGA (hardware) will decide whether or not to produce a shock. So, in order to achieve this I should produce the heart beat signal sequence as output from matlab. Not sure if I'm making sense.

If you understood what I have just explained and and you know exactly what I'm looking for, would it be possible to change your function to produce this output. Again in simple words my hardware will read the heart beat signal sequence from the matlab output continuously and will produce a shock when missing beat was detected and last for at least 4 seconds.

Thanks Rumi

[Star Strider]

Correct. It will be the R-wave from Lead II (or Lead V_3 in a complete EKG, but not here since as I mentioned before, the precordial ‘V’ leads are in error in these records). The EKG looks different in different leads, so it is necessary to specify the lead you are monitoring.

The full 30-second heart beat sequence is the EKGf array. It is the bandpass-filtered EKG signal for all 12 leads.

The ‘RwvIdxAll’ array has the indices of all the R-waves the algorithm detects. (Because of the atrial fibrillation, I used a window of 0.25s to be sure I was detecting only R-waves.) Figure 4 indicates it detects all of the R-waves, and nothing else.

If you want the ‘RwvIdxAll’ array in seconds in addition to index values, multiply it by the sampling interval Ts=1/Fs, with Fs=512Hz. I can do that in the function and return it as an output of the function if you want.

These EKGs all have some degree of tachycardia, and none of the R-R intervals is as long as 4 seconds. I can add statements to the function to output the instantaneous rate and R-R interval if you want. (Those vectors will be one less than the length of the R-wave vector, by definition.)

[Mohamed]

Well, then do you mean to say that I cannot use these ECG data files for my project since the data from leads has errors? How do I specify the lead I want to monitor? If you think I can still use these files, which one out of the four patient data is more appropriate?

Well, the idea is to remove (delete) some data, lets say 4 seconds worth of data to simulate missing beat and then to demonstrate the hardware detects it and takes action (produces shock) without delay.

I know you're an expert in this area and you have the capacity to produce any output I want :)

All I need is an output of the 'R' sequence. For example '1' bit when R was detected and continuous '0' bits all other time. This can then be read by my hardware and I can calculate the interval and the rate in my micro-controller and hence make instantaneous decision when things go wrong.

Please advise me if you think I should follow a different approach.

Cheers, Rumi

[Star Strider]

These are certainly appropriate files if you are dealing with sick hearts! You can use the first six leads (Leads I-aV_F), but the precordial ( ‘V’ ) leads are not reliable in these tracings.

In the EKGf matrix, the second dimension is the lead (Lead II for patient k is EKGf(:,2,k), and Lead aV_F is EKGf(:,6,k)).

The indicator (0,1) vector indicating the presence of an R-wave is easy to create. I’ll add it to the output list.

All the difficult work I’ve finished, and everything you want I can easily create from the data I’ve already calculated.

Anything in particular you want as input arguments?

----------

(As for expertise, I am a physician (M.D.), Diplomate of the American Board of Internal Medicine, and I have a M.S. in Biomedical Engineering.)

[Mohamed]

OK, I think it all makes more sense to me now. So from the knowledge you have shared with me so far what I've understood is, R wave can be detected or read either from Lead-II or from pre-cordial. But, due to the error in the 'V' leads you mentioned earlier in this situation it is not possible to the pre-cordial readings. So, you have used the Lead-II (I-aV_F) to detect the 'R' wave. Is my understanding correct?

OK, so how about setting the 'patient id' and the signal/wave ('letter') as the input arguments? Do you think they are appropriate and valid input arguments in this context?

So for example, when (2,Q) are passed as input arguments, it means to process the patient 2's data and return the Q wave as output. In this manner is it possible to specify the letters P,Q,R,S,T as the output signal?

Sorry, if I look too fussy :). If this is too much work and/or you think makes no sense, leave 'R' as the fixed output signal.

Finally, yes so if you could output the vector (0,1) to indicate the presence/absence of the signal/wave would be really great.

P.S: There is no question, all the knowledge you've shared with me so far is going to help me a lot in my project.

I'm glad you're apart of this community which is helping learners like us. Keep it up.

Hopefully, after your next answer I will be done and will accept your answer. YOU'RE A LEGEND!

Cheers, Rumi

[Star Strider]

Yes. I use Lead II because it is the standard lead for R-wave detection. (It is aligned with the normal QRS axis in the frontal plane.) I would prefer Lead V_3 but it is not reliable in these records.

The only input argument is ‘Patient Number’. I added a bit of logic to be sure it’s valid (1-4) here. The binary (0,1) vector is now the only output. You can easily add others if you want. The most relevant would likely be ‘Tv’. It is the vector of time (seconds) for all events in every record.

I have removed the previous version and attached the new completed version it its place

Call it as:

Rwave01 = EKG_Explore( PtNr );

The ‘Rwave01’ (or whatever you want to call it in your application) is the vector of R-wave occurrences where the presence of an R-wave is ‘1’ and it is ‘0’ elsewhere.

Signal/wave is impossible. You are pretty much stuck with the R-wave (ventricular action potential) in these and most other tracings. Particularly, P-waves are absolutely nonexistent in all these records because all the patients have atrial fibrillation (AF). Q-waves are always less than 1 millivolt unless there are pathological Q-waves indicative of a transmural myocardial infarction in the leads opposite the infarct. Also, they are going to be extremely difficult to detect in the presence of AF. The S-wave is also difficult to detect in most leads; it is the end of the intrinsicoid deflection, although it can be seen in some of these tracings as a component of the RSR’ pattern and indicative of right bundle branch block, right ventricular hypertrophy, and some other conditions, alone or in combination. (Considering that mitral stenosis can cause both AF and right ventricular hypertrophy, that would be important to rule out in these patients.) It is also essentially impossible to determine the T-wave here because of the AF. However it is extremely important to detect the T-wave in your application (or infer its approximate location) to avoid a stimulus spike on the descending slope of the T-wave (‘ventricular vulnerable period’), since stimulating here will usually result in the onset of ventricular tachycardia that unless treated quickly will progress to torsade-de-pointes and frequently ventricular fibrillation. That is counterproductive for a pacemaker!

I very much appreciate your enthusiastic support of my efforts! This was actually a fun project for me.

[Mohamed]

Ok, I think you have made a huge effort in getting this right and made valuable explanation as necessary to make the whole process much clear.

One last question out of interest. I'm sure I will stop at this point. When we say heart BPM whihc segment are we counting is it the QRS or just the 'R' wave.

And which corresponds to pulse rate? Or are both BPM and pulse rate the same?

Cheers, MI

[Star Strider]

BPM is by definition the number of QRS complexes per minute (particularly the R-wave, Lead II being most representative) because the QRS complex represents ventricular depolarisation, and that generates the arterial pulse wave. (It is possible in extreme situations to see normal cardiac electrical activity but no corresponding arterial pulse. This is called ‘electromechanical dissociation’ and is almost always an irreversible, agonal event.) So in most situations, QRS complexes per minute correspond to arterial pulse waves per minute.

Some of the data I generate are beats/second and beats/minute, calculated on a beat-to-beat basis. (The first element in those arrays is ‘NaN’ to make the vector lengths come out the same as the others.) You can add those to the output arguments if you want them. They are clearly defined in the documentation.

The part of the QRS complex you detect depends on the lead you are considering. Normally, this is usually the Lead II R-wave because Lead II is most closely aligned with the normal QRS axis.

The EKG is actually the cardiac electrophysiologic cycle projected on the various EKG leads. Think of it as a vector with the origin at the SA node, varying in magnitude and orientation as it rotates through the heart over time. The VKG captures the entire loop, and illustrates the variation of the QRS complex with respect to the lead you are using:

In these patients, it looks a bit strange because of the presence of atrial fibrillation. Normal VKGs have three distinct loops: atrial depolarisation (P-R), ventricular depolarisation (QRS), and ventricular repolarisation (S-T). The presence of atrial fibrillation obliterates the P-wave and obscures the T wave.

[Mohamed]

Thanks for the detailed explanation and appreciate your time and thanks for sharing your valuable knowledge. I would like to keep in touch so that I can seek assistance from you if I need any help in future. Would it be possible for you to share your email id if that is OK with you?

Finally I have accepted your answer. Well done!

Cheers, Rumi

[Star Strider]

My pleasure!

Cardiac electrophysiology — and cardiovascular physiology in general — takes time and effort to learn and understand, so I don’t mind getting you started on the correct track. For details, I refer you to the classic text on the subject, Braunwald’s Heart Disease, 7th Ed., Ch. 9, Electrocardiography, available free online. The endocardial EKG your transvenous pacemaker catheter electrodes ‘see’ is going to be different, so I suggest you search PubMed for those references.

You can find my contact information in my File Exchange profile. That may change over time, so I won’t post it here.

[Mohamed]

Great and thanks!

[Star Strider]

Again, my pleasure!

[Mohamed]

Hi,

Could you clarify the below when you get a chance pleasE?

1. Why you have used a function like tf2sos in the script an what are sos and g represent? I don't think I understand the text available in the help very well

2. Why have you decided to use 5Hz and 45Hz as the low and high as the band in the bandpass filter?

3. What does the following statement mean in the formatting section? Line 88: TLimV = fix(Fs*0.2):fix(Fs*2.7); Line 91: idx12 = min((3*(k1-1)+1), mod((3*(k1-1)+1),12)) + (k1>4) + (k1>8);

4. Could you explain me logic below

% DETECT ‘ACCEPTABLE’ LEAD II QRS COMPLEXES & CREATE TIME WINDOW TO EXCLUDE ANY REMAINING NOISE

RwvIdx0 = find( (abs(EKGf(:,2)) >= mean(abs(EKGf(:,2)))+1*std(abs(EKGf(:,2)))) );

5. On what basis you have selected 0.25*Fs for detecting the QRS complex?

RwvIdxAll = RwvIdx0(find(diff([0; RwvIdx0]) > 0.25*Fs));

Cheers, Rumi

[Star Strider]

I’ll do my best!

1. Transfer functions can have stability problems when applied to filters, especially with higher orders or narrow passbands. Second-order-section representations (sos) generally do not.
2. Much of biomedical signal processing is in some sense empirical or heuristic. The lower frequency 5 Hz limit eliminates the baseline drift, The higher frequency 45 Hz gets all the important detail in the EKG while eliminating the possibility of mains-frequency (50-60 Hz) interference.
3. That has to do with plotting a few QRS complexes of the EKG in the 12-lead scalar representation. It chooses a 2.5 second section of the EKG to plot. I did that primarily for my convenience, since that is the sort of time window most 12-lead EKG tracings provide. It’s only used in figure(3).
4. In order to detect the first deflection in the QRS complex of quite abnormal EKGs, I used the threshold of one standard deviation from the mean as the criterion. I used the absolute value of the QRS complex to be certain I detected the initial deflection, since the orientation varies between leads and different cardiac pathologies.
5. The 0.25*Fs value (actually 0.25 sec) is a time window (in samples) to be sure noise in that interval did not detect a QRS onset event. Note that this is significantly less than the S-T interval (usually less than 0.4 sec) and twice what is considered to be a pathological QRS complex (greater than 0.12 sec). Such time windowing is a common practice in biomedical signal processing.

Did this answer your questions?

[Mohamed]

Yes, I think you have explained it very well.

I now need to send this output to my hardware for further processing. I found some examples where the serialPort can be used for this purpose.

I can certainly sort out the code to write the R vector values {0,1} to the serialPort. My question is, how can I write this R vector values in a timely manner? I mean the detection of the R ('1' or '0') should reflect the actual time in occurrence, shouldn't it? Only then I can take action if R was missing for a period of time.

Is the time interval between two value in the R vector be 0.25 sec?

Sorry for throwing too many questions at you and bothering you.

Cheers, Rumi

[Star Strider]

Programming FPGAs is not in my areas of expertise. I cannot help you with the hardware programming or communication. I interfaced a USB fingertip pulse-oximeter using only the core MATLAB functions a few years ago, but that is the limit of my hardware interfacing experience.

The detection of the R-wave should definitely reflect the time it actually occurs. How long you decide to wait before you generate a pacemaker stimulus is an important design consideration. You certainly need to take into account the underlying cardiac rate so that you know when one or more consecutive beats have been ‘missed’. The normal resting heart rate can be as low as 50 beats/min in well-conditioned individuals.

the 0.25 second criterion was a detection window in my QRS detection algorithm. I would not use it in the design of your pacemaker, since you would not want to pace at such a rapid rate (that equates to 240 beats/min) anyway. If you sense a rate that high, your pacemaker should not do anything.

If you can locate a copy of Braunwald, it probably has a chapter on cardiac pacing. That could help you in your design.

Have fun with this!

[Mohamed]

Not a problem. I understand it very well. Anyway thank you very much for helping on this so far.

Cheers, Rumi

[Star Strider]

Again, my pleasure!

[Mohamed]

One final... last question.

I noticed in one of the figures where you are marking the detection of R wave with a '+' sign. The question is can this be both +/-. I see the marking are +/- in the figure. That's all. I won't bother you for sure on this question again.

Thank you. Rumi

[Star Strider]

I plotted the first detected deflection in each QRS complex with a red ‘+’ sign simply to make it visible. The ‘+’ sign is a marker, and not a sign indicator here.

[Mohamed]

Sorry for not being very clear. My question is not about the '+' sign itself, but the value (voltage) represented by the marking on the wave? Can it be + or minus on the y axis?

Sorry, if my question is not making sense to you.

I understand the fact that a sign wave can have values on the 'y' axis which is either '+' or '-'. I wanted to clarify if that applies to 'R' as well.

Cheers, Rumi

[Star Strider]

Yes, it can be + or - as plotted on the actual EKG.

The EKGs in your data not only have atrial fibrillation, they also have varying conduction system abnormalities, so I used the absolute value of the QRS complex to be certain I detected the initial deflection, regardless of the direction. (I needed the initial deflection so I could estimate the P-R, QRS, and S-T segments appropriately.)

[Mohamed]

Thank you very much for your help throughout and for the useful explanations. I will leave you in peace.

Finally, I have to appreciate the fact that service of professionals like you will have a great impact on the community.

Good stuff! Keep it up :)

Cheers, Rumi

[Star Strider]

My pleasure, as always.