Manually discretise a Second Order low pass filter

I need to manually create a second order digital filter for a set of data. I have carried out FFT on the data to identify the 4 sin waves the data is made up of and have identified the highest frequency sin wave as the one I wish to filter out.

So far I have understood that a second order filter can be represented as the transfer function:

y= wn^2/s^2+(2*zeta*wn*s)+wn^2

where wn = natural frequency of highest freq sin curve

zeta=damping ratio s=differential operator

I believe this equation can be used to identify the coefficients that can be input into Matlab's own bilinear function, although I also need to carry out manual bilinear transformation in order to provide a comparisons.

I have also read that to discretise, the transfer function should be in terms of z^-1 and the relationship between s and z^-1 is:

s=2(1-z^-1)/Ts(1+Z^-1)

At this stage, I have tried to substitute the equation for s back into my original second order equation and rearrange to get coefficients into a format that can be compared to the discrete time transfer function.

This is however proving difficult and I am not 100% on the methodology. Any advice, pointers of where to find additional information would be appreciated.

Thanks in advance,

Mike

Answer 1

[Star Strider]

It’s just as straightforward as you describe it. It’s tedious algebra to do it manually, but if you have the Symbolic Math Toolbox, it’s easy:

syms wn s zeta Ts z
y= wn^2/s^2+(2*zeta*wn*s)+wn^2;
s2z=2*(1-1/z)/(Ts*(1+1/z));     % Bilinear Transform, ‘s’ To ‘z’
y = subs(y, s, s2z);            % Convert Continuous To Discrete
[yn,yd] = numden(y);            % Numerator, Denominator
yn = collect(yn, z)             % Numerator Polynomial In ‘z’
yd = collect(yd, z)             % Denominator Polynomial In ‘z’
yncf = coeffs(yn, z)            % Numerator Coefficients
ydcf = coeffs(yd, z)            % Denominator Coefficients

Answer 2

[Mike Scott]

Great thanks. Unfortunately I don't have the Toolbox, but I'll have a crack at it manually and see how I get on.

Thanks again

[Star Strider]

My pleasure!

Since you don’t have the toolbox, I’ll make life a bit easier for you:

yn =

wn*(wn*Ts^3 + 4*wn*Ts + 16*zeta)*z^3 - wn*(- 3*wn*Ts^3 + 4*wn*Ts + 48*zeta)*z^2 + wn*(3*wn*Ts^3 - 4*wn*Ts + 48*zeta)*z + wn*(wn*Ts^3 + 4*wn*Ts - 16*zeta)

yd =

4*Ts*z^3 - 4*Ts*z^2 - 4*Ts*z + 4*Ts

yncf =

[ wn*(wn*Ts^3 + 4*wn*Ts - 16*zeta), wn*(3*wn*Ts^3 - 4*wn*Ts + 48*zeta), -wn*(- 3*wn*Ts^3 + 4*wn*Ts + 48*zeta), wn*(wn*Ts^3 + 4*wn*Ts + 16*zeta)]

ydcf =

[ 4*Ts, -4*Ts, -4*Ts, 4*Ts]

(Command Window output from my code.)

[Chiao-Ting Li]

Hi,

I recently have a similar question and stumbled into this post.

The original question posed the transfer function *without* brackets in the denominator (by an innecent mistake I think).

The correct formula for a second order filter should be:

y= wn^2/ ( s^2+(2*zeta*wn*s)+wn^2 )

And, therefore, the answer provided by @Star should be modified to be the following (the methodology is correct though):

yn =
Ts^2*wn^2*z^2 + 2*Ts^2*wn^2*z + Ts^2*wn^2
yd =
(Ts^2*wn^2 + 4*zeta*Ts*wn + 4)*z^2 + (2*Ts^2*wn^2 - 8)*z + Ts^2*wn^2 - 4*zeta*Ts*wn + 4
yncf =
[ Ts^2*wn^2, 2*Ts^2*wn^2, Ts^2*wn^2]
ydcf =
[ Ts^2*wn^2 - 4*zeta*Ts*wn + 4, 2*Ts^2*wn^2 - 8, Ts^2*wn^2 + 4*zeta*Ts*wn + 4]

thought I'd give an update on this post, in case someone just copy-and-paste the origial answer and get wiered results like I did.

Regards,

Chiao-Ting

[Piotr]

Hello Chiao-Ting,

thank you for the update. I have a question as I do not fully understand the units behind each of the parameters used in the main formulae.

In my problem I have information provided for the damped frequency which is in Hz and damping factor which is in Mneper/s. What units do you use? Is Ts a sampling rate [s]?

BR, Piotr