scipy.signal.

minimum_phase#

scipy.signal.minimum_phase(h, method='homomorphic', n_fft=None, *, half=True)[源码]#

将线性相位 FIR 滤波器转换为最小相位滤波器

参数:
harray

线性相位 FIR 滤波器系数。

method{‘hilbert’, ‘homomorphic’}

提供的方法有:

‘homomorphic’(默认)

该方法 [4] [5] 对具有奇数个抽头的滤波器效果最佳。当 half=True(默认)时,生成的最小相位滤波器抽头数减半,幅度响应近似于原滤波器幅度响应的平方根;当 half=False 时,抽头数与原滤波器相同,幅度频谱与原滤波器一致。

‘hilbert’

该方法 [1] 旨在用于具有单位或零增益区域的等波纹滤波器(例如,由 remez 设计的滤波器)。

n_fftint

用于 FFT 的点数。应至少比信号长度大几倍(参见“注意”)。

halfbool

如果为 True,则创建一个长度为原滤波器一半、幅度频谱为原滤波器平方根的滤波器。如果为 False,则创建一个长度与原滤波器相同、幅度频谱与原滤波器匹配的滤波器(仅在 method='homomorphic' 时受支持)。

版本 1.14.0 中新增。

返回:
h_minimumarray

滤波器的最小相位版本。当 half is True 时,长度为 (len(h) + 1) // 2,否则为 len(h)

另请参阅

firwin
firwin2
remez

附注

希尔伯特(Hilbert)方法 [1] 或同态(homomorphic)方法 [4] [5] 都需要选择 FFT 长度来估计滤波器的复倒谱。

在希尔伯特方法中,与理想频谱的偏差 epsilon 与阻带零点数 n_stop 和 FFT 长度 n_fft 的关系如下:

epsilon = 2. * n_stop / n_fft

例如,对于 100 个阻带零点和 2048 的 FFT 长度,epsilon = 0.0976。如果我们保守地假设阻带零点数比滤波器长度少一个,我们可以将 FFT 长度取为满足 epsilon=0.01 的下一个 2 的幂次:

n_fft = 2 ** int(np.ceil(np.log2(2 * (len(h) - 1) / 0.01)))

这对于希尔伯特方法和同态方法都能提供合理的结果,也是当 n_fft=None 时所使用的值。

创建最小相位滤波器还存在其他实现方式,包括零点反转 [2] 和谱因式分解 [3] [4]。有关更多信息,请参阅 此 DSPGuru 页面

数组 API 标准支持

minimum_phase 除了支持 NumPy 之外,还对 Python Array API 标准兼容的后端提供实验性支持。请考虑通过设置环境变量 SCIPY_ARRAY_API=1 并提供 CuPy、PyTorch、JAX 或 Dask 数组作为数组参数来测试这些功能。支持以下后端和设备(或其他功能)的组合。

CPU

GPU

NumPy

不适用

CuPy

不适用

PyTorch

JAX

⚠️ 无 JIT

Dask

⚠️ 计算图

不适用

有关更多信息,请参阅 对数组 API 标准的支持

参考文献

[1] (1,2)

N. Damera-Venkata and B. L. Evans, “Optimal design of real and complex minimum phase digital FIR filters,” Acoustics, Speech, and Signal Processing, 1999. Proceedings., 1999 IEEE International Conference on, Phoenix, AZ, 1999, pp. 1145-1148 vol.3. DOI:10.1109/ICASSP.1999.756179

[2]

X. Chen and T. W. Parks, “Design of optimal minimum phase FIR filters by direct factorization,” Signal Processing, vol. 10, no. 4, pp. 369-383, Jun. 1986.

[3]

T. Saramaki, “Finite Impulse Response Filter Design,” in Handbook for Digital Signal Processing, chapter 4, New York: Wiley-Interscience, 1993.

[4] (1,2,3)

J. S. Lim, Advanced Topics in Signal Processing. Englewood Cliffs, N.J.: Prentice Hall, 1988.

[5] (1,2)

A. V. Oppenheim, R. W. Schafer, and J. R. Buck, “Discrete-Time Signal Processing,” 3rd edition. Upper Saddle River, N.J.: Pearson, 2009.

示例

创建一个过渡带为 [0.2, 0.3] 的最佳线性相位低通滤波器 h(假设奈奎斯特频率为 1)

>>> import numpy as np
>>> from scipy.signal import remez, minimum_phase, freqz, group_delay
>>> import matplotlib.pyplot as plt
>>> freq = [0, 0.2, 0.3, 1.0]
>>> desired = [1, 0]
>>> h_linear = remez(151, freq, desired, fs=2)

将其转换为最小相位

>>> h_hil = minimum_phase(h_linear, method='hilbert')
>>> h_hom = minimum_phase(h_linear, method='homomorphic')
>>> h_hom_full = minimum_phase(h_linear, method='homomorphic', half=False)

比较这四个滤波器的脉冲响应和频率响应

>>> fig0, ax0 = plt.subplots(figsize=(6, 3), tight_layout=True)
>>> fig1, axs = plt.subplots(3, sharex='all', figsize=(6, 6), tight_layout=True)
>>> ax0.set_title("Impulse response")
>>> ax0.set(xlabel='Samples', ylabel='Amplitude', xlim=(0, len(h_linear) - 1))
>>> axs[0].set_title("Frequency Response")
>>> axs[0].set(xlim=(0, .65), ylabel="Magnitude / dB")
>>> axs[1].set(ylabel="Phase / rad")
>>> axs[2].set(ylabel="Group Delay / samples", ylim=(-31, 81),
...             xlabel='Normalized Frequency (Nyqist frequency: 1)')
>>> for h, lb in ((h_linear,   f'Linear ({len(h_linear)})'),
...               (h_hil,      f'Min-Hilbert ({len(h_hil)})'),
...               (h_hom,      f'Min-Homomorphic ({len(h_hom)})'),
...               (h_hom_full, f'Min-Homom. Full ({len(h_hom_full)})')):
...     w_H, H = freqz(h, fs=2)
...     w_gd, gd = group_delay((h, 1), fs=2)
...
...     alpha = 1.0 if lb == 'linear' else 0.5  # full opacity for 'linear' line
...     ax0.plot(h, '.-', alpha=alpha, label=lb)
...     axs[0].plot(w_H, 20 * np.log10(np.abs(H)), alpha=alpha)
...     axs[1].plot(w_H, np.unwrap(np.angle(H)), alpha=alpha, label=lb)
...     axs[2].plot(w_gd, gd, alpha=alpha)
>>> ax0.grid(True)
>>> ax0.legend(title='Filter Phase (Order)')
>>> axs[1].legend(title='Filter Phase (Order)', loc='lower right')
>>> for ax_ in axs:  # shade transition band:
...     ax_.axvspan(freq[1], freq[2], color='y', alpha=.25)
...     ax_.grid(True)
>>> plt.show()
../../_images/scipy-signal-minimum_phase-1_00_00.png
../../_images/scipy-signal-minimum_phase-1_00_01.png

脉冲响应和群延迟图描述了线性相位滤波器 h 的 75 个样本延迟。相位在阻带也应该是线性的——由于数值较小,那里的数值噪声占主导地位。此外,这些图表清楚地显示最小相位滤波器在通带和过渡带具有减小的(负)相位斜率。图表还说明了参数为 method='homomorphic', half=False 的滤波器具有与线性滤波器 h 相同的阶数和幅度响应,而其他最小相位滤波器只有一半的阶数和幅度响应的平方根。