C99相当于MATLAB“filter”？

为什么MATLAB和C版本会产生不同的结果？

MATLAB：

[B_coeffs, A_coeffs ] = butter(4, 100/(16000/2), 'high'); state = zeros( 4, 1 ); input = zeros( 64,1 ); for i=1:64 input(i)=i; end [filtered_output, state] = filter( B_coeffs, A_coeffs, input, state ); 

C：

 int main(...) { for(int test=0; teststate[0], z1 = hpf->state[1], z2 = hpf->state[2], z3 = hpf->state[3]; unsigned int N = 64; unsigned int i = 0; for(i=0;iB_coeffs[0] * Xi + z0; z0 = hpf->B_coeffs[1] * Xi + z1 - hpf->A_coeffs[1] * Yi; z1 = hpf->B_coeffs[2] * Xi + z2 - hpf->A_coeffs[2] * Yi; z2 = hpf->B_coeffs[3] * Xi + z3 - hpf->A_coeffs[3] * Yi; z3 = hpf->B_coeffs[4] * Xi - hpf->A_coeffs[4] * Yi; Yout[i] = (float) Yi; } hpf->state[0] = z0; hpf->state[1] = z1; hpf->state[2] = z2; hpf->state[3] = z3; return; } 

哪里

 typedef struct { float A_coeffs[5]; float B_coeffs[5]; float state[4]; } t_high_pass_filter; void high_pass_filter_init( t_high_pass_filter* hpf) { hpf->A_coeffs[0] = 1.0000; hpf->A_coeffs[1] = -3.8974; hpf->A_coeffs[2] = 5.6974; hpf->A_coeffs[3] = -3.7025; hpf->A_coeffs[4] = 0.9025; hpf->B_coeffs[0] = 0.9500; hpf->B_coeffs[1] = -3.7999; hpf->B_coeffs[2] = 5.6999; hpf->B_coeffs[3] = -3.7999; hpf->B_coeffs[4] = 0.9500; hpf->state[0] = 0.0; hpf->state[1] = 0.0; hpf->state[2] = 0.0; hpf->state[3] = 0.0; } 

**产出是：**

 MATLAB: C: ---------------------------- 0.9500 0.9500 1.8025 1.8026 2.5625 2.5631 3.2350 3.2369 3.8247 3.8292 4.3360 4.3460 4.7736 4.7930 5.1416 5.1767 5.4442 5.5035 5.6854 5.7807 5.8691 6.0156 5.9991 6.2162 6.0788 6.3909 6.1119 6.5487 6.1016 6.6989 6.0511 6.8514 5.9637 7.0167 5.8421 7.2057 5.6894 7.4298 5.5083 7.7009 5.3013 8.0314 5.0710 8.4342 4.8199 8.9225 4.5501 9.5101 4.2640 10.2110 3.9637 11.0399 3.6511 12.0115 3.3281 13.1412 2.9965 14.4443 2.6582 15.9368 2.3146 17.6347 1.9674 19.5543 1.6180 21.7122 1.2677 24.1250 0.9179 26.8095 0.5698 29.7829 0.2245 33.0621 -0.1169 36.6643 -0.4535 40.6067 -0.7842 44.9066 -1.1084 49.5812 -1.4251 54.6477 -1.7336 60.1232 -2.0333 66.0249 -2.3236 72.3697 -2.6039 79.1746 -2.8738 86.4562 -3.1327 94.2313 -3.3804 102.5161 -3.6164 111.3270 -3.8405 120.6801 -4.0524 130.5911 -4.2520 141.0756 -4.4390 152.1490 -4.6134 163.8264 -4.7750 176.1225 -4.9239 189.0520 -5.0600 202.6291 -5.1832 216.8677 -5.2938 231.7814 -5.3917 247.3837 -5.4771 263.6875 -5.5501 280.7055 -5.6108 298.4500 

---

前几个值是相同的（或类似的），但随后它们发散。 此外，在第一次迭代后，filter状态完全不同。

我究竟做错了什么？

Answer 1

在第二次编辑之后，您的问题变得清晰： 混乱的行为 。

第一次，您似乎刚刚将MATLAB命令窗口中的系数复制到C函数中。 但是，MATLAB的format似乎设置为short ，因为C函数中的系数四舍五入到小数点后4位。 这个舍入（以及第一次在C函数中使用float ）是你的问题。

这是我这次做的：

1. 复制你的MATLAB脚本
2. 复制您的C代码，并将其转换为MATLAB MEX格式，以便更轻松地比较事物
3. 调整C代码，使其接受
   1. 什么都没有（在这种情况下，它使用像以前一样的“内置”， 圆形版本
   2. 与MATLAB脚本中使用的系数相同（带有附加数字）
4. 运行脚本并进行比较。

结论：您看到对初始值的敏感度非常高。

TL; DR ：此代码：

 clc [B_coeffs, A_coeffs] = butter(4, 100/(16000/2), 'high'); state = zeros(4, 1); input = 1:64; % MATLAB version [filtered_output0, state0] = filter(B_coeffs, A_coeffs, input, state); mex filter_test.c % MEX, using built-in constants (of which only the first few digits are equal) [filtered_output1, state1] = filter_test([],[], input, state, 0); % MEX, using the exact same constants [filtered_output2, state2] = filter_test(B_coeffs, A_coeffs, input, state, 1); % Compare! [filtered_output0.' filtered_output1.' filtered_output2.'] [state0 state1 state2] 

filter_test.c包含：

 #include  #include "mex.h" #define N ( 64u) #define C ( 5u) #define S (C-1u) /* helper struct for HPF */ typedef struct { double A_coeffs[C]; double B_coeffs[C]; double state[S]; } t_high_pass_filter; /* "Default" values (note that these are ROUNDED to 4 digits only) void high_pass_filter_init(t_high_pass_filter* hpf) { hpf->A_coeffs[0] = 1.0000; hpf->A_coeffs[1] = -3.8974; hpf->A_coeffs[2] = 5.6974; hpf->A_coeffs[3] = -3.7025; hpf->A_coeffs[4] = 0.9025; hpf->B_coeffs[0] = 0.9500; hpf->B_coeffs[1] = -3.7999; hpf->B_coeffs[2] = 5.6999; hpf->B_coeffs[3] = -3.7999; hpf->B_coeffs[4] = 0.9500; hpf->state[0] = 0.0; hpf->state[1] = 0.0; hpf->state[2] = 0.0; hpf->state[3] = 0.0; } /* the actual filter */ void high_pass_filter_do(t_high_pass_filter* hpf, const double *Xin, double *Yout) { double Xi, Yi; double z0 = hpf->state[0], z1 = hpf->state[1], z2 = hpf->state[2], z3 = hpf->state[3]; unsigned int i = 0u; for(i=0; iB_coeffs[0] * Xi + z0; z0 = hpf->B_coeffs[1] * Xi + z1 - hpf->A_coeffs[1] * Yi; z1 = hpf->B_coeffs[2] * Xi + z2 - hpf->A_coeffs[2] * Yi; z2 = hpf->B_coeffs[3] * Xi + z3 - hpf->A_coeffs[3] * Yi; z3 = hpf->B_coeffs[4] * Xi - hpf->A_coeffs[4] * Yi; Yout[i] = Yi; } hpf->state[0] = z0; hpf->state[1] = z1; hpf->state[2] = z2; hpf->state[3] = z3; } /* Wrapper between MATLAB MEX and filter function */ void filter(const double *B_coeffs, const double *A_coeffs, const double *input, const double *state, double *filtered_output, double *state_output) { unsigned int i = 0u; t_high_pass_filter hpf; /* Use built-in defaults when coefficients * have not been passed explicitly */ if (B_coeffs == NULL) { high_pass_filter_init(&hpf); } /* Otherwise, use the coefficients on the arguments */ else { for (i=0u; i 

给出这个输出：

 % MATLAB C with rounding C without rounding %=============================================================================== filtered values: 9.499817826393004e-001 9.500000000000000e-001 9.499817826393004e-001 1.802482117980145e+000 1.802630000000000e+000 1.802482117980145e+000 2.562533610562189e+000 2.563140162000000e+000 2.562533610562189e+000 ... (58 more values) ... ... -5.477082906502112e+000 2.637900416547026e+002 -5.477082906502112e+000 -5.550054712403821e+000 2.808186934967214e+002 -5.550054712403821e+000 -5.610782439991141e+000 2.985749768301545e+002 -5.610782439991141e+000 states after filter run: -6.740826417997180e+001 2.553210086702404e+002 -6.740826417997180e+001 2.006572641218511e+002 -7.151404729138806e+002 2.006572641218511e+002 -1.991114894348111e+002 6.686913808328562e+002 -1.991114894348111e+002 6.586237103693912e+001 -2.086639165892145e+002 6.586237103693912e+001 

下一次，确保你在复制粘贴之前使用format long e 🙂

Answer 2

当我运行你的代码时，我得到了这个：

 MATLAB C ================================ 0.9500 0.950000 1.8026 1.802630 2.5631 2.563139 ... (58 more values) ... 263.7900 263.687500 280.8187 280.705475 298.5750 298.450012 

但是，当我改变C代码使用double而不是float ，我得到了这个：

 MATLAB C ================================ 0.9500 0.950000 1.8026 1.802630 2.5631 2.563140 ... (58 more values) ... 263.7900 263.790042 280.8187 280.818693 298.5750 298.574977 

所以，@ Richard是正确的：

1. 您在MATLAB方面做错了，我无法重现您的原始值
2. 您的C代码是正确的，但与MATLAB版本不同，因为MATLAB使用double而不是float 。