Re: freqz inconsistency ?

[David Bateman]

According to Pascal A. Dupuis <address@hidden> (on 01/12/05):
> Hello,
> 
> I'm trying to understand some issues with invfreqz, and found some
> inconsistency in freqz:
> 
> pdi = 0.97;
> [Zz1, Wz1]=freqz((1-pdi), [1 -pdi], 128);
> [Zz2, Wz2]=freqz((1-pdi), [1 -pdi], Wz1);
> 
> norm(Wz1-Wz2)
> ans = 0
> 
> norm(Zz1-Zz2)
> ans = 0.35010

Under Matlab v7.0 sp1, I get

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

ans =

     0

>> norm(Zz1-Zz2)

ans =

   1.7417e-15

>> Zz1a = Zz1; Zz2a = Zz2; save -v6 temp.mat Zz1a Zz2a

Loading them into octave and comparing with the same values calculated by
octave gives an interesting result, seen in the attached files. It seems that
octave is in error in the second form. Looking at the last values in the

octave:25> Zz2a(end), Zz2(end)
ans = 0.01522846 - 0.00018404i
ans = -1.5228e-02 - 1.8974e-04i

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

Regards
David







> 
> as both frequency responses are evaluated at the same points, they
> should give the same result, shouldn't they ?
> 
> I tried cheating a bit with 
> [Zz2, Wz2]=freqz((1-pdi), [1 -pdi], const*Wz1);
> but couldn't find a constant such that both frequency responses are
> the same. Any idea ?
> 
> 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
> -------------------------------------------------------------

-- 
David Bateman                                address@hidden
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

Zz2_vs_Zz2a.eps
Description: PostScript document

Zz1_vs_Zz1a.eps
Description: PostScript document