1. 首頁
  2. 問答
  3. C99等價於MATLAB「過濾器」？

C99等價於MATLAB「過濾器」？

2016-11-18 164 views
1

爲什麼MATLAB和C版本產生不同的結果？C99等價於MATLAB「過濾器」？

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

前幾個值是相同的（或相似），但然後他們分歧。另外，第一次迭代後，濾波器狀態完全不同。

我在做什麼錯？

來源

2016-11-18 Danijel

+0

默認情況下，Matlab使用雙精度。如果您將C代碼更改爲在任何地方使用雙精度，會發生什麼？ – Richard

+0

不能輕易改變雙精度，使用一些FFT庫等，所以需要我很多時間。 – Danijel

回答

3

經過第二次編輯後，問題變得清晰：chaotic behavior。

第一次，你似乎剛把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 <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.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之前複製粘貼:)

來源

2016-11-21 11:21:49

+0

就是這樣，最後。謝謝。 MATLAB'format long'會做到這一點。算法對於濾波器係數的微小差異非常敏感，這是非常令人意想不到的，從來沒有預料到！再次感謝。 – Danijel

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。

來源

2016-11-18 13:50:36

+0

謝謝，這有助於很多。但是，仍然沒有找到爲什麼我的MATLAB不能正常工作。有些事情出錯了......： -/ – Danijel

+0

@Danijel鍵入'過濾器'來查看是否有陰影功能或者其他東西......它應該是內置的 –

+0

@Danijel請注意，'filter'重新調整輸入if 'A（1）≠1'，這不是在你的C版本中完成的...但是，你的'A（1）= 1'，所以......這不應該是問題 –

相關問題

 相關問題