我想对具有任意采样率的时间序列进行加权。
模拟 A 加权滤波器由 IEC 61672-1 准确定义。但是没有数字滤波器的定义。一种方法是使用双线性变换 (BLT) 将模拟滤波器转换为数字滤波器（如此处应用 A 加权所做的）。然而，这种方法在 nyquist 附近存在极端扭曲（即使模拟极点/零点被预先扭曲）：
图 1：A 加权频率响应比较，采样率为 $25600\textrm{ Hz}$。
相反，我正在考虑使用一种算法来设计具有任意频率响应的数字 IIR 滤波器并插入模拟 A 加权滤波器的频率响应。
我已经研究了 MATLAB yulewalk，但我需要一个相应的 Python 实现来尝试。我还在几个地方遇到过 Berchin 的 FDLS 方法，例如这个问题，但所有链接似乎都已损坏。
一个常见的误解是，数字滤波器对模拟滤波器的逼近一定很差，接近奈奎斯特。这个想法可能来自无处不在的双线性变换，通常情况确实如此。当然，Nyquist 对离散时间滤波器的频率响应有一定的限制，但它们不一定会导致模拟滤波器在该频率范围内的近似错误。接近奈奎斯特的近似质量取决于几个因素，其中包括模拟滤波器的频率响应特性，以及是否只需要近似模拟滤波器的幅度或相位这一事实。
我设计了一个 6 阶 IIR 滤波器，近似于模拟 A 加权滤波器的频率响应，定义如下：
$$H(s)=\frac{k\cdot s^4}{(s+129.4)^2(s+676.7)(s+4636)(s+76655)^2}\tag{1}$$
$k=7.39705×10^9$。
我选择了 $48$ kHz 的采样频率。设计过程是我前段时间提出的启发式迭代过程。这是基于方程误差法的最小二乘近似，我可能有一天会写下所有细节。
下面是设计结果图。请注意，这两个图都上升到 Nyquist ($24$ kHz)。左图显示，在模拟和数字频率响应的对数图之间看不到任何差异。右图显示了近似误差，定义为幅度响应差的绝对值。该错误显示了典型的最小二乘行为。
您可以自己检查过滤器。以下是系数：
b = 0.169994948147430 0.280415310498794 -1.120574766348363 0.131562559965936 0.974153561246036 -0.282740857326553 -0.152810756202003 一个= 1.00000000000000000 -2.12979364760736134 0.42996125885751674 1.62132698199721426 -0.96669962900852902 0.00121015844426781 0.04400300696788968
取决于你需要多花哨。48 kHz 的良好实现如下
b = [0.234301792299513 -0.468603584599026 -0.234301792299513 ...
0.937207169198054 -0.234301792299515 -0.468603584599025 0.234301792299513];
a = [1.000000000000000 -4.113043408775871 6.553121752655047 ...
-4.990849294163381 1.785737302937573 -0.246190595319487 0.011224250033231];
如果您需要不同的采样率，您可以简单地计算 IIR 滤波器在 48 kHz 的脉冲响应，截断并重新采样到您想要的任何采样率。在 48 kHz 时，256 个抽头可让您降至约 50 Hz。1024 分接头低至 20 Hz（如果您愿意的话），而 2048 与您所需要的一样好。
重采样总是会影响奈奎斯特频率的频率响应，但这对于所有信号处理任务来说都是不争的事实。
原始过滤系数由 Ir Christophe Couvreur Faculte Polytechnique de Mons 博士提供
删除了用于垃圾邮件保护的详细信息，如果您需要联系信息，请告诉我。
以下是我使用单独的 HPF (c2d)) 和 LPF (MIM) 准备的 5 阶和 6 阶滤波器 (fs=48kHz) 的系数：
四阶 HPF，一阶 LPF：
[0.588402730019084 -2.114578549207340 2.574286896638522 -0.919416694862366 -0.367726753456895 0.239032370868995]
[1.0 -4.193437345479311 6.853084418397484 -5.397323870680988 2.009149234848459 -0.271472430592242]
四阶 HPF，二阶 LPF：
[[0.5884027300190836 -1.6384145461593100 0.8580879826180272 1.1837389307393180 -1.1416401766192728 -0.0386320187694699 0.1884570981716236]
[1.0 -3.384188885614779 3.453610828634393 0.171625555633360 -2.392296872830998 1.376227862809258 -0.224978476938553]
以及与他们的模拟模型进行比较的图表：
（绿色 = 模拟模型，红色 = 5 阶滤波器，黑色 = 6 阶滤波器，小窗口显示接近 Nyqvist 频率时的响应差异）
编辑2：
这是构建过滤器的 Octave 源代码：
% Octave packages -------------------------------
pkg load signal
% other requirements:
% Magnitude Invariance method (MIM) -implementation ( one implementation can be found from https://soar.wichita.edu/handle/10057/1564 )
clf;
format long;
%Sampling Rate
Fs = 44100;
% A-weighting filter frequencies according to IEC/CD 1672.
% Source: https://www.dsprelated.com/showcode/214.php
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
% Analog model -----------------------------------------------------
A1000 = 1.9997;
NUM = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0];
DEN = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
DEN = conv(conv(DEN,[1 2*pi*f3]),[1 2*pi*f2]);
Ds = tf(NUM, DEN);
[Ab, Aa, T] = tfdata(Ds, 'v');
% Digital filter --------------------------------------------------
%LPF1 (MIM method)
lporder = 1; % 1 = 5th order A-Weighting filter, 2= 6th order ....
w0 = 2 * pi * f4;
bC = 1;
aC = [1 w0];
AD0 = tf(bC, aC);
LP2 = c2dn(AD0^2, 1/Fs, 'mim', lporder, 4096^2, 'lowpass');
% HPF (BLT method)
% HPF1
w0 = 2 * pi * f1;
bC = [w0 0];
aC = [1 w0];
AD1 = tf(bC, aC);
HP1 = c2d(AD1^2, 1/Fs, 'tustin');
% HPF2
w0 = 2 * pi * f2;
bC = [w0 0];
aC = [1 w0];
AD2 = tf(bC, aC);
HP2 = c2d(AD2, 1/Fs, 'tustin');
% HPF3
w0 = 2 * pi * f3;
bC = [w0 0];
aC = [1 w0];
AD3 = tf(bC, aC);
HP3 = c2d(AD3, 1/Fs, 'tustin');
% Combine filters
FLT = (LP2*HP1*HP2*HP3);
% Adjust 0dB@1kHz
GAINoffset = abs(freqresp(FLT,1000*2*pi));
FLT = FLT/GAINoffset;
% get coefficients
[ad, bd, T] = tfdata(FLT, 'v');
ad
bd
% plot
fs2 = Fs/2;
nf = logspace(0, 5, fs2);
[mag, pha] = bode(Ds,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'g', 'linewidth', 4.0, 'linestyle', '-');
hold on;
[mag, pha] = bode(FLT,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'k', 'linewidth', 2.0, 'linestyle', '--');
title('A-Weighting filter');
xlabel('Hz');ylabel('dB');
axis([1 fs2+5000 -200 10]);
legend('Analog model', 'Digital (MIM+BLT)', 'location', 'southeast');
grid on;
我没有测量误差曲线（不知道如何根据模拟模型计算它）。
EDIT1：显示低采样率（4、8、16、32 kHz）的高频响应的图：
很久以前，我制作了自己的 matlab 例程，为我提供了用于 A 加权的 IIR 滤波器系数。从一个现已失效的链接中，我得到了表单的传递函数：
ka*s^4
------------------------------------------------
(s+129.4)^2 * (s+676.7) * (s+4636) * (s+76655)^2
从这个出发点，我使用 matlab 构建了一个离散版本的过滤器。这很丑陋，但在这里：
function [b,a]=makeAweightingFilter(fs);
%build the continuous-time transfer function
num=poly([0 0 0 0]); %numerator of transfer function (TF)
den=poly(-129.4336); %one of the poles of the TF
den=conv(den,poly(-129.4336)); %another pole of the TF
den=conv(den,poly(-676.6991)); %another pole of the TF
den=conv(den,poly(-4636.3624)); %another pole of the TF
den=conv(den,poly(-76654.86075)); %another pole of the TF
den=conv(den,poly(-76654.86075)); %another pole of the TF
a_weight=tf(num,den); %form the transfer function
%Scale transfer function to yield 0dB response at 1000 Hz
offset = abs(freqresp(a_weight,1000*2*pi));
num=1/offset*poly([0 0 0 0]); %scale the numerator
a_weight=tf(num,den); %re-form the transfer function
%Convert to discrete time system
a_w_discrete = c2d(a_weight,1/fs,'zoh');
%get filter coefficients from the numerator
%and denominator of the discrete transfer function
b=a_w_discrete.num{1};
a=a_w_discrete.den{1};