# Zero-frequency response of an IIR transfer function

I have a question about the zero-frequency (DC) response of a IIR filter. The coefficients of the 3rd order IIR filter are as follows,

% polynomial coefficients for the numerator of the IIR transfer function
b = [0.0066324, -0.0130292, 0.0063988];

% polynomial coefficients for the denominator of the IIR transfer function
a = [1, -1.4158855, 0.4158913];


Note that the roots of a are all less than unity, but one is quite close to the unity. I am trying to understand why MATLAB's filter() function produces the following results:

rng default;
x = rand(15,1);
y = filter(b,a,x);
mean(y)/mean(x)

>> 0.0012


In the last statement, I computed the ratio of the input and output mean, MATLAB returns 0.0012.

However, from the zero-frequency of the transfer function $H(z)=0.3448$ (i.e., sum(b)/sum(a)=0.3448), I have expected that mean(y)/mean(x)=0.3448. Checking this,

[H,W] = freqz(b,a);
H(1)

>> 0.3448


Have I misunderstood something here about MATLAB's filter() function, or other misconception of the zero-frequency response of a transfer function? Thanks to anyone who is able to provide an insight to this problem.

Because of the pole near the unit circle the impulse response doesn't really settle until after tens of thousands of samples. In comparison, your input of 15 samples is very short. Try making it longer:

x = rand(15E6,1);
y = filter(b,a,x);
mean(y)/mean(x)

>> 0.34252


Also, the filter starts in a state that corresponds to an all-zeros input history. If you precondition it with a run of the actual input, the result will converge faster:

x = rand(15E6,1);
y = filter(b,a,x);
mean(y(5E6:length(y)))/mean(x(5E6:length(x)))

>> 0.34483


Your basic assumption is correct: if you convolve two sequences

$$y[n]=(x\star h)[n]\tag{1}$$

then their $\mathcal{Z}$-transforms are related by

$$Y(z)=X(z)H(z)\tag{2}$$

from which

$$Y(1)=X(1)H(1)\tag{3}$$

follows. In the case of a recursive filter, $H(z)=B(z)/A(z)$ and

$$H(1)=\frac{B(1)}{A(1)}=\frac{\sum_nb_n}{\sum_na_n}$$

where $a_n$ and $b_n$ are the denominator and numerator coefficients, respectively.

The problem with your implementation is that the function filter() does not compute the complete convolution of the input signal with the filter's impulse response, but it only computes as many output samples as there are input samples. Note that the length of the actual convolution result is always longer than the input sequence. And in your case it is very long, actually infinitely long (because you have an IIR filter).

Olli Niemitalo's solution is OK, but it does not address all details of the problem, because making the input signal longer does not completely solve the problem. It only decreases the relative contribution of the transient due to the filter's memory. The transient response of the filter after the input signal is switched off is still not taken into account.

Equation (3) must of course hold for any input signal, also for very short signals. So what you need to do is append zeros to the input signal such that you also get the transient response after the input signal is switched off. Theoretically you would need to append infinitely many zeros because your filter's memory is (theoretically) infinitely long. And since one of the filter's poles is very close to the unit circle, your filter's memory is especially large, even for practical purposes.

In Matlab/Octave you could do something like this:

b = [0.0066324, -0.0130292, 0.0063988];
a = [1, -1.4158855, 0.4158913];
x = rand(15,1);
x0 = [x;zeros(1e6,1)];
y0 = filter(b,a,x0);
[sum(y0)/sum(x0),sum(b)/sum(a)]    % should be (almost) the same


Part A: Bilinear Transformation Method

1. Design lowpass IIR filter with the following specifications:

Filter order = 2, Butterworth type

Cut-off frequency=800 Hz

Sampling rate =8000 Hz

Design using the bilinear z-transform design method

Print the lowpass IIR filter coefficients and plot the frequency responses using MATLAB.

MATLAB>>freqz(bLP,aLP,512,8000); axis([0 4000 –40 1]); % sampling rate=8000 Hz