C99 equivalent to MATLAB “filter”?

Question

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

where

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

**The outputs are: **

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.

What am I doing wrong?

Answer 1

When I run your code, I get this:

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

But, when I change the C code to use double instead of float, I get this:

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 

So, @Richard is correct:

1. You did something wrong on the MATLAB side, I can't reproduce your original values
2. Your C code is correct, but differs from the MATLAB version because MATLAB uses double instead of float.

Answer 2

After your second edit, it is becoming clear what your problem is: chaotic behavior.

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.

Here's what I did this time:

1. copy your MATLAB script
2. Copy your C code, and cast it into MATLAB MEX format to compare things more easily
3. Adjust the C code such that it accepts either
   1. nothing (in which case it uses the "built-in", rounded versions like before
   2. The same coefficients as used in your MATLAB script (with additional digits)
4. Run scripts and compare.

Conclusions: you're seeing a very high sensitivity to initial values.

TL;DR: this code:

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]

Where filter_test.c contains:

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

gives this 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

Next time, make sure you use format long e before copy-pasting like that :)