I am fairly new to DSP, and I'm trying to implement the Pan Tompkins algorithm for QRS detection of ECG signals in MATLAB. The first stage of the algorithm consists of a second-order low-pass filter. The difference equation of the filter is

$$ y(n)=2y(n-1)-y(n-2)+x(n)-2x(n-6)+x(n-12) $$

Here is the MATLAB code I have written.

load 108m.mat; sig_108 = val(1,:); clear val;
lowPassResult = lowPassFilter(sig_108);


function y = lowPassFilter(x)
    y = zeros(1,21600);
    y(12) = 0;
    y(11) = 0;
    for n = 13:length(x)
        y(n) = 2*y(n-1)-y(n-2)+x(n)-2*x(n-6)+x(n-12);

Output of MATLAB code

I am using ECG record 108 from the MIT-BIH arrhythmia database available at PhysioBank ATM. One problem I am having is that the output of the lowpass filter appears to grow without bound towards $+\infty$ as can be seen in the figure above. Where have I gone wrong in my code that could be causing this issue?

  • 1
    $\begingroup$ well the coefficients for the feedback paths correspond to two poles right on the unit circle at $z=1$. to be stable they need to be inside the unit circle meaning $|z|<1$. $\endgroup$ – robert bristow-johnson Oct 2 '19 at 3:36
  • $\begingroup$ Yes, I actually did look into this by finding the roots of the characteristic polynomial in the time domain and found that the equation has a double root at $\lambda = 1$. $\endgroup$ – brlauwer324 Oct 2 '19 at 4:29

I am going to explain Robert's answer below.

Apply Z-transform to both sides of the equation. You get the following: \begin{align} Y(z) [1 - 2z^{-1} + z^{-2}] = X(z)[1 -2z^{-6} + z^{-12}] \end{align}

The transfer function is given as \begin{align} H(z) = \frac{Y(z)}{X(z)} = \frac{1 - 2z^{-1} + z^{}-2}{1 - 2z^{-6} + z^{-12}} = \frac{(1-z^{-6})^{2}}{(1-z^{-1})^{2}} \end{align}

The denominator represents the poles of the filter $H(z)$. In this case, there exists two poles, both of which lie on the unit circle at $z=1$. For a filter to be stable, all the poles should lie inside the unit circle, i.e., the roots of the denominator polynomial of the transfer function $H(z)$ should lie inside the unit circle.

Perhaps, you can also try a different MATLAB implementation of the above filter:

a1 = [ 1 −2 1 ];
b1 = [ 1 0 0 0 0 0 −2 0 0 0 0 0 1 ];
lp_sig = filter( b1, a1, x );
  • $\begingroup$ I tried the code you provided and the correct output was generated. However, I can't use the built-in filter function, but now I know that my error lies in my code. I'm guessing it likely is due to the initial conditions of $y(n)$ that I have set. $\endgroup$ – brlauwer324 Oct 2 '19 at 4:50
  • $\begingroup$ Did you look at stackoverflow.com/questions/50588879/…? $\endgroup$ – Maxtron Oct 2 '19 at 5:17
  • 1
    $\begingroup$ I did not, but I found that $x(n)$ is relaxed for $n<0$, so $x(n)=0$ for all $n<0$. Thus, I need to re-index the input by adding 12 zeros at the beginning of $x(n)$ to prevent any errors regarding negative indices. $\endgroup$ – brlauwer324 Oct 2 '19 at 5:39

Your Answer

