0
$\begingroup$

My problem is this: Given a frequency magnitude response $|\hat{h}(\omega)|$ of a minimum phase system, how can one reliably recover its phase response?

In my case, the magnitude is quite short (~30-50 samples). Here is an example:

frequency magnitude plot of example signal, 31 samples

I don't want to use any approach requiring polynomial root finding since they are very sensitive to noise (in my experience, at least).

My current approach is to use the Hilbert transform:

$$ \operatorname{arg} \left[\hat{h}(\omega) \right] = - \mathcal{H} \left[ \ln | \hat{h}( \omega ) | \right] $$

This is how I've coded it in Python (I've padded the log magnitude with zeros as suggested by this post):

def min_phase_fir(freq_mag, zero_pad): logmag = np.log(freq_mag) logmag = np.pad(logmag, (0, zero_pad), 'constant', constant_values=(0, 0)) phase_response = -np.imag(hilbert(logmag)) return phase_response[:len(freq_mag)]

My questions are:

  1. How do I verify that this phase response is indeed the correct answer?
  2. Should I window the frequency magnitude at the point where I add zero-padding to ensure it's (somewhat) continuous? If yes, how?

Edit: Here is the impulse response from which the above magnitude was obtained, in case you need reference. Note it is not minimum-phase: impulse response

Improve this question
$\endgroup$
2
  • $\begingroup$ Do you have the time domain data magnitude? $\endgroup$ Commented Jul 25, 2024 at 15:36
  • $\begingroup$ @Baddioes I do have the time domain data that generated example magnitude data. However it's not minimum phase. I've added it to the question for reference. $\endgroup$ Commented Jul 25, 2024 at 15:55

2 Answers 2

Reset to default
2
$\begingroup$

How do I verify that this phase response is indeed the correct answer?

Visual inspection works pretty well here. Things that indicate problems are

  1. pre-ringing, i.e. unusual activity at the end of the impulse repsone
  2. echoes: attenuated copy of the impulse response somewhere in the middle.

Should I window the frequency magnitude at the point where I add zero-padding to ensure it's (somewhat) continuous? If yes, how?

That would definitely help here. You need to interpolate the spectrum. You can do this by

  1. Applying a zero phase
  2. Go back to the time domain
  3. Circulate the impulse response to the peak is in the middle (typically half the length).
  4. Pad zeros equally in front and rear
  5. Go back into the frequency domain. Discard the phase.

Then proceed with a Hilbert transform approach. You cannot implement an ideal Hilbert Transformer since it's infinitely non-causal. Hence there is always some approximation required and I typically roll my own so I have full control of the trade-offs.

In my case, the magnitude is quite short (~30-50 samples).

This will be your main problem. The original impulse response looks truncated which means you don't have a dense enough frequency grid and you are lacking information to fully characterize the system.

Improve this answer
$\endgroup$
1
  • $\begingroup$ What do you mean by applying zero phase? Computing $|\hat{h}(\omega)|^2$? $\endgroup$ Commented Jul 25, 2024 at 16:58
1
$\begingroup$

Here is a baby add-on to @Hilmar answer to your first question.

How do I verify that this phase response is indeed the correct answer?

You need to check if the reconstructed filter is indeed a minimum phase filter. Towards that aim, combine the magnitude and phase responses to reconstruct the frequency response and obtain the impulse response of the filter. To ensure the FIR filter is minimum phase (MP), check if all the zeros of the impulse response lie inside the unit circle in the z-plane. The python snippet code below can check if the reconstructed filter is MP or not

# Find the zeros of the polynomial represented by the impulse response zeros = np.roots(impulse_response) # Check if all zeros are inside the unit circle np.all(np.abs(zeros) < 1)

Make sure that you remove trailing zeros in impulse_response before computing the roots. Since this is an FIR filter, you can skip testing for stability. In case of an IIR filter, you also need to check for stability of the system (equivalent to check if the denominator polynomial is a minimum phase or not).

In MATLAB, you can use the isminphase function to check for minimum phase.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.