C99 эквивалент MATLAB "фильтр"?

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

С:

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

После вашего второго редактирования становится ясно, в чем ваша проблема: chaotic behavior.

Первый раз, похоже, вы только что скопировали коэффициенты из командного окна MATLAB в свою функцию C. Тем не менее, format MATLAB, по-видимому, установлен в short, так как коэффициенты в функции C равны округленным до 4 знаков после запятой. Это закругление (а также использование float в функции C в первый раз) является вашей проблемой.

Вот что я сделал на этот раз:

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 перед тем, как копировать такие файлы :)

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 

Итак, @ Ричард прав:

1. Вы сделали что-то неправильно на стороне MATLAB, я не могу воспроизвести ваши исходные значения
2. Ваш код C верен, но отличается от версии MATLAB, поскольку MATLAB использует double вместо float.