matlab filter函数C实现

https://groups.google.com/d/msg/comp.soft-sys.matlab/TpI55V__870/OecRMp3TCcIJ

/* 
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 <[email protected]> 
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 
[email protected] 
*/ 

#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 

/*matlab test
>> a = [1 0 0]
>> b = [1, 2, 3]
>> x = [1, 2, 3, 4, 8, 7, 6, 5, 1, 4, 2, 3]
>> y = filter(b, a, x)
y =
1     4    10    16    25    35    44    38    29    21    13    19
*/

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 */