From help-request at octave dot org Wed Jan 12 12:25:46 2005
Subject: Re: freqz inconsistency ?
From: "Pascal A. Dupuis" <Pascal dot Dupuis at esat dot kuleuven dot ac dot be>
To: help at octave dot org
Date: Wed, 12 Jan 2005 10:19:16 -0600

Some more tests:

pdi = 0.90000
%# try system as (1-pdi)/(1 - pdi z^-1)

[Zz1, Wz1]=freqz((1-pdi), [1 -pdi], 128);
[Zz2, Wz2]=freqz((1-pdi), [1 -pdi], Wz1);

[(1-pdi)/(1-pdi*exp(-j*Wz1(2))) Zz1(2)]
ans =

  0.95115 - 0.20951i  0.95115 - 0.20951i
OKAY (MatLab agrees whith this )

[(1-pdi)/(1-pdi*exp(-j*Wz2(2))) Zz2(2)]
ans =

  0.95115 - 0.20951i  0.94572 - 0.23279i
WRONG.

%# This can be inverted as:

We1=-imag(log(1/pdi*(1-(1-pdi)./Zz1)));
We2=-imag(log(1/pdi*(1-(1-pdi)./Zz2)));

[Wz1(1:5); We1(1:5)]
ans =

   0.000000   0.024544   0.049087   0.073631   0.098175
  -0.000000   0.024544   0.049087   0.073631   0.098175

aggreement is perfect

[Wz2(1:5); We2(1:5)]
ans =

  0.00000   0.02454   0.04909   0.07363   0.09817
  -0.00000   0.02725   0.05439   0.08131   0.10791

It fails !

The frequency increments are:
octave>pi/(Wz1(2)-Wz1(1))
ans = 128
octave> pi/(We1(2)-We2(1))
ans = 128.00
octave> pi/(We2(2)-We2(1))
ans = 115.28

So the frequency seems incorrect. Looking into the code, there is
nothing wrong with frequency, BUT the computation is performed with
fft in case of scalar n, and with polyval in case of vector n. 

Could someone please sort things out ? One of the problem is linked
with polynomial ordering: MatLab docs states that the polynomial is
written a(1) + a(2)z^-1 + ..., does octave follows the same convention ?

The other one is linked with the discrepancy in the fft<->polyval
computation, I guess the right way is to compute :
polyval(a(na:-1:1), exp(-j*w))

TIA

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



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