Re: freqz inconsistency ?

[Pascal A. Dupuis]

On Wed, Jan 12, 2005 at 04:35:49PM +0100, David Bateman wrote:
> I don't know enough about freqz to really fix the bug though after a 
> quick look. Expect to note that in the first case extent is a positive
> integer and in the second extent is zero and therefore there are two
> different ways of calculating the frequency response used. That is 
> 
>     hb = fft (postpad (b, extent));
>     hb = hb(1:n);
> 
> or
> 
>     hb = polyval (postpad (b, k), exp (j*w));
> 
> I imagine its a difference at this point causing the problem
> 


I located the problem.

Starting point: H(z)=1/(1-az^-1), encoding: B = 1; A = [1 -a];

Given 0 <= w < pi, the frequency response is H = 1./(1 - a exp(-j*w).
Let's say we compute polyval([1 -a], exp(-j*w)). The problem is that
we'll be computing exp(-j*w) - a: polynomial arguments are used in the
reverse order of 'logical' ordering.

First fix: compute 
z = exp(-j*w);
H = polyval(b(nb:-1:1), z)./polyval(a(na:-1:1), z);

Second fix: multiply A and B by the inverse of the highest power of
z^-n, i.e. consider H(z) = z/(z-a) instead of 1/(1-az^-1), so:
k = max(na, nb);
b = postpad(b, k); %# will give [1 0]
a = postpad(a, k); %# will give [1 -a]
z = exp(j*w);
H = polyval(b, z)./polyval(a, z);

The advantage of the first form is that there is no expansion of the
filters coefficients, so if one of them is of size 1, the resulting
computation may be reduced to avoid dividing two complex numbers.

The actual bug is that, actually, freqz tries to use the second form,
but do not augment polynomials before testing their length. Patch
attached to revert to the first form.

Best regards

Pascal Dupuis

-- 
Dr. ir. Pascal Dupuis
K. U. Leuven, ESAT/ELECTA (formerly ELEN):  http://www.esat.kuleuven.ac.be/
Kasteelpark Arenberg, 10; B-3001 Leuven-Heverlee, Belgium
Tel. +32-16-32 10 21 -- Fax +32-16-32 19 85

freqz.m.patch
Description: Text document