如何将Matlab fdesign.bandpass参数转换为scipy.signal.buttord对应参数

核心疑问梳理 #

当前复现过程存在两处待解决的技术问题：

• Matlab fdesign.bandpass 接口的第一阻带衰减Ast1、通带波纹Ap、第二阻带衰减Ast2参数，如何对应scipy.signal.buttord仅支持的通带最大损耗gpass、阻带最小衰减gstop参数？
• Matlab中群延迟（group delay）的计算与应用逻辑如何用Python实现？

原参考Matlab实现代码如下：

%eegFS = 2000;  % Signal freq, Hz
%lcut = 5;  % low cut freq, Hz
%hcut = 100;  % high cut freq, Hz
%attenHz = 4;  % transition band
%attendB = 40; % attenuation

nyq = round(eegFS/2);

%make bandpass
Fstop1 = (lcut - attenHz)  / nyq;  % First Stopband Frequency
Fpass1 = lcut  / nyq;  % First Passband Frequency
Fpass2 = hcut / nyq;  % Second Passband Frequency
Fstop2 = (hcut + attenHz) / nyq;  % Second Stopband Frequency
Astop1 = attendB;    % First Stopband Attenuation (dB)
Apass  = 1;     % Passband Ripple (dB)
Astop2 = attendB;    % Second Stopband Attenuation (dB)
h = fdesign.bandpass('fst1,fp1,fp2,fst2,ast1,ap,ast2', Fstop1, Fpass1, ...
                 Fpass2, Fstop2, Astop1, Apass, Astop2);
Hd = design(h, 'kaiserwin');   
b = Hd.Numerator;

%group delay
[a,f] = grpdelay(b,1,nyq,eegFS);
k = f >= lcut & f <= hcut;
gd = fix(mean(a(k)));

% apply filter
eegF = filter(b, 1, eegData, [], 2);
eegF = cat(2,eegF(:,gd+1:end),zeros(size(eegF,1),gd));

问题详解 #

1 参数对应关系 #

首先需要注意：你提供的Matlab代码最终采用凯撒窗设计FIR滤波器，若要完全对齐效果无需使用巴特沃斯滤波器的buttord接口，直接使用scipy.signal.kaiserord计算滤波器阶数即可，参数对应关系完全匹配，无需额外转换。
如果确实需要使用巴特沃斯滤波器实现，参数对应规则如下：

• scipy.signal.buttord的gpass参数直接对应Matlab的通带波纹Apass，即通带内允许的最大衰减值，单位为dB
• scipy.signal.buttord的gstop参数对应Matlab的Astop1/Astop2：如果两个阻带衰减要求一致，直接填入统一值即可；如果Astop1和Astop2数值不同，取两者的最小值作为gstop，保证两个阻带的衰减要求都能被满足。

2 群延迟的Python实现 #

Matlab的grpdelay接口对应scipy.signal.group_delay，计算逻辑完全对齐，具体实现步骤如下：

1. 调用scipy.signal.group_delay传入滤波器系数，得到所有频率点的群延迟数值
2. 筛选通带范围内的群延迟值，取均值后向下取整得到需要补偿的偏移量gd
3. 滤波完成后裁去数组前gd个偏移的采样点，末尾补对应数量的0即可对齐原数组长度。

完整Python复现代码 #

import numpy as np
from scipy import signal

# 滤波器参数配置
eegFS = 2000  # 采样频率，单位Hz
lcut = 5      # 通带下限频率，单位Hz
hcut = 100    # 通带上限频率，单位Hz
attenHz = 4   # 过渡带宽度，单位Hz
attendB = 40  # 阻带最小衰减，单位dB
nyq = eegFS / 2

# 归一化频率计算
Fstop1 = (lcut - attenHz) / nyq
Fpass1 = lcut / nyq
Fpass2 = hcut / nyq
Fstop2 = (hcut + attenHz) / nyq
Apass = 1

# 凯撒窗FIR滤波器设计，对齐原Matlab实现
transition_width = min(Fpass1 - Fstop1, Fstop2 - Fpass2)
numtaps, beta = signal.kaiserord(attendB, transition_width)
# 确保滤波器阶数为奇数，避免群延迟出现半采样点偏移
numtaps = numtaps if numtaps % 2 == 1 else numtaps + 1
taps = signal.firwin(numtaps, [Fpass1, Fpass2], window=('kaiser', beta), pass_zero='bandpass')

# 群延迟计算
freq, gd_vals = signal.group_delay((taps, 1), fs=eegFS)
passband_mask = (freq >= lcut) & (freq <= hcut)
gd = int(np.mean(gd_vals[passband_mask]))

# 滤波与延迟补偿，假设eegData为形状（通道数, 采样点数）的EEG数据
eegF = signal.lfilter(taps, 1, eegData, axis=-1)
# 裁剪偏移点并补零对齐长度
eegF = np.concatenate([eegF[:, gd:], np.zeros((eegF.shape[0], gd), dtype=eegF.dtype)], axis=-1)

内容的提问来源于stack exchange，提问作者Tristan Stöber