1
$\begingroup$

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);

figure;
subplot(2,1,1)
plot(lowPassResult);
title('Lowpass')
subplot(2,1,2)
plot(sig_108);
title('Original')

function y = lowPassFilter(x)
    %y(n)=2y(n-1)-y(n-2)+x(n)-2x(n-6)+x(n-12)
    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);
    end
end

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?

Improve this question
$\endgroup$
2
  • 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
    Commented Oct 2, 2019 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
    Commented Oct 2, 2019 at 4:29

2 Answers 2

Reset to default
1
$\begingroup$

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 );
Improve this answer
$\endgroup$
3
  • $\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
    Commented Oct 2, 2019 at 4:50
  • $\begingroup$ Did you look at stackoverflow.com/questions/50588879/…? $\endgroup$
    – Maxtron
    Commented Oct 2, 2019 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
    Commented Oct 2, 2019 at 5:39
0
$\begingroup$

This is an old question, but I recently ran into this problem myself and spent a lot of time trying to figure it out. @brlauwer324 provides the answer in the comment: x(1:12) should be initialized to zero. Here's the MATLAB code:

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

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

figure;
subplot(2,1,1)
plot(lowPassResult);
title('Lowpass')
subplot(2,1,2)
plot(sig_108);
title('Original')

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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