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