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 True 时,长度为 (len(h) + 1) // 2,否则为 len(h)

另请参阅

firwin
firwin2
remez

注释

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

在希尔伯特方法的情况下,与理想频谱的偏差 epsilon 与阻带零点的数量 n_stop 和 FFT 长度 n_fft 相关,如下所示

epsilon = 2. * n_stop / n_fft

例如,对于 100 个阻带零点和 2048 的 FFT 长度,epsilon = 0.0976。如果我们保守地假设阻带零点的数量比滤波器长度少 1,则可以将 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 页面

参考文献

[1] (1,2)

N. Damera-Venkata 和 B. L. Evans,“实数和复数最小相位数字 FIR 滤波器的优化设计”,《声学、语音和信号处理,1999 年。论文集,1999 年 IEEE 国际会议》,菲尼克斯,亚利桑那州,1999 年,第 1145-1148 页,第 3 卷。DOI:10.1109/ICASSP.1999.756179

[2]

X. Chen 和 T. W. Parks,“通过直接因式分解设计最佳最小相位 FIR 滤波器”,《信号处理》,第 10 卷,第 4 期,第 369-383 页,1986 年 6 月。

[3]

T. Saramaki,“有限脉冲响应滤波器设计”,《数字信号处理手册》第 4 章,纽约:Wiley-Interscience,1993 年。

[4] (1,2,3)

J. S. Lim,《信号处理高级主题》。新泽西州恩格尔伍德悬崖:Prentice Hall,1988 年。

[5] (1,2)

A. V. Oppenheim、R. W. Schafer 和 J. R. Buck,《离散时间信号处理》,第 3 版。新泽西州上马鞍河: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 具有相同的阶数和幅度响应,而其他最小相位滤波器只有一半的阶数和幅度响应的平方根。