2nd Order transfer function to difference equations

I'm wondering if someone could check to see if my conversion of a standard second order transfer function to a difference equation is correct, and maybe also help with doing a computer implementation. Starting Equation: $$\frac{Y(s)}{R(s)}=\frac{\omega_n^2}{s^2 + 2\zeta\omega_ns + \omega_n^2}$$

Using the backwards-difference equation, $$s = \frac{1-z^{-1}}{T}$$

Result: $$a = 2(1-\zeta\omega_nT)$$ $$b = 2\zeta\omega_nT + \omega_n^2T^2 + 1$$ $$y(k) = \frac{r(k-2) - ar(k-1) + br(k)}{\omega_n^2T}$$

Now, how would someone implement in an algorithm to filter data (since this is just a 2nd order low pass filter)?

% Coded in MATLAB
time = 1:1024;
wn = 0.5;
zwn = 1.414;
T = 1;
a = 2*zwn*T + (T^2)*wn^2 + 1;
b = T*wn^2;
c = 2*(1-zwn*T);
output = 10*sin(0.05*time + 10) + 4*sin(1*time+180);
y(1) = output(1);
y(2) = output(2);

for i=3:length(time)
y(i) = (1/b)*(a*output(i) - output(i-2) + c*output(i-1));
end

I would expect this to filter the output's 1 Hz noise, but all this does is reduce the amplitude of the entire graph. Where am I going wrong?

• You might be better off using the bilinear transform - this preserves DC gain: $s\rightarrow \frac{2}{T}\frac{z-1}{z+1}$. What are the $\zeta$ and $\omega _n$ values? – Chu May 30 '17 at 17:12
• Are you really using a value of zeta of 2.828? – Andy aka May 30 '17 at 17:46
• I'll try to use the bilinear transform and see where that takes me. And yes Andy, I am using a zeta of 2.8. The actual values I'm using are not that important to me at this point in time, more so asking if my implementation and difference equation solution are correct. – Collaptic May 30 '17 at 18:37
• Your final equation for $\small y(k)$ is wrong. The equation should be recursive and the terms in $\small r(k-1)$ and $\small r(k-2)$ should be in $\small y(k-1)$ and $\small y(k-2)$. I haven't double checked your coefficients, but my first run through indicates they are wrong. – Chu May 30 '17 at 20:56