通过希尔伯特变换返回值的最小相位相位

信息处理 Python 阶段 希尔伯特变换 最小相位 hrtf
2022-01-30 12:48:35

继我之前的问题：HRIR 最小相位之后，我设法计算了 FIR 滤波器（在我的特定情况下是 HRTF 滤波器）的最小相位相位。

但是我不确定我的函数返回的相位值。

Python函数如下。它期望一个任意的 HRIR 作为输入。我正在使用 44100 长 FFT 使频率箱在 44100Hz 采样率的整数频率下完美对齐（例如，1000Hz 的正弦曲线的幅度正好为 1，没有频谱泄漏）

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()

对于下面的图，我使用的是KEMAR MIT数据集，特别是已经为 HRIR 计算了图[azimuth=144°, elevation=30°, distance=1.4m]。选择的 HRIR 如下（每只耳朵一个）： 高分辨率红外

震级

一旦我使用上面的 Python 函数执行我的最小相位计算，这就是我得到的（单耳的相位）： 阶段

但是，看情节，我无法理解返回的值

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

返回的阶段是否展开？还是我做错了什么？一旦超过 +3.14 或低于 -3.14 (+/- pi)，原始相位就会按预期“环绕”，而最小相位则不会。

2个回答

我计算并比较了最小相位 HRIR 和原始相位。

这是我的最终代码：

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

从情节来看，我得到的结果似乎是合理的： 中相脉冲响应

并且幅度与我预期的相同（原始阶段和最小阶段一完全重叠）：

震级

但是，我并不完全确定相位行为（如我的问题中所述）及其返回值。

抱歉，很抱歉，但我认为解决这个问题可能会引起普遍的兴趣，因为最小阶段是一件很酷的事情。

展开阶段：

这似乎正是如此。函数虚部和实部之间的潜在Kramers-Kroning关系是因果信号的一般属性。不涉及相位缠绕。

有事吗：

我可以验证你的实现是正确的，虽然首先考虑虚部是令人费解的，当从最小阶段的简单定义开始时Φ0以及频谱的幅度，您将其定义为|HRIRfft|：

Φ0=−H(log(|HRIRfft|)

查看希尔伯特的 Scipy 文档将揭开这个谜团。

其它你可能感兴趣的问题
上一篇如何优化归一化滤波器系数？ 下一篇傅里叶变换作为参数乘以齐次波的积分