A question regarding MATLAB's 'invfreqz()' command

Started by Rick Lyons 5 years ago6 replieslatest reply 5 years ago273 views
Hi. When I execute the following MATLAB code:
  b = [1 0 0 0 -1]
  a = [1, -1]
  [h, w] = freqz(b, a, 64);
  b_Order = length(b) -1; a_Order = length(a) -1;
  [bb, aa] = invfreqz(h, w, b_Order, a_Order)

I expect to produce

  bb =

     1  0  0  0  -1
  aa =

       1   -1

But instead of the above 'bb' and 'aa' sequences my code produces:

  bb =

     NaN   NaN   NaN   NaN   NaN
  aa =

       1   NaN

Can someone here tell me what's going on here? Thanks.
[ - ]
Reply by anamariatomeOctober 23, 2019

if you change a=[1,-0.9] everything works as you expect. With your example the array h has NaN in the first position  I guess that is the problem. 

[ - ]
Reply by Rick LyonsOctober 23, 2019

Hi anamariatome.
Thanks for pointing out that h(1) = NaN condition.
(I hadn't noticed that.) To correct for that h(1) = NaN condition, I replaced my '[h,w] = freqz(b, a, 1024);' command with:

  [h,w] = freqz(b, a, 1024);
  h(1) = sum(ImpResp);

But with that change in my code I still obtained wildly incorrect values for variables 'bb' and 'aa'.

Following your suggestion, going back to my original code, when I redefined the denominator 'a' coefficients as:

  a = [1, -.9999]

to move the pole very slightly inside the unit circle I obtained correct 'bb' and 'aa' results! So it's not the h(1) = NaN condition causing me problems, it's the pole at z = 1 causing a divide by zero condition inside MATLAB's invfreqz() command. Moving that pole to very slightly inside the unit circle is the thing to do.

I think this all means that MATLAB's invfreqz() command has a "bug."

Thanks again.

[ - ]
Reply by kazOctober 23, 2019

Hi Rick,

I would think the dc gain of your filter is to blame.

iir dc gain = sum(b)/sum(a) = 0/0 = NAN because of divide by zero... if I am right...

[ - ]
Reply by omersayliOctober 23, 2019

Hi Mr Lyons,

First value of h is NaN as I have checked. Maybe that is the reason?

Edit: (As noted by anamariatome also, that is the problem indeed, pole at z=1(f=0) causing NaN value) I tried with other a,b values

[ - ]
Reply by Rick LyonsOctober 23, 2019

Hi omersayli. Please call me "Rick." My friends do.

Thanks for your thoughts, and you are correct. Check out my Oct. 22 Reply to anamariatome.

[ - ]
Reply by omersayliOctober 23, 2019

You are right Rick , 

  b = [1 0 0 0 -1]

  a = [1 -0.9]

  freqz_points = 512

  [h, w] = freqz(b, a, freqz_points);

  [imp,t] = impz(b,a);

  h(1) = sum(imp)

  b_Order = length(b)-1; a_Order = length(a)-1;

  [bb, aa] = invfreqz(h, w, b_Order, a_Order)

trying this, it works (I substituted h(1) with impz result as you suggested also). When trying   a = [1 -1], you get different results for different "freqz_points" parameter! That is really weird

EDIT: Interestingly, When making h(1) = 0 , the 'invfreqz' function works quite good, at least for the a, b arrays I have tried.