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

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.