From help-request at octave dot org Thu Jan 13 09:06:55 2005
Subject: Re: freqz inconsistency ?
From: "Pascal A. Dupuis" <Pascal dot Dupuis at esat dot kuleuven dot ac dot be>
To: David Bateman <David dot Bateman at motorola dot com>
Cc: help at octave dot org
Date: Thu, 13 Jan 2005 02:44:37 -0600


--vkogqOf2sHV7VnPd
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline

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

--vkogqOf2sHV7VnPd
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="freqz.m.patch"

--- freqz.m.orig	2005-01-12 16:41:23.000000000 +0100
+++ freqz.m	2005-01-13 09:40:48.000000000 +0100
 at @ -68,6 +68,10 @@
 ## than returning them.
 ##  at end deftypefn
 
+## 12 JAN 2005: modified by <Pascal dot Dupuis at esat dot kuleuven dot ac dot be>: 
+## -avoid recomputing length(a), length(b) when la, lb are available
+## -fix (inverted the polynomial order) such that 
+## 	[H1, W1]= freqz(B, A); [H2, W2]=freqz(B, A, W1) give the same result 
 ## Author: jwe ???
 
 function [h_r, w_r] = freqz (b, a, n, region, Fs)
 at @ -125,6 +129,7 @@
     endif
     n = length (n);
     extent = 0;
+    q = exp(-j*w);
   elseif (strcmp (region, "whole"))
     w = 2 * pi * (0:n-1) / n;
     extent = n;
 at @ -133,8 +138,8 @@
     extent = 2 * n;
   endif
 
-  if (length (b) == 1)
-    if (length (a) == 1)
+  if (1 == lb )
+    if (1 == la)
       hb = b * ones (1, n);
     else
       hb = b;
 at @ -143,15 +148,15 @@
     hb = fft (postpad (b, extent));
     hb = hb(1:n);
   else
-    hb = polyval (postpad (b, k), exp (j*w));
+    hb = polyval (b(lb:-1:1), q);
   endif
-  if (length (a) == 1)
+  if (1 == la)
     ha = a;
   elseif (extent >= k)
     ha = fft (postpad (a, extent));
     ha = ha(1:n);
   else
-    ha = polyval (postpad (a, k), exp (j*w));
+    ha = polyval (a(la:-1:1), q);
   endif
   h = hb ./ ha;
   w = Fs * w / (2*pi);

--vkogqOf2sHV7VnPd--



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------