本周遇到了一个真正的挠头问题。我正在 C# 中实现一个 IIR 过滤器，所以我直接从 matlab 源代码复制以filter它们的时域过滤器函数(直接形式 II 转置)：
// direct form ii transposed
for (int i = 0; i < data.Length; i++)
{
Xi = data[i];
Yi = b[0] * Xi + Z[0];
for (int j = 1; j < order; j++)
Z[j - 1] = b[j] * Xi + Z[j] - a[j] * Yi;
Z[order - 1] = b[order] * Xi - a[order] * Yi;
output[i] = Yi;
}
return output;
奇怪的是，当我用脉冲测试滤波器时，我得到的值与 Matlab 报告的值略有不同。我也从 Matlab 获得滤波器系数。这是代码：
[b,a] = butter(3, [.0360, .1160], 'bandpass');
x = zeros(50,1);
x(1) = 1.0;
y = filter(b,a,x)
我在 C# 代码中使用b和a中的值作为系数。
Matlab 报告的y的前几个值：
>> y(1:13)
ans =
0.0016
0.0084
0.0216
0.0368
0.0487
0.0537
0.0501
0.0382
0.0194
-0.0038
-0.0286
-0.0519
-0.0713
由于这与我的 C# 端口不同，我直接将代码从filter复制到 C 文件，并使用相同的系数在那里运行它。输出与我在 C# 实现中获得的脉冲响应版本完全相同，略有不同：
[0] 0.0016 double
[1] 0.0086161600000000012 double
[2] 0.022182403216000009 double
[3] 0.038161063110481647 double
[4] 0.051323531488129848 double
[5] 0.05827273642334313 double
[6] 0.057456579295617483 double
[7] 0.048968543791003127 double
[8] 0.034196988694833064 double
[9] 0.015389667539999874 double
[10] -0.0048027826346631469 double
[11] -0.023749640330880527 double
[12] -0.039187648694732236 double
[13] -0.04946710058803272 double
我仔细查看了源以filter，我没有看到任何在计算滤波器输出之前按摩系数的证据。filter仅在a[0]不等于 1 的情况下对前馈系数进行归一化(在这种情况下，它肯定会这样做)。除此之外，我希望从 Matlab 的 C 代码中看到与从 Matlab 中输出完全相同的过滤器输出。
我真的很想知道这种差异来自哪里，因为我需要确信我的过滤器完全正确(我们不是都......我已经检查并重新检查了我的过滤系数。它们在 C/C# 和 Matlab 之间是相同的。
我用来获取这些"错误"值的完整 C 文件如下。我尝试了作为固定数量的状态实现的过滤器(在本例中为 6 个)和 N 个状态的一般情况(此处注释掉)。两者都来自 Matlab 源代码，并且都产生相同的"错误"输出：
# define order 6
# define len 50
int main(void){
// filter coeffs
float a[7] = { 1.0000, -5.3851, 12.1978, -14.8780, 10.3077, -3.8465, 0.6041 };
float b[7] = { 0.0016, 0.0, -0.0047, 0.0, 0.0047, 0, -0.0016 };
float a1 = a[1], a2 = a[2], a3 = a[3], a4 = a[4], a5 = a[5], a6 = a[6];
// input, output, and state arrays
float X[len];
float Y[len];
float Z[order];
float Xi;
float Yi;
float z0, z1, z2, z3,z4, z5;
// indeces
int i,j;
// initialize arrays
for(i=0;i<len;i++) {
X[i] = 0.0;
Y[i] = 0.0;
}
X[0] = 1.0;
for(i=0;i<order;i++)
Z[i] = 0.0;
z0 = Z[0];
z1 = Z[1];
z2 = Z[2];
z3 = Z[3];
z4 = Z[4];
z5 = Z[5];
i = 0;
while (i < len) {
Xi = X[i];
Yi = b[0] * Xi + z0;
z0 = b[1] * Xi + z1 - a1 * Yi;
z1 = b[2] * Xi + z2 - a2 * Yi;
z2 = b[3] * Xi + z3 - a3 * Yi;
z3 = b[4] * Xi + z4 - a4 * Yi;
z4 = b[5] * Xi + z5 - a5 * Yi;
z5 = b[6] * Xi - a6 * Yi;
Y[i++] = Yi;
}
//// copied from matlab filter source code
//i=0;
//while (i < len) {
// Xi = X[i]; // Get signal
// Yi = b[0] * Xi + Z[0]; // Filtered value
// for (j = 1; j < order; j++) { // Update conditions
// Z[j - 1] = b[j] * Xi + Z[j] - a[j] * Yi;
// }
// Z[order - 1] = b[order] * Xi - a[order] * Yi;
//
// Y[i++] = Yi; // Write to output
// }
}
如果我使用滤波器系数的所有数字，我会得到 Matlab 答案。
具体说来
float a[7] = {1,
-5.3850853610906082025167052052,
12.1978301571107792256043467205,
-14.8779557262737220924009307055,
10.3076512098041828124905805453,
-3.84649525959781790618308150442,
0.604109699507274999774608659209};
float b[7] = {0.00156701035058826889344307797813, 0,
-0.0047010310517648064634887994373, 0,
0.0047010310517648064634887994373, 0,
-0.00156701035058826889344307797813};
(使用 Scipy 获得的滤波系数：[b,a] = scipy.signal.butter(3, [.0360, .1160], 'bandpass')因为我不是由 . 💸
有了这个，你的C代码可以打印出来：
[0] = 0.0015670103020966053009033203125
[1] = 0.00843848474323749542236328125
[2] = 0.02162680588662624359130859375
[3] = 0.036844909191131591796875
[4] = 0.048709094524383544921875
[5] = 0.05368389189243316650390625
[6] = 0.05014741420745849609375
[7] = 0.0382179915904998779296875
[8] = 0.0194064676761627197265625
[9] = -0.003834001719951629638671875
与您的 Matlab 输出相匹配。
与浮点/双倍无关。在这种情况下，当将实现从一种语言移植到另一种语言时，请努力确保输入位精确，并使用相对误差(而不是复制粘贴和比较打印输出)来比较输出。
(请注意，由于通带的对称性，b的所有其他元素都是 0。这可用于减少所需的翻牌数量！