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?