ECG processing

The trace of a heartbeat rarely arrives clean, and each kind of noise needs its own filter. Watch the solid trace below settle onto the dashed clean one as the error falls, stage by stage.

One filter for each kind of noise

Put electrodes on someone’s chest and they pick up a few millivolts from the heart. Each heartbeat draws the same shape. There is a small bump, then a sharp spike with a dip on each side, then a broad hump. This trace is an electrocardiogram, or ECG.

The small bump is the P wave. The spike and its two dips are the QRS complex: Q is the first dip, R the tall spike, S the second dip. The broad hump is the T wave. The R spike is the easiest part to find, so the time of a beat means the time of its R peak.

A real ECG rarely arrives clean. Three kinds of noise ride on it, and each lives in its own range of frequencies:

• Baseline wander is a slow drift of the whole trace, below about 0.5 Hz. Breathing and electrodes that shift on the skin cause it.
• Mains hum is the 50 Hz pickup from the power supply that you met in “Notch out the mains hum” of Resonators, notches and combs (17.4).
• Muscle noise comes from the electrical activity of other muscles. It is a fuzz spread over a wide range of frequencies.

Since the three noises sit in different bands, one filter can take out each. Think of sieving sand through three meshes, each made for one size of grain. This page builds the three filters, then finds the beats and reads the heart rate.

Everything here is a teaching model, not a clinical tool. The ECG is synthetic, built from formulas, and the page makes no diagnosis of any kind.

A synthetic ECG

I build each beat from five bell-shaped waves, one for each of P, Q, R, S and T. Each wave of beat $k$ is

where $t_k$ is the beat’s R time in seconds, $a$ the wave’s height, $c$ its centre measured from $t_k$, and $b$ its width. This is the bell of Random variables for signals (24.1), used here as a shape. The five waves are:

The time from one R peak to the next is $T_\text{RR}$, in seconds. A real heart does not keep perfect time, so I let it vary a little. The first beat is at 0.4 s, and each next one comes

later, where $u$ is a seeded random number between 0 and 1 (seed 311). I keep adding beats while they land before 9.6 s, 0.4 s from the end, which gives 12 beats in the 10 s record. The intervals run from 0.772 s to 0.826 s.

The record is sampled at $f_s=360$ Hz, the rate of the MIT-BIH recordings that ECG methods are often tested on, so it holds 3600 samples.

The three noises are formulas too. The baseline wander is

a breath every 4 s and a slower drift. The hum is $0.1\sin(2\pi\cdot50\,t+0.3)$ mV. The muscle noise is 0.03 mV times a seeded Gaussian draw for each sample (seed 3110). Its RMS, the size of How big is a signal (1.3), comes out at 0.0298 mV.

Measuring what is left

To score each filter I compare its output with the clean ECG, which only a synthetic record can give me. The score is the RMS of the difference, over seconds 1 to 9. I leave out the first and last second, where the filters start and stop.

The differences are small, so I quote them in microvolts, µV. One microvolt is a thousandth of a millivolt. The picture at the top of the page scores each stage this way.

Watch the error fall at every stage, from 244 µV raw to 111, 86 and 82 µV. Then notice the last two numbers: 82 µV is left, and the filters move the clean ECG by 80 µV on their own. Only 18 µV of noise is left, so most of the 82 µV is the filters’ own change to the ECG.

The high-pass, run both ways

Baseline wander sits below 0.5 Hz. Apart from their average, the beats have almost nothing there, so a high-pass at 0.5 Hz removes the wander; the average comes back later. I use the second-order Butterworth of Analog prototype filters (20.1), turned into a high-pass as in Frequency transformations (20.4).

A high-pass run in real time bends the phase of the slow parts of each beat, so the P and T waves shift and change shape. A stored ECG does not need real time. So I use the forward–backward filtering of “Run it forwards, then backwards” in Choosing FIR or IIR (20.6): no delay, no phase, and the gain squared. On the clean ECG, one pass changes the shape by 30 µV, both ways by 1.9 µV.

The squared gain moves the band edge. One pass is −3 dB at 0.50 Hz. The pair is −3 dB where one pass is −1.5 dB, since

That happens a little higher, at 0.62 Hz. So a forward–backward filter’s band edge is not the one you designed: measure it on $\lvert H\rvert^2$.

The breathing wave at 0.25 Hz keeps $\lvert H\rvert^2=0.059$ of its size, so 0.018 mV of its 0.3 mV is left. The 0.05 Hz drift keeps a ten-thousandth.

The notch, run both ways

For the hum I use a notch at 50 Hz, as in 17.4. There the notch’s width came from its pole radius $r$. SciPy’s iirnotch asks instead for $Q$, the centre frequency divided by the −3 dB width. It is the same idea as the $Q$ of Audio equalisers and biquads (20.5): a larger $Q$, a narrower band (this $Q$ is a number, not the Q wave).

With $Q=30$ at 50 Hz, one pass is $50/30=1.67$ Hz wide. Its poles sit at $r=0.986$, and 17.4’s rule $(1-r)f_s/\pi$ gives nearly the same width, 1.65 Hz. Run both ways, the −3 dB band widens to 48.72 Hz to 51.31 Hz. The hum is left with less than 1 µV.

The low-pass, and what it cannot remove

Here the muscle noise is spread evenly over every frequency up to $f_s/2=180$ Hz. The beats need far less: even the sharp QRS has almost all of its energy below 40 Hz. So a low-pass at 40 Hz, a fourth-order Butterworth, run both ways, removes the muscle noise above the beats’ band. As with the high-pass, the pair’s −3 dB point moves, here down to 36.09 Hz.

The muscle noise inside the beats’ own band stays. No filter can tell it apart from the beats, because it sits at the same frequencies. Of its 30 µV, 14 µV is left.

What each filter costs

Now to the 80 µV change that the filters make to the clean ECG. Almost all of it comes from the high-pass, and it is not rounding. The clean ECG’s waves are mostly above the line, so over seconds 1 to 9 it averages 0.079 mV. A high-pass passes nothing at 0 Hz, so it removes that average: its output averages nearly 0.

So the high-pass moves the whole trace down by about that much. The flat stretch between beats, at 0 mV in the clean ECG, sits at −0.078 mV after cleaning. The shape barely changes: about its own average, the change is 1.9 µV.

The low-pass does round the beats, but only a little. On its own it changes the clean ECG by 2.9 µV and lowers the R peaks from 0.99 mV to 0.97 mV. The notch changes it by less than 1 µV.

That is what the first picture’s last stage says. Its R peaks of 0.90 mV are the clean 0.99 mV, moved down by the 0.08 mV shift and rounded by the low-pass. Measured from the flat stretch before each beat, an R peak still stands 0.97 mV tall.

The errors add up as powers do. In “A longer average: less noise, more delay” of Simple smoothing filters (18.2), the powers of unrelated noises added. Here each part of the difference is nearly unrelated to the others, so the error is close to the square root of the sum of their squares:

Every entry is an RMS in µV over seconds 1 to 9. Each filter removes one column, and the high-pass adds the ECG change. In the last row the 80 µV dominates: $\sqrt{80^2+13^2+1^2+14^2}=82$.

Square, smooth, threshold: a beat

The clean trace is nice to look at, but a monitor wants numbers: when did each beat happen, and how fast is the heart going? For that it must find each QRS. Height alone is not enough: in this raw record, one T wave reaches 0.82 mV, above the lowest R peak, 0.75 mV.

In 1985 J. Pan and W. J. Tompkins published a detector that is still widely used and taught. The Pan–Tompkins detector turns each QRS into one clear bump, then catches the bumps with a threshold. It runs forwards only, so it works in real time. Here are its stages, in a simplified version.

Band-pass

The QRS is the fastest part of a beat: in this model, 55 % of its energy lies between 5 and 15 Hz, against 3 % of the T wave’s. The baseline wander is far slower and the hum faster. So the first stage is a band-pass from 5 to 15 Hz, run forwards only. It starts from a second-order Butterworth prototype, and 20.4’s substitution splits each of its two poles in two.

Its gain is 1.00 at 10 Hz and 0.707 at 5 and 15 Hz. At 1 Hz it is 0.018, and at 50 Hz 0.038. After it, each QRS swings by up to 0.54 mV, while the T waves stay under 0.09 mV.

Slope

A QRS is steep and a T wave is gentle, so the next stage takes the slope, as in “A slope from samples” of Special FIR filters (19.4). Here the input $x[n]$ is the band-passed trace. Pan and Tompkins use five taps:

On a straight line this gives exactly the line’s slope per second: the outer pair spans 4 samples, the inner pair 2, and $4+2\cdot2=8$. The result is the slope 2 samples back, at the taps’ centre.

Square

Next, square each slope value. That removes the sign, so a falling edge counts as much as a rising one. It also makes big values much bigger than small ones. Outside the QRS, the squared slope never passes 3.2 % of its largest value.

Moving-window integration

The squared slope of one QRS has several sharp spikes, one for each steep edge. The next stage averages them into one bump. It is the moving average of “A longer average: less noise, more delay” in 18.2, over 54 samples, which is 150 ms. Pan and Tompkins call it moving-window integration.

As 18.2 promised, the average costs delay. Together with the band-pass and the slope, each bump peaks about 137 ms after its R peak.

Threshold

The last stage is the decision of “Where to set the threshold” in Matched filters and detection (26.1): call it a beat when the bump passes a threshold. I set it at 30 % of the integrated trace’s largest value in the first 2 s, so the detector learns the record’s scale from its start.

A bump’s top is not perfectly smooth: its small ripples give 28 local peaks above the threshold for 12 bumps. So each beat opens a 200 ms window, its refractory time. A higher peak inside that window replaces the beat, and a lower one is ignored, so the detector keeps the highest peak in each window. Heart muscle cannot fire again that soon after a beat, so no real beat is lost.

The bump comes late, so it does not mark the R peak itself. To place the R peak, I look back over the 150 ms before each bump’s peak and take the largest sample of the cleaned ECG from the first instrument. That cleaned trace was filtered both ways, so this last step works on the stored record.

Watch how small the T waves are in the first panel, then watch the bumps rise in the third, each well above the 30 % line. Notice that every triangle lines up with a true tick: the largest miss is 1.33 ms. That is less than half a sample, 1.39 ms, because each found R peak is the sample nearest the true R time.

Heart rate

The heart rate counts beats per minute. With $T_\text{RR}$ in seconds,

For a whole record I use the mean interval. The 12 detected beats give 11 intervals with a mean of 0.796 s, so $\mathrm{HR}=75.4$ beats per minute. The true beat times give the same rate to that rounding.

The real detector does more than this page’s version. It keeps two thresholds that adapt as it goes, one for signal peaks and one for noise peaks. When no beat has come for too long, it searches back with a lower threshold. And its filters use whole-number coefficients, so that it ran in real time on the small processors of 1985.

Worked example

Heart rate from the beats. The detector puts the first R peak at sample 144 and the last at sample 3296. Their times are $144/360=0.400$ s and $3296/360=9.156$ s.

The mean of the 11 intervals is the total span divided by 11, since the times in between cancel: $(9.156-0.400)/11=0.796$ s. Then $\mathrm{HR}=60/0.796=75.4$ beats per minute.

A notch’s width, both ways. Take the notch with $Q=30$ at 50 Hz. One pass is $50/30=1.67$ Hz wide at −3 dB, from 49.17 Hz to 50.84 Hz.

Run forwards and backwards, the gain is squared. Its −3 dB edges lie where one pass is −1.5 dB, which is farther from 50 Hz. Root-finding on $\lvert H\rvert^2$ gives 48.72 Hz and 51.31 Hz, a band 2.59 Hz wide.

The high-pass’s shift. The clean ECG averages 0.079 mV over seconds 1 to 9. The high-pass removes that average, so it alone puts about 80 µV into the final error, whatever the other filters do.

Where you’ll meet this

Bedside monitors in hospitals clean the ECG and count beats with stages like these. So do ambulatory monitors that a patient wears at home for a day or more. Fitness bands and smartwatches find beats too, from an ECG or from light shone through the skin.

Automated defibrillators analyse the heart’s rhythm before they advise a shock, and that analysis also starts from a filtered ECG. In the Americas the mains runs at 60 Hz, so the notch moves to 60 Hz. Brain signals are quieter and their rhythms slower: EEG and rhythms (31.2) takes them up next.

I left out the adaptive thresholds and the search back, the classification of irregular rhythms, and the geometry of the 12-lead ECG. For more, see J. Pan and W. J. Tompkins, “A real-time QRS detection algorithm” (IEEE Trans. Biomed. Eng., 1985), and Sörnmo and Laguna, Bioelectrical Signal Processing in Cardiac and Neurological Applications (2005), chapter 7. The synthetic ECG follows the idea of P. E. McSharry and others, “A dynamical model for generating synthetic electrocardiogram signals” (IEEE Trans. Biomed. Eng., 2003).

Reference card