C99はMATLAB "filter"に相当しますか？

Question

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; 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 

---

最初のいくつかの値は同じ（または類似していますが）、それらは乖離します。また、最初の反復の後では、フィルタの状態は全く異なります。

私は間違っていますか？私はあなたのコードを実行すると

Answer 1

2回目の編集後、問題がどのようなものであるかが明確になっています。chaotic behavior。

MATLABコマンドウィンドウの係数をC関数にコピーしたばかりです。ただし、MATLABのformatは、shortに設定されているように見えます.C関数の係数はで、小数点以下4桁にはを四捨五入しています。この丸め（最初にC関数でfloatを使用するだけでなく）も問題です。ここで

は、私は、この時間をやったことだ：

1. は
2. あなたのMATLABスクリプトをコピーし、あなたのCコードをコピーし、
3. より簡単に物事を比較するために、MATLAB MEX形式にそれをキャストしているCのコードは、このような調整それは、「ビルトイン」を使用した場合のいずれかで
   1. 何も（受け入れない、は前
   ようバージョンを丸め0
   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が含まれています

% 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.086702404e+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を使用してください：

#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); 
    } 
} 

はこの出力を提供しますそのようにコピー貼り付ける前に：

Answer 2

、私はこれを取得：だから、

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 

@：私は代わりにfloatのdoubleを使用するようにCコードを変更すると、

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 

をしかし、私はこの取得しますリチャードは正しいです：

1. 何か MATLAB側で間違っていますが、元の値を再現できません
2. MATLABはfloatの代わりにdoubleを使用するため、Cコードは正確ですがMATLABのバージョンとは異なります。