minimum_phase#
- scipy.signal.minimum_phase(h, method='homomorphic', n_fft=None, *, half=True)[源代码]#
将线性相位 FIR 滤波器转换为最小相位
- 参数:
- harray
线性相位 FIR 滤波器系数。
- method{‘hilbert’, ‘homomorphic’}
提供的方法有
- n_fftint
用于 FFT 的点数。应该至少比信号长度大几倍(请参阅“注释”)。
- halfbool
如果
True
,则创建一个长度为原始滤波器一半的滤波器,其幅度频谱是原始滤波器的平方根。如果False
,则创建一个与原始滤波器长度相同的滤波器,其幅度频谱旨在与原始滤波器匹配(仅在method='homomorphic'
时支持)。1.14.0 版本新增。
- 返回:
- h_minimumarray
滤波器的最小相位版本,当
half 为 True
时,长度为(len(h) + 1) // 2
,否则为len(h)
。
注释
希尔伯特 [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 年。
示例
创建一个过渡带为 [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()
脉冲响应和群延迟图描绘了线性相位滤波器 h 的 75 个样本延迟。相位在阻带中也应该是线性的——由于幅度很小,数值噪声在那里占主导地位。此外,这些图表明,最小相位滤波器在通带和过渡带中清晰地显示出减小(负)的相位斜率。这些图还说明了参数为
method='homomorphic', half=False
的滤波器与线性滤波器 h 具有相同的阶数和幅度响应,而其他最小相位滤波器只有一半的阶数和幅度响应的平方根。