3
$\begingroup$

Following my previous question: HRIR Minimum phase I managed to compute the minimum-phase phase of a FIR filter (in my particular case, HRTF filters).

However I am not sure of the phase values returned by my function.

The Python function is the following. It expects an arbitrary HRIR as input. I am using a 44100 long FFT to have the freqency bins perfectly aligned at integer frequencies at 44100Hz sampling rate (so for example a sinusoid at 1000Hz will have a magnitude of exactly 1 with no spectral leakage)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
from scipy.fft import fft, ifft, fftfreq

def min_phase_coversion(HRIR):

    HRIR_fft = fft(HRIR,44100)

    xf = fftfreq(len(HRIR_fft), 1/44100)

    phase=np.angle(HRIR_fft)

    minimum_phase = np.imag(-hilbert(np.log(np.abs(HRIR_fft))))


    plt.figure()
    plt.plot(xf, phase)
    plt.plot(xf,minimum_phase)
    plt.plot()
    plt.grid()
    plt.show()

For the following plots I am using the KEMAR MIT dataset, in particular the plot have been computed for the HRIR at [azimuth=144°, elevation=30°, distance=1.4m]. The selected HRIR is the following (one per ear): HRIR

Magnitude

Once I perform my minimum-phase computation using the above Python function, this is what I get (phase for a single ear): Phase

However, looking at the plot, I can't understand the values returned by

np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

Is the returned phase unwrapped? Or am i doing something wrong? The original phase "wraps" as expected once it gets over +3.14 or under -3.14 (+/- pi), while the minimum one doesn't.

$\endgroup$

2 Answers 2

1
$\begingroup$

I computed and compared the minimum phase HRIR and the original one.

This is my final code:

def min_phase_coversion(HRIR):
    '''

    :param HRIR: the desired HRIR impulse response to convert into minimum phase
    :return: the minimum phase version of the original HRIR
    '''

    HRIR_fft = fft(HRIR,44100)

    #computing magnitude, tested with sinusoid, it works.
    xf = fftfreq(len(HRIR_fft), 1/44100)

    fft_magnitude = np.abs(HRIR_fft/len(HRIR_fft))

    phase=np.angle(HRIR_fft)

    minimum_phase = np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

    plt.figure()
    plt.plot(xf, fft_magnitude) #computing and plotting magnitude
    plt.plot(xf, phase, label='original phase')
    plt.plot(xf,minimum_phase, label='minimum phase')
    plt.legend()
    plt.plot()
    plt.grid()
    plt.show()

    min_phase_HRIR = irfft(np.abs(HRIR_fft)*np.exp(1j*minimum_phase), 44100) # |H(w)|*e^(jPhi(w))

    HRIR_original = irfft(HRIR_fft, 44100)


    plt.figure()
    plt.plot(min_phase_HRIR[0:512], label="Min_phase", linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.plot(HRIR_original[0:512], label="original phase", linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.legend()
    plt.show()


    plt.figure()
    min_phase_fft = fft(min_phase_HRIR, 44100)
    min_phase_mangnitude =  np.abs(min_phase_fft/len(min_phase_fft))
    plt.plot(xf,fft_magnitude)
    plt.plot(xf,min_phase_mangnitude)
    plt.show()

    return min_phase_HRIR

The results I get seem reasonable looking at the plots: minphase impulse response

and the magnitude is the same as I expected (the original one and the min phase one are perfectly overlapping):

magnitude

However, I am not totally sure of the phase behaviour (as mentioned in my question) and its returned values.

$\endgroup$
2
  • $\begingroup$ Hi, I am very interested in this question. Finally did you opt for this implementation? $\endgroup$
    – Chutlhu
    Jan 31 at 9:06
  • $\begingroup$ Yes, looks reasonable, but to this day im still not 100% sure $\endgroup$ Jan 31 at 11:02
0
$\begingroup$

Sorry for necobumping, but I think resolving this may be of general interest, because the minimum phase is quite a cool thing.

unwrapped phase:

This seems to be exactly the case. The underlying Kramers-Kroning relation between imaginary part and real part of a function is a general property of causal signals. No phase wrapping involved.

something wrong:

I can verify that you implementation is correct, although taking the imaginary part is puzzling at first, when starting with the simple definition of the minimum phase $\Phi_0$ and the magnitude of the spectrum, which you defined as $|\mathrm{HRIRfft}|$:

$$ \Phi_0=-\mathscr{H}(log(|\mathrm{HRIRfft}|) $$

Looking at the Scipy docs for hilbert will lift that mystery.

$\endgroup$
1
  • $\begingroup$ So at the end of the day it was correct? $\endgroup$ Jun 27, 2023 at 12:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.