多通道滤波器截止校正
5.00/5 (9投票s)
澄清多通道滤波器截止校正的使用
引言
由于多次滤波导致截止频率的衰减，零相位移（零延迟）IIR滤波器的截止频率应进行校正。 这也适用于MATLAB的filtfilt函数，但在文档中没有直接提及。
似乎很少有人意识到这种校正，而意识到的人往往使用不正确的实现。 因此，我选择发布我关于这个问题的发现，以使人们意识到它并指向正确的参考。 我使用了MATLAB作为示例，但它适用于所有实现。
理论
正如David A. Winter在他的书《人体运动的生物力学和运动控制，第四版》（第69页），以及D. Gordon E. Robertson和他的同事在他们的书《生物力学研究方法，第二版》（第288页）中所描述的那样，必须对零相位移Butterworth滤波器的截止频率进行校正。 这是因为为了补偿每次不均匀滤波通过引入的相移，滤波器再次在反向时间方向上运行，并且对于每次滤波通过，-3db截止频率将被推低。
Winter展示了如何使用校正因子来补偿二阶Butterworth滤波器的截止频率衰减，该校正因子应用于调整后的角截止频率（而不是直接应用于角截止频率）。 Winter用于调整后的角截止频率的变量名有点误导，因此我使用了与Robertson等人使用的名称相似的变量名（仅用U代替Ω）。
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter Butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
请注意，这些作者的旧版本书籍和论文可能不使用此更新的校正方法。
背景
如果你在互联网上搜索这个截止调整，你会发现许多人错误地实现了这个方法。 它不是将校正因子应用于调整后的角截止频率 (Uc)，而是将其应用于截止频率。 例如，在一些情况下，对于双通滤波器，校正就像这样完成：f_corrected = f_cutoff / 0.802。 这将给出不正确的响应，尽管它在较低频率下相差不大。
为了研究如果直接应用于截止频率，校正因子会是什么样子，可以从校正后的参数Un中推导出截止频率（参见结果部分，了解如何完成）。 事实证明，它是频率相关的。 让我们将这个频率相关的因子称为Cf = f_cutoff / f_corrected。 如果有人知道如何直接获得这个因子，请告诉我。
结果
根据Winter的说法，获得二阶低通Butterworth系数的正确方法是
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
k1 = sqrt(2) * Un;
k2 = Un^2;
a0 = k2 / (1 + k1 + k2);
a1 = 2 * a0;
a2 = a0;
k3 = a1 / k2;
b1 = k3 - a1;
b2 = 1 - a1 - k3;
b = [ a0 a1 a2 ];
a = [ 1 -b1 -b2 ];
data_filtered = myfilter(data,b,a,'low',filter_passes);
如果要使用仅将截止频率作为输入的现有函数（如MATLAB butter(..)函数），则可以将校正后的参数Un转换回频率。
filter_passes = 2; % filtfilt use 2 passes
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
f_corrected = atan(Un)*sampling_frequency/pi; % corrected cutoff frequency
[b, a] = butter(2, f_corrected / (sampling_frequency / 2), 'low');
data_filtered = filtfilt(b, a, data);
现在我们可以测试参考单通滤波器、未校正的多通滤波器、INCORRECT频率校正的多通滤波器和Winter校正的多通滤波器之间的差异。 为了测试多次滤波通过会发生什么，使用了以滤波通过次数作为输入的滤波器（参见附件）。 所示的图形是滤波后的噪声频谱。
通过两次滤波，未校正的截止频率变为134Hz而不是150Hz，错误校正的截止频率变得太高，而Winter提出的校正后的截止频率正如以下图中所示的那样准确。
随着滤波通过次数的增加，差异变得更大（在4次滤波通过时，错误校正的截止频率在fs/3处已经达到fs/2）。
同样的原理也适用于高通滤波器。
要获得正确的二阶高通Butterworth系数，您可以使用Robertson等提出的方法
cutoff_frequency = sampling_frequency/2 - cutoff_frequency; % Robertson high-pass modification
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc / C; % David A. Winter correction
k1 = sqrt(2) * Un;
k2 = Un^2;
a0 = k2 / (1 + k1 + k2);
a1 = 2 * a0;
a2 = a0;
k3 = a1 / k2;
b1 = k3 - a1;
b2 = 1 - a1 - k3;
% Robertson high-pass modification
a1 = -a1;
b1 = -b1;
b = [ a0 a1 a2 ];
a = [ 1 -b1 -b2 ];
data_filtered = myfilter(data,b,a,'high',filter_passes);
您可以像低通滤波器一样，通过乘以校正因子（而不是除以它）来获得校正后的高通截止频率。
filter_passes = 2; % filtfilt use 2 passes
C = (((2^(1/filter_passes))-1)^(1/4)); % David A. Winter butterworth correction factor
Wc = 2*pi*cutoff_frequency; % angular cutoff frequency
Uc = tan(Wc/(2*sampling_frequency)); % adjusted angular cutoff frequency
Un = Uc * C; % multiply by C for highpass correction
f_corrected = atan(Un)*sampling_frequency/pi; % corrected cutoff frequency
[b, a] = butter(2, f_corrected / (sampling_frequency / 2), 'high');
data_filtered = filtfilt(b, a, data);
如果我们看一下高通滤波器的差异，我们发现不正确的频率校正不像低通滤波器那样糟糕，但在更高的截止频率和更多的滤波通过次数下会变得更糟。
讨论
重要的是要注意，这个特定的校正因子 (C) 仅适用于二阶 Butterworth 系数，不适用于更高阶数。 由于更锐利的截止曲线，更高阶的滤波器不会受到截止频率衰减的相同程度的影响。 然而，它仍然受到影响，需要考虑。 如果有人知道考虑滤波器阶数的校正因子，请告诉我。
历史
- 首次发布