1
$\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$
0
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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