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

Question

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}}$

Answer 1

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.

Answer 2

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.