# Low pass filter output showing instability

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


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?

• 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$. – robert bristow-johnson Oct 2 '19 at 3:36
• 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$. – 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 );

• 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. – brlauwer324 Oct 2 '19 at 4:50
• Did you look at stackoverflow.com/questions/50588879/…? – Maxtron Oct 2 '19 at 5:17
• 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. – brlauwer324 Oct 2 '19 at 5:39