IIR filter from diff. equation

9 次查看（过去 30 天）
Jakob Sørensen
Jakob Sørensen 2012-9-28
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 :-)

请先登录，再回答此问题。

采纳的回答

Wayne King
Wayne King 2012-9-28
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)
  2 个评论
Jakob Sørensen
Jakob Sørensen 2012-9-28
编辑：Jakob Sørensen 2012-9-28
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
Wayne King 2012-9-28
编辑：Wayne King 2012-9-28
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.

请先登录，再进行评论。

更多回答（2 个）

Wayne King
Wayne King 2012-10-1
编辑：Wayne King 2012-10-1
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.
  1 个评论
Jakob Sørensen
Jakob Sørensen 2012-10-1
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 :-)

请先登录，再进行评论。

Jakob Sørensen
Jakob Sørensen 2012-10-1
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!

类别

在 Help Center 和 File Exchange 中查找有关 Frequency Transformations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by