Re: Using invfreqz() without knowing the phase response?

[Sergei Steshenko]

----- Original Message -----
> From: "OCuanachain, Oisin (Oisin)" <address@hidden>
> To: "address@hidden" <address@hidden>
> Cc: 
> Sent: Tuesday, July 17, 2012 4:39 PM
> Subject: Using invfreqz() without knowing the phase response?
> 
> Hi, 
> 
> I would like to use invfreqz() to create a digital filter [b(z),a(z)] as a 
> model 
> approximating some empirical data. The problem is that I only have the 
> amplitude 
> response of the system not the phase response. This shouldn't really be a 
> problem as I am only interested in the amplitude response however to apply 
> invfreqz() I need a phase response. Eg
> Amplitude response = A(w)
> Phase response     = phi(w)
> Then 
> H                  = A(w).*exp(-1*j*phi(w))
> Then [b(Z), a(z)] = invfreqz(H,w);
> 
> So what I really want to know is, what is an 'easy' dummy phase response 
> I can use for phi(w) which will maximize the probability of invfreqz 
> generating 
> digital coefficients which will match my desired amplitude response well ?  
> Or 
> to put it another way what phi(w) can I use to ensure that invfreqz() 
> doesn't run into problems. 
> 
> I realise I could probably work this out myself by reading the source code 
> for 
> invfreqz() to see exactly what algorithm it uses, but if someone else already 
> knows the answer it would save me the trouble :)
> 
> Thanks, 
> 
> Oisín. 
> 

The "standard" way is to simply assume that you need a fit for a minimal phase 
system.

For a minimal phase system there is one to one correspondence between magnitude 
and phase responses, i.e. one can be calculated from the other. See 
https://en.wikipedia.org/wiki/Minimum_phase -> 
https://en.wikipedia.org/wiki/Minimum_phase#Relationship_of_magnitude_response_to_phase_response
 .

The formula is:

min_phase = unwrap(-imag(hilbert(log(abs([frequency_response, 
fliplr(frequency_response(2:end))])))), pi); # thanks to Ben Abbot for 'fliplr' 
part

- in this case 'hilbert' is Octave/Matlabe function;

imag(hilbert(foo))

is true Hilbert transform of 'foo'.

Then, when you have your phase response calculated, use it with magnitude 
response in order to create teh complex spectrum to be fit.

...

You have to realize that for various reasons implementation of Hilbert 
transform is approximate. There was a recent thread with "Hilbert transform" 
subject in this list, you may find it interesting.

Regards,
  Sergei.