C99是否等同于MATLAB“过滤器"? [英] C99 equivalent to MATLAB "filter"?
本文介绍了C99是否等同于MATLAB“过滤器"?的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
为什么MATLAB和C版本会产生不同的结果?
Why do the MATLAB and C versions produce different results?
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; test<64;test++)
Xin[test]=test+1;
...
high_pass_filter_init(...)
high_pass_filter_do(...)
}
// Do the filtering
void high_pass_filter_do( t_high_pass_filter* hpf, float *Xin, float *Yout )
{
double Xi, Yi;
double z0 = hpf->state[0],
z1 = hpf->state[1],
z2 = hpf->state[2],
z3 = hpf->state[3];
unsigned int N = 64;
unsigned int i = 0;
for(i=0;i<N;i++)
{
Xi = Xin[i];
Yi = hpf->B_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
前几个值相同(或相似)，但随后却有所不同.同样，第一次迭代后，滤波器状态完全不同.
The first few values are the same (or similar), but then they diverge. Also, the filter state is totally different after first iteration.
我在做什么错了?
推荐答案
第二次编辑后，您的问题已经很清楚了:混乱行为.
After your second edit, it is becoming clear what your problem is: chaotic behavior.
第一次，您似乎刚刚将系数从MATLAB命令窗口复制到了C函数中.但是，MATLAB的format似乎已设置为short，因为C函数中的系数四舍五入到小数点后四位.这种舍入(以及第一次在C函数中使用float)是您的问题.
First time round, you seem to have just copied the coefficients from the MATLAB command window into your C function. However, MATLAB's format appears to have been set to short, as the coefficients in the C function are rounded to 4 decimal places. This rounding (as well as using float in the C function the first time around) is your problem.
这是我这次所做的:
- 复制您的MATLAB脚本
- 复制您的C代码，并将其转换为MATLAB MEX格式，以更轻松地进行比较
- 调整C代码，使其可以接受
- copy your MATLAB script
- Copy your C code, and cast it into MATLAB MEX format to compare things more easily
- Adjust the C code such that it accepts either
- 什么都没有(在这种情况下，它使用像以前一样的内置"，四舍五入版本
- 与MATLAB脚本中使用的系数相同(带有附加数字)
结论:您看到对初始值的灵敏度非常高.
Conclusions: you're seeing a very high sensitivity to initial values.
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 <stdio.h>
#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; i<N; ++i)
{
Xi = Xin[i];
Yi = hpf->B_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<C; ++i) {
hpf.B_coeffs[i] = B_coeffs[i];
hpf.A_coeffs[i] = A_coeffs[i];
}
for (i=0u; i<S; ++i)
hpf.state[i] = state[i];
}
/* Apply filter function */
high_pass_filter_do(&hpf,
input,
filtered_output);
/* Assign output state explicitly */
for (i=0u; i<S; ++i)
state_output[i] = hpf.state[i];
}
/* MATLAB/MEX Gateway function */
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
unsigned int i = 0u;
double *B_coeffs = mxGetPr(prhs[0]),
*A_coeffs = mxGetPr(prhs[1]),
*input = mxGetPr(prhs[2]),
*state = mxGetPr(prhs[3]);
int flag = *mxGetPr(prhs[4]);
double *filtered_output,
*state_output;
/* filtered_output, state */
plhs[0] = mxCreateDoubleMatrix(1, N, mxREAL);
plhs[1] = mxCreateDoubleMatrix(S, 1, mxREAL);
filtered_output = mxGetPr(plhs[0]);
state_output = mxGetPr(plhs[1]);
if (flag == 0)
{
filter(NULL,
NULL,
input,
state,
filtered_output,
state_output);
}
else
{
filter(B_coeffs,
A_coeffs,
input,
state,
filtered_output,
state_output);
}
}
提供以下输出:
% 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 :)
Next time, make sure you use format long e before copy-pasting like that :)
这篇关于C99是否等同于MATLAB“过滤器"?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
