filtfilt C语言实现，可直接运行验证

#include 
#include 
#include 
#include 

/*滤波总体流程：
1、确定延时系数zi
2、对数据进行镜像，消除边缘效应
3、正向滤波，反向滤波
4、将数据从滤波后的数据中截取出来*/

//filter函数
void filter（const float* x float* y int xlen float* a float* b int nfilt float* zi）//nfilt为系数数组长度，zi为延时系数
{
#define EPS 0.000001
	float tmp;
	int i j;


	//判断a[0]是否为1，若不是，归一化下
	if （（*a - 1.0 > EPS） || （*a - 1.0 < -EPS））
	{
		tmp = *a;
		for （i = 0; i < nfilt; i++）
		{
			b[i] /= tmp;
			a[i] /= tmp;
		}
	}


	memset（y 0 xlen * sizeof（float））;//将y清零，以双浮点为单位

    //此处赋值a[0]=0可以直接一步循环到位，不需要讲b[0]*x[0]单独计算
	a[0] = 0.0;
	for （i = 0; i < xlen; i++）
	{
		for （j = 0; i >= j && j < nfilt; j++）
		{
			y[i] += （b[j] * x[i - j] - a[j] * y[i - j]）;
		}
		//若滤波器阶数大于2，需加上延时系数
		if （zi&&i < nfilt - 1） y[i] += zi[i];
	}
	a[0] = 1.0;//还原滤波系数
}




//矩阵乘法 m为a的行数，n为a的列数数，k为b的行数，第一个矩阵列数必须要和第二个矩阵的行数相同，输出结果存放在c中
void trmul（float *a float *b float *c int m int n int k）
{
	int i j l u;
	for （i = 0; i <= m - 1; i++）
		for （j = 0; j <= k - 1; j++）
		{
			u = i * k + j; c[u] = 0.0;
			for （l = 0; l <= n - 1; l++）
				c[u] = c[u] + a[i*n + l] * b[l*k + j];
		}
	return;
}


//求逆矩阵，当返回值为0时成功，a变为逆矩阵
int rinv（float *a int n）//逆矩阵
{
	int *is *js i j k l u v;
	float d p;
	is = （int *）malloc（n * sizeof（int））;
	js = （int *）malloc（n * sizeof（int））;
	for （k = 0; k <= n - 1; k++）
	{
		d = 0.0;
		for （i = k; i <= n - 1; i++）
			for （j = k; j <= n - 1; j++）
			{
				l = i * n + j; p = fabs（a[l]）;
				if （p > d） { d = p; is[k] = i; js[k] = j; }
			}
		if （is[k] != k）
			for （j = 0; j <= n - 1; j++）
			{
				u = k * n + j; v = is[k] * n + j;
				p = a[u]; a[u] = a[v]; a[v] = p;
			}
		if （js[k] != k）
			for （i = 0; i <= n - 1; i++）
			{
				u = i * n + k; v = i * n + js[k];
				p = a[u]; a[u] = a[v]; a[v] = p;
			}
		l = k * n + k;
		a[l] = 1.0 / a[l];
		for （j = 0; j <= n - 1; j++）
			if （j != k）
			{
				u = k * n + j; a[u] = a[u] * a[l];
			}
		for （i = 0; i <= n - 1; i++）
			if （i != k）
				for （j = 0; j <= n - 1; j++）
					if （j != k）
					{
						u = i * n + j;
						a[u] = a[u] - a[i*n + k] * a[k*n + j];
					}
		for （i = 0; i <= n - 1; i++）
			if （i != k）
			{
				u = i * n + k; a[u] = -a[u] * a[l];
			}
	}
	for （k = n - 1; k >= 0; k--）
	{
		if （js[k] != k）
			for （j = 0; j <= n - 1; j++）
			{
				u = k * n + j; v = js[k] * n + j;
				p = a[u]; a[u] = a[v]; a[v] = p;
			}
		if （is[k] != k）
			for （i = 0; i <= n - 1; i++）
			{
				u = i * n + k; v