# Behaviour of difference equation does not match that of its z-domain transfer function

I have obtained the following z-domain transfer function:

$$\frac{Y(z)}{U(z)}=\frac{3.3641×10^{-7}×z^{6}+1.5584×10^{-5}×z^{5}+6.7263×10^{-5}×z^{4}+ 5.5016×10^{-5}×z^{3}+8.525×10^{-6}×z^{2}+1.2303×10^{-7}×z+9.7492×10^{-20}}{z^{7}-5.5998×z^{6} +13.5548×z^{5}-18.3685×z^{4}+15.0393×z^3-7.4352×z^2+2.0541×z-0.2245}$$

This is a 7-th order Butterworth filter with a corner frequency of 498.79 Hz.
Now, to obtain the difference equation so I can implement it on a microcontroller, I am multiplying the right side by $$\frac{z^{-7}}{z^{-7}}.$$
I thus obtain:
$$\frac{Y(z)}{U(z)}=\frac{3.3641×10^{-7}×z^{-1}+1.5584×10^{-5}×z^{-2}+6.7263×10^{-5}×z^{-3}+ 5.5016×10^{-5}×z^{-4}+8.525×10^{-6}×z^{-5}+1.2303×10^{-7}×z^{-6}+9.7492×10^{-20}×z^{-7}}{1-5.5998×z^{-1} +13.5548×z^{-2}-18.3685×z^{-3}+15.0393×z^{-4}-7.4352×z^{-5}+2.0541×z^{-6}-0.2245×z^{-7}}$$

I have done this because it allows me to work with values from the past, rather than try and read the future (both input and output). After re-arranging, I get the difference equation:
$$y[n]=5.5598×y[n-1]-13.5548×y[n-2]+18.3685×y[n-3]-15.0393×y[n-4]+7.4352×y[n-5]-2.0541×y[n-6]+0.2245×y[n-7]+3.3641×10^{-7}×u[n-1]+1.5584 ×10^{-5}×u[n-2]+6.7263×10^{-5}×u[n-3]+5.5016×10^{-5}×u[n-4]+8.525×10^{-6}×u[n-5]+1.2303×10^{-7}×u[n-6]+9.7492×10^{-20}×u[n-7]$$

I implemented this difference equation, but I am not getting the same results as for its z-domain transfer function. I couldn't identify any mistake in what I did. Could someone please help me identify it? I am 90% sure there is some mistake in my workings somewhere. I tried to do the same procedure for another (simpler) transfer function and the behaviour of the discrete difference equation matched exactly that of the simulated z-domain transfer function for the same input.

UPDATE 1:

The poles of the z-domain are: $$0.8894±0.2805j, 0.7979±0.1995j, 0.747±0.1022j, 0.731$$ and they are all stable.

Here's more proof it's stable (z domain): I first constructed the following continuous transfer function, which I used together with the MATLAB c2d() function to get the z-domain transfer function I mentioned earliler. The method was "impulse" and a sampling frequency of 10 kHz. The continuous form is:

$$\frac{Y(s)}{U(s)}=\frac{2.9696×10^{24}}{s^{7}+1.4084×10^{4}×s^6 +9.9181×10^7×s^5+4.4917×10^{11}×s^4+1.4077×10^{15}×s^3+3.053×10^{18}×s^2+ 4.2582×10^{21}×s+2.9696×10^{24}}$$

• I'm going to answer your question anyway, but it would be good if you edit your question to say what you are seeing that leads you to say "not the same results". First, because I may go back and refine my answer, and second, because anytime you say "not the same" or "doesn't work" you should show what is happening and what you expected to happen. Jan 21 at 16:38
• Are you using single-precision floating-point numbers or double-precision?
– Ben
Jan 21 at 16:39
• Factoring your polynomial as written I see roots at $1.2857 \pm j0.2660,\ 0.9261 \pm j0.5875,\, 0.4620 \pm j0.4647,\, 0.2521$. Was it your intent to have a strongly unstable system? Jan 21 at 16:51
• (Given the degree to which your filter poles are unstable I'm going to hold off on an answer -- are you sure you got the coefficients right?) Jan 21 at 17:00
• @DanielTork: Try to use the coefficients in your question, WITH the number of decimals that you've shown us. Then you'll see that at least those coeffs will give you an unstable system. But it's good news that your actual coefficients result in a stable system. So either there's an error in your transfer function (the one in the question text), or you got a lesson in coefficient sensitivity. Especially for higher order filters (and 7 is already high), the pole locations are extremely sensitive to small errors in the coefficients. Jan 21 at 18:43

If you want a discrete-time Butterworth filter you should design it in the discrete domain. If you use Matlab (or Octave), you can just use the function butter(). I didn't check your continuous-time transfer function, and I don't know how exactly you transformed it to the discrete domain, but the result is definitely unstable. You can see this by computing the roots of your denominator polynomial. Assuming that the denominator coefficients are stored in the vector a, the following Matlab command computes the roots:

roots(a)
ans =

1.28574 + 0.26605i
1.28574 - 0.26605i
0.92606 + 0.58748i
0.92606 - 0.58748i
0.46204 + 0.46475i
0.46204 - 0.46475i
0.25212 + 0.00000i


Four of these roots are outside the unit circle, i.e., have magnitudes greater than $$1$$, which makes the filter unstable (in analogy to poles in the right half-plane for continuous-time systems).

A direct design in the discrete domain is straightforward and more successful (using the cut-off frequency stated in your question, and a sampling frequency of $$10$$ kHz:

[b,a] = butter(7, 498.79/5000);


You can try this in order to see how far or close your coefficients are to the optimal ones. In practice, you should use the pole and zero locations as output arguments, and compute the coefficients of the second-order sections from the poles and zeros, because that's numerically more stable:

[z,p,g] = butter(7, 498.79/5000);
sos = zp2sos(z,p,k);


Now you can use the coefficients in sos to implement the filter as a cascade of second-order sections. This is a numerically much more stable method than computing the transfer function directly.

In sum, there are two problems with your approach. First, you got the wrong coefficients, and second, you were trying to implement the transfer function directly instead of using second-order sections.

UPDATE:

If your filter is actually stable then the discrepancy between what we saw lies in the high sensitivity of the pole locations to coefficients errors. Using the coefficients given in your question (with $$4$$ decimals) does result in an unstable filter.

The reason why your implementation still doesn't work (apart from possible bugs in your code) is the fact that you directly implemented the transfer function as a high-order difference equation instead of splitting it up into second-order sections.

• Thanks for the patience, I appreciate it! Jan 21 at 19:09

Probably the problems that you are seeing stem from trying to do sensible math on $$7^{th}$$-order polynomials without using insanely high numerical precision.

There are some nearly-absolute dos and don'ts to IIR filter design (either digtial or analog), all having to do with keeping the filter order as small as possible.

• Don't use filter sections of more than 2nd-order*. As the filter order increases, the sensitivty of the poles (and zeros) to the coefficients (or resistor & capacitor values) goes up with filter order. A 2nd-order section is the smallest section you can implement when you have resonant poles -- so keep it to that.
• If at all possible, do design your filter as a collection of poles and zeros. Note that all the old-style "big name" filters such as Butterworth, Chebychev, Gaussian, etc., are all available as tables of coefficients, or in these degenerate days of using math packages, there'll be a flag you can set in the call.
• If you're going to design the filter in the Laplace domain then convert to $$z$$, do it as 2nd-order sections. Or if you're just getting it from a math package, ask for the $$z$$-domain version.

If, for some peculiar reason, you must design your filter as a ratio of polynomials of order greater than 2, be painfully aware of the impact this will have on your filter performance. Use whatever extended-precision data types your math package has to offer, and sufficient precision to get a decent answer. While you're doing that, as yourself if it wouldn't be less trouble to design in poles and zeros, or 2nd-order sections.

* This doesn't apply to passive LC or resonant transmission-line filters because you get into a lot less trouble with no amplifying elements.

• Thanks a lot for the patience and everything. If I could accept two answers, yours would be on the list. Jan 21 at 19:11
• You should give it to Tim. Matt doesn't need the rep, he breaks the bank anyway. Jan 21 at 23:15
• Tim is just trying to give good answers, not chasing reputation points. Jan 22 at 18:49