IIR filter from diff. equation

Hey,

I got a rather silly problem. I got the difference equation of a transfer function, and I need to create a filter from this. My first thought was to take the a's and b's directly from it, but if I do that, and use

 bode(b,a);

I get what looks like a highpass filter, from a difference equation that is supposed to be of a lowpass. The difference equation is

 y(n) = 2*y(nT-T) - y(nT-2T) + x(nT) - 2x(nT-6T) + x(nT-12T)

I'm not gonna tell you what I took a and b to be, since I think i got it wrong :-)

Answer 1

[Wayne King]

Perhaps you are not representing the system correctly in bode()?

 A = [1 -2 1];
 B = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
 fvtool(B,A)

[Jakob Sørensen]

fvtool(B,A) shows a lowpass filter (as I hoped), but when i try bode(B,A), it shows a highpass. Is it supposed to be bode(A,B), and if so, why?

Edit: And how come A = [1 -2 1];?

Edit to the edit: Nevermind about the part about A, I just needed another cup of coffee to see that one for myself :p

[Wayne King]

No, it's not bode(A,B), the problem is that bode() is operating on the premise that you have positive powers of the variable (not negative as you have), so it's interpreting:

 A = [1 -2 1];

as z^2-2z^1+1 for example.

Answer 2

[Wayne King]

This is a lowpass filter:

 aLowpass = [1 -2 1];
 bLowpass = [1 0 0 0 0 0 -2 0 0 0 0 0 1];

As you can see with:

fvtool(bLowpass,aLowpass,'Fs',100)

but it is not a very good one. And the highpass filter is not particular good either.

Since you have the Signal Processing Toolbox, why not design your filters with that software?

For example, say you want a lowpass Butterworth (IIR) filter for data sampled at 100 Hz and you want to lowpass everything below 10 Hz. I'll make the attenuation in the stopband 40 dB and the passband ripple 0.5 dB.

D = fdesign.lowpass('Fp,Fst,Ap,Ast',10,15,0.5,40,100);
filtobj = design(D,'butter');
fvtool(filtobj)

Now you can filter your data with:

lowpass_ecg = filter(filtobj,ecg_chan);

By the way, you are supposed to accept answers when people have answered your question.

[Jakob Sørensen]

I see I forgot to mention an important thing, the fact that the assignment I'm working on is to implement an existing algorithm. This algorithm is made by a guy called Tompkins and some friends, back in 1985, and is made for realtime QRS detection. This is why I'm trying to make this hopeless thing work :-)

Answer 3

[Jakob Sørensen]

Oka, I have been looking at it, over and over again, and I think I need some more help...

The algorithm I'm trying to implement is from an article about realtime QRS (ECG spikes) detection. In the article, two filters (a lowpass and a highpass) is given as the following difference equations:

Lowpass with cutoff freq. of around 11 Hz and order 2:

 y(nT) = 2*y(n-T) - y(nT-2T) + x(nT) - 2x(nT-6T) + x(nT-12T)

And a highpass with cutoff freq of 5 Hz and order ?:

 y(nT) = 32*x(nT-16T) - [y(nT-t) + x(nT) - x(nT - 32T)]

This, I've implemented like this:

% Lowpass filtering
aLowpass = [1 -2 1];
bLowpass = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
lowpass_ecg = filter(bLowpass, aLowpass, ecg_chan);

% Highpass filtering
aHighpass = [1 1];
bHighpass = [-1 zeros(1,15) 32 zeros(1,14) 1];
bandpass_ecg = filter(bHighpass, aHighpass, lowpass_ecg);

% Use fvtool to illustrate filters
fvtool(bLowpass, aLowpass, bHighpass, aHighpass);

But when I look at the result (fvtools in the last line), it doesn't look right at all. The sampling frequency is 100 Hz, if that is relevant somehow.

Now that was a long question, I know, but I would really appriciate it. Thanks!