filter.m to C - algorithm problem

[yllerena]

I'm trying to create my own C implementation of matlab filter.m funtion to down-sample a signal around 2 kHz to speed up autocorrelation process. But I have some inconsistences on results. That's way I need some help on this.

To test the veracity of my implementation and under assuntion of:
y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na)

In filte.c (attached) file are the test set and result set and also the matlab result set for comparition)
==================================
#include <stdio.h>
#include <stdlib.h>

void filter(double *b, int ord, double a, double *x, int np, double *y) {
y[0]=b[0]*x[0];
int n,k;
for ( n=0; n<ord; n++) {

for (k=0; k<n+1; k++)
y[n] += b[k]*x[n-k];

for (k=0; k<n; k++)
y[n] -= a*y[n-k-1];
}

for (n=ord; n<np; n++) {

for (k=0; k<ord; k++)
y[n] += b[k]*x[n-k];

for (k=0; k<ord; k++)
y[n] -= a*y[n-k-1];
}
}

void main() {
double b[24] = {0.000000000000000, 0.000000333087156, 0.000001496873275, 0.000003958627788, 0.000008120976691, 0.000014289216468, 0.000022641741288, 0.000033205011502,
0.000045834367448, 0.000060201817729, 0.000075791709912, 0.000091904930112, 0.000107671983786, 0.000122074992086, 0.000133978306167, 0.000142167106464,
0.000145393026160, 0.000142425528942, 0.000132107491488, 0.000113413201241, 0.000085506789149, 0.000047798983375, -0.000000000000000,-0.000057833614664};
double x[72] = {0.018035888671875, -0.019561767578125, -0.021911621093750, -0.022125244140625, -0.022583007812500, -0.023925781250000, -0.024719238281250, -0.025207519531250,
-0.026062011718750, -0.027343750000000, -0.028350830078125, -0.030029296875000, -0.030975341796875, -0.031372070312500, -0.032012939453125, -0.032806396484375,
-0.033569335937500, -0.035156250000000, -0.035186767578125, -0.035461425781250, -0.036071777343750, -0.037292480468750, -0.037322998046875, -0.038085937500000,
-0.039215087890625, -0.039245605468750, -0.039672851562500, -0.039916992187500, -0.040374755859375, -0.039184570312500, -0.038635253906250, -0.039123535156250,
-0.040679931640625, -0.040679931640625, -0.039825439453125, -0.039672851562500, -0.039276123046875, -0.039093017578125, -0.038635253906250, -0.038543701171875,
-0.038879394531250, -0.039001464843750, -0.038238525390625, -0.037597656250000, -0.036773681640625, -0.036376953125000, -0.035766601562500, -0.034881591796875,
-0.034667968750000, -0.034088134765625, -0.033111572265625, -0.032562255859375, -0.031585693359375, -0.030548095703125, -0.029327392578125, -0.029144287109375,
-0.028320312500000, -0.027709960937500, -0.026672363281250, -0.026062011718750, -0.025482177734375, -0.025115966796875, -0.024230957031250, -0.023498535156250,
-0.023101806640625, -0.022247314453125, -0.021301269531250, -0.020904541015625, -0.020080566406250, -0.018554687500000, -0.018188476562500, -0.017486572265625};

double *y;
y = (double *) malloc(72 * sizeof(double));

filter(b,24,1,x,72,y);
int i;
for (i=0; i<72; i++) {
printf("%.15f ",y[i]);
}

}

/* Result
0.000000000000000 0.000000006007523 0.000000014474143 0.000000014335737 -0.000000005954683 -0.000000057385410 -0.000000149938358 -0.000000291793276
-0.000000488407714 -0.000000741937453 -0.000001050908125 -0.000001409802937 -0.000001809000273 -0.000002234576123 -0.000002668108635 -0.000003086976954
-0.000003464981861 -0.000003773101643 -0.000003980694477 -0.000004056219216 -0.000003968281081 -0.000003687129048 -0.000003186343097 -0.000002443992560
0.000000815087300 -0.000001254812358 -0.000001366742322 -0.000001307769569 -0.000001235503455 -0.000001253497965 -0.000001256025278 -0.000001244436353
-0.000001286030533 -0.000001418951286 -0.000001591368623 -0.000001904559844 -0.000002231430220 -0.000002538024524 -0.000002870260838 -0.000003215156027
-0.000003540022139 -0.000003926858378 -0.000004115407794 -0.000004192924874 -0.000004156223620 -0.000004035167979 -0.000003625825633 -0.000003061956385
-0.000002331996243 0.000000987735378 -0.000001061928367 -0.000001152228361 -0.000001110064527 -0.000000878691016 -0.000000726927706 -0.000000638802812
-0.000000696999764 -0.000000729513200 -0.000000760205157 -0.000000869949092 -0.000001091533397 -0.000001339949343 -0.000001538558821 -0.000001797326394
-0.000002141919099 -0.000002498693914 -0.000002859401337 -0.000003030116043 -0.000003069301805 -0.000003037576816 -0.000002916745842 -0.000002482098645
*/

/*
Matlab results
0.000000000000000, 0.000038728700701, 0.000212468158825, 0.000669952090429, 0.001603869555242, 0.003239789754726, 0.005821697041210, 0.009595712388848,
0.014788450415421, 0.021584557421138, 0.030107784610860, 0.040399973744600, 0.052400321296955, 0.065929790887386, 0.080678650284441, 0.096198113468023,
0.111898106348139, 0.127052412074315, 0.140811853415091, 0.152223263872654, 0.160255890238859, 0.163831872700465, 0.161860570084496, 0.153278230462050,
0.151654132385857, 0.150283032845994, 0.148819248112176, 0.147502321307337, 0.146569989138610, 0.146029414393046, 0.146376158011450, 0.146483290451214,
0.146723609633208, 0.147979917324537, 0.149170865285480, 0.150130608131086, 0.151282454133807, 0.152053641843408, 0.152245396299540, 0.151768031227638,
0.150816292130522, 0.149814912956444, 0.148571612894952, 0.147797061447995, 0.146939840157426, 0.145917255327909, 0.144987225361455, 0.143622055161194,
0.141875587930011, 0.140309739204180, 0.138772362906510, 0.137599749478392, 0.136113778987935, 0.135812276863310, 0.136298453733806, 0.136624098564743,
0.136895645597126, 0.137684533116548, 0.138179364639825, 0.139324399180288, 0.140602894249594, 0.141566998378014, 0.141782363390037, 0.141445743997878,
0.140566082677064, 0.139567963532260, 0.138825406303969, 0.137300249608882, 0.135515699671258, 0.134191265772710, 0.132725082619978, 0.131163927795268
*/
=====================================
Best regards

[Thang Ngo]

/*
http://mechatronics.ece.usu.edu/yqchen/filter.c/FILTER.C
FILTER.C
An ANSI C implementation of MATLAB FILTER.M (built-in)
Written by Chen Yangquan <ele...@nus.edu.sg>
1998-11-11

Example of FILTER.C that translate from filter.m matlab
test main and comments are written by NgoTranDucThang - Danang University of Technology, VietNam
tngo...@gmail.com
*/

#include<stdio.h>
#define ORDER 2 /* ORDER is the number of element a[] or b[] minus 1;*/
#define NP 11/* NP is the number of output or input filter minus 1;*/


void filter(int,float *,float *,int,float *,float *);
/*
Y = filter(B,A,X)
a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

- a(2)*y(n-1) - ... - a(na+1)*y(n-na)

*/
//this program is corresponding to command "y = filter(b,1,x)" in matlab
void main(void){

float a[3]={1,0,0};
float b[3]={1,2,3};
float x[12]={1,2,3,4,8,7,6,5,1,4,2,3};
float y[12];

filter(ORDER,a,b,NP,x,y);

return;
}
void filter(int ord, float *a, float *b, int np, float *x, float *y)
{
int i,j;

y[0]=b[0]*x[0];

for (i=1;i<ord+1;i++)
{
y[i]=0.0;
for (j=0;j<i+1;j++)
y[i]=y[i]+b[j]*x[i-j];
for (j=0;j<i;j++)
y[i]=y[i]-a[j+1]*y[i-j-1];
}
/* end of initial part */
for (i=ord+1;i<np+1;i++)
{
y[i]=0.0;
for (j=0;j<ord+1;j++)
y[i]=y[i]+b[j]*x[i-j];
for (j=0;j<ord;j++)
y[i]=y[i]-a[j+1]*y[i-j-1];
}
return;
} /* end of filter */

Cause it's has no response so I just comment for others

[Madi]

Hi thanks for the code.. but just a bi confused:

if I want to user for N = 3 with Butterworth filter and got this coefficient:

A:
1
-5.63705390875799
13.4684474638276
-17.4440131274997
12.9142339871711
-5.18269945149946
0.881578059262132

B:
2.93906303448239e-05
0
-8.81718910344717e-05
0
8.81718910344717e-05
0
-2.93906303448239e-05

From the code given, where can I put this value??

and for:

#define NP 11 /* NP is the number of output or input filter minus 1;*/

now what number of NP should I put? It is should according to my number of sample data?

Thank You

[Jan Simon]

Hi yllerena,

The direct form II approach is much faster. Here a simple M-implementation, which creates exactly the same results as Matlab's FILTER:
http://www.mathworks.com/matlabcentral/answers/9900#answer_13623

An equivalent C-implementation:

void CoreDoubleN(double *X, mwSize MX, mwSize NX, double *a, double *b,
mwSize order, double *Z, double *Y)
{
// Direct form II transposed method for general filter length.
// Implemented as time domain difference equations.
// INPUT:
// X: Double array. Operation happens of 1st dimension.
// MX: Number of elements in the 1st dimension
// NX: Number of columns, considers mutliple dimensions.
// a, b: Double vector of filter parameters. Both have nParam elements.
// The first element of a is 1.
// Z: DOUBLE array, initial conditions.
// order: Number of filter parameters, order of filter + 1.
// OUTPUT:
// Z: DOUBLE array, final conditions.
// Y: Double array, allocated by the caller.

double Xi, Yi;
mwSize i, j, R;

i = 0;
while (NX--) { // Next slice
R = i + MX; // End of the column
while (i < R) {
Xi = X[i]; // Get signal
Yi = b[0] * Xi + Z[0]; // Filtered value
for (j = 1; j < order; j++) { // Update conditions
Z[j - 1] = b[j] * Xi + Z[j] - a[j] * Yi;
}
Z[order - 1] = b[order] * Xi - a[order] * Yi;

Y[i++] = Yi; // Write to output
}
Z += order; // Next condition vector
}

To my surprise this is up to 3 times faster than Matlab's FILTER. See:
http://www.mathworks.de/matlabcentral/fileexchange/32261-filterm

Kind regards, Jan

[TJ Sartain]

This is very intriguing; especially about being faster than the actual MATLAB. I am, however, having trouble implementing this in Obj-C. I'm using C-style arrays as part of an object called Matrix; hence the "->items" notation. Here's my translation:

- (Matrix *) filter: (NSArray *)BA z: (NSArray *)zi
{
Matrix *Y = [Matrix zeros:rows cols:1];
Matrix *B = [[Matrix alloc] initWithNSArray:BA[0]];
Matrix *A = [[Matrix alloc] initWithNSArray:BA[1]];
Matrix *Z = [[Matrix alloc] initWithNSArray:zi];
long double Xi, Yi;
int i, j, R;
int order = [B height] ;
int MX = rows;
int NX = cols;

i = 0;
while (NX--) {

R = i + MX;

while (i < R)
{
Xi = items[i];
Yi = B->items[0] * Xi + Z->items[0];

for (j = 1; j < order; j++)
{

Z->items[j - 1] = B->items[j] * Xi + Z->items[j] - A->items[j] * Yi;
}
Z->items[order - 1] = B->items[order] * Xi - A->items[order] * Yi;

Y->items[i++] = Yi;
}
}
return Y;
}

You'll notice I'm missing the Z += order. That's because I don't know what that is.

Sometimes I get malloc error (which is weird that it's sometimes). But when it does execute all the way, I don't get MATLAB-matching values. Anyone see any glaring errors.