Matlab中的滤波函数filter的代码

jtll521

impulse 发表于 2012-11-7 15:04

                               
登录/注册后可看大图

我的经过检验没问题，把你的代码贴出来，后面的不对，我估计是你的初始值赋值有问题，另外注意这个函数与 ...

首先，谢谢您及时回复，下面的是我自己写的filter 主要就是想实现matlab 里面的filter作用。

ord 是得到的滤波器系数的个数，
a,b,分别是滤波器的两个系数，np是数据个数，
x是输入数据
y是结果
void Filter(int ord, double *a, double *b, int np, double *x, double *y)
{
    int i,j;
    y[0]=b[0]*x[0];
    for (i=1;i<ord;i++)
    {
        y=0.0;
        for (j=0;j<i+1;j++)
            y=y+b[j]*x[i-j];
        for (j=0;j<i;j++)
            y=y-a[j+1]*y[i-j-1];
    }
    for (i=ord;i<np;i++)
    {
        y=0.0;
        for (j=0;j<ord;j++)
            y=y+b[j]*x[i-j];
        for (j=0;j<ord;j++)
            y=y-a[j+1]*y[i-j-1];
    }
}

下面这一段是我调用的你贴出的代码，我输入了20个数据，在x里面。
&vecData[0]  &vecData1[0]这两个是用matlab 获得的滤波器系数，个数为11个
      double* y1 = NULL;
      double* zi = new double[20];
       memset(zi,0,20*sizeof(double));
   filter(x, 20, &vecData[0], 11,  &vecData1[0], 11, zi, 20, y1, hhg);

用你这个代码和我上面写的代码得到的结果是一样的，但是和matlab 里面出来的结果后面的部分不太一样，
希望您看一下谢谢。
另外我也是东南大学的，是生物医学工程的。看来您应该是师兄了。呵呵

impulse

jtll521 发表于 2012-11-7 15:26

                               
登录/注册后可看大图

首先，谢谢您及时回复，下面的是我自己写的filter 主要就是想实现matlab 里面的filter作用。

ord 是得 ...

zi显然错了，zi在matlab中的个数好像max（a，b）-1，我的函数里和b的长度一样，我在函数说明里还特意注明了这一点。
我是东大3系的，05年毕业的。

jtll521

impulse 发表于 2012-11-7 16:00

                               
登录/注册后可看大图

zi显然错了，zi在matlab中的个数好像max（a，b）-1，我的函数里和b的长度一样，我在函数说明里还特意注明 ...

嗯，谢谢您回复。但是我运行之后确实有问题啊，你帮我看一下可以吗？纠结了很久了这个问题

[B, A]=butter(5,[0.002,0.016])
X=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
format long;
filter(B,A,X)
下面这段代码是我用上面的语句生成的滤波器系数，相当于是个5阶的带通滤波器吧，但是不管用你的函数还是自己写的filfer
计算出来的和用上面的matlab出来的结果后面几个不一样，前面的一半左右是完全正确的。您能有时间帮我看一下吗？这个纠结了几天了，一直以为可能是最后的filter有问题，但是又不知道在哪里出了问题。谢谢
    double y[20]={0};
    double x[20]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};

    double ddd[11] = {0.000000026659795524309, 0, -0.000000133298977621545, 0,
    0.000000266597955243090 , 0, -0.000000266597955243090, 0, 0.0000000133298977621545,0
    -0.000000026659795524309    };    //导入用matlab生成的滤波器系数

    double dd1[11] = {1.0,-9.8561167706614,43.7168612494972,-114.9140602054998,  198.2401376041890,
         -234.5190776628390,192.6758570892651, -108.5537178561892, 40.1380766719147, -8.7952830647879,0.8673229451114 };
     
    vector<double> vecData1;
    vector<double> vecData;
    for (k=0; k<11; k++)
    {
      vecData.push_back(ddd[k]);
    }

    for (k=0; k<11; k++)
    {
        vecData1.push_back(dd1[k]);
    }
       int hhg;
       double* y1 = NULL;
       double* zi = new double[11];
       memset(zi,0,11*sizeof(double));
       filter(x, 20, &vecData[0], 11,  &vecData1[0], 11, zi, 11, y1, hhg);

impulse

jtll521 发表于 2012-11-8 01:47

                               
登录/注册后可看大图

嗯，谢谢您回复。但是我运行之后确实有问题啊，你帮我看一下可以吗？纠结了很久了这个问题

=butter(5 ...

这两天实在忙，周日我来看看

jtll521

impulse 发表于 2012-11-8 09:59

                               
登录/注册后可看大图

这两天实在忙，周日我来看看

嗯，好的。谢谢了

redplum

多谢了，大家的讨论很有启发

impulse

jtll521 发表于 2012-11-9 01:32

                               
登录/注册后可看大图

嗯，好的。谢谢了

结果没错呀，和matlab结果误差为E-15。

代码如下：

1. double* y=NULL;
   
2. double x[20]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
   
3. double b[11] ={0.000000026659795524309, 0, -0.000000133298977621545, 0, 0.000000266597955243090, 0, -0.000000266597955243090, 0, 0.0000000133298977621545, 0, -0.000000026659795524309}; //导入用matlab生成的滤波器系数
   
4. double a[11] ={1.0, -9.8561167706614, 43.7168612494972, -114.9140602054998, 198.2401376041890, -234.5190776628390, 192.6758570892651, -108.5537178561892, 40.1380766719147, -8.7952830647879, 0.8673229451114};
   
5. int len_y=0;
   
6. double* z=new double[11];
   
7. memset((double*)z, 0, sizeof(double)*11);
   
8. filter(x, 20, b, 11, a, 11, z, 11, y, len_y);
   
9. //saveData("Y.txt", y, len_y ,2);
   
10. delete[] y;

复制代码

jtll521

impulse 发表于 2012-11-9 22:34

                               
登录/注册后可看大图

结果没错呀，和matlab结果误差为E-15。

代码如下：

嗯，我又仔细检测了一下代码，我里面有个地方写错了。
不好意思，浪费您的时间了，真的很感谢您每次都能及时回答。
您好像也做过信号方面的东西，我最近也在做这个，希望能和您多交流。
还有很多不懂的东西。

不知道您处理信号时有没有滤波，我想既然写了filter估计也有滤波函数，我
现在用的是网上找的一个低通滤波器，是c++写的，但是我感觉效果不是很好，
不知道您有没有可以用的滤波器函数，能发一份给我吗？十分感谢
这是毕设的内容，我对信号方面的不太懂，我是学图像的，希望指导一下，当然是
在您有时间的情况下，谢谢啦。

songwenshuai

你好 能给我一个完整的 C工程吗?我不太理解参数都是什么，谢谢

songwenshuai

// filter.cpp: 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void filter(double *x, int len_x, double *coeff_b, int len_b, double *coeff_a, int len_a, double* zi, int len_zi, double* &filter_x, int& len_filter_x);
int main(int argc, char *argv[])
{
    int i;

    // 第一个参数
    double x[20] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20 };
    // 第二个参数
    int  len_x = 20;
    // 第三个参数
    double coeff_b[18] = { -0.000000,0.000000,0.000000,0.000000,0.000001,0.000020,0.000153,0.000539,0.000957,0.000888,0.000430,0.000105,0.000012,0.000001,0.000000,0.000000,-0.000000,0.000000 };
    // 第四个参数
    int len_b = 18;
    // 第五个参数
    double coeff_a[18] = { 1.000000,-8.507642,38.059387,-116.740268,272.044806,-507.528390,781.139627,-1009.730407,1107.339849,-1034.816858,823.630706,-555.293117,313.568811,-145.476731,53.711375,-14.938188,2.813064,-0.272917 };
    // 第六个参数
    int len_a = 18;
    // 第七个参数
    double* zi = new double[100];
    memset((double*)zi, 0, sizeof(double) * 100);
    // 第八个参数
    int len_zi = 100;
    // 第九个参数
    double* filter_x = NULL;
    // 第十个参数
    int len_filter_x = 0;
    filter(x, len_x, coeff_b, len_b, coeff_a, len_a, zi, len_zi, filter_x, len_filter_x);
    printf("Uint Impulse Response \n");
    for (i = 0; i < len_zi; i++) {
            printf("%f\n",zi[i]);
    }
    system("pause");
    return 0;
}

/*函数使用说明：
1、
2、对于len_a=1的情况，函数内部补充数组a的个数与b相同；
3、函数与matlab的filter函数区别在于zi的个数不同，本函数为len_b，matlab则为len_b-1
filter(x, 20, &vecData[0], 11,  &vecData1[0], 11, zi, 20, y1, hhg);
*/
/**************************************************************************************************
* 函 数     : void filter(b,a,m,n,x,len,px,py)
* 调 用     : filter
* 参 数     :
*             x                 一个数组 double x[20]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}; 长度是
*             len_x             x 的长度
*             coeff_b           系数
*             len_b             滤波器分母系数个数
*             coeff_a           系数
*                      len_a             滤波器分子系数个数
*             zi                zi 既是输入参数(initial conditions),又是输出参数(final conditions),因此,zi需要赋值,首次调用赋0；
*             len_zi            zi 长度
*             filter_x
*             len_filter_x
*************************************************************************************************/
void filter(double *x, int len_x, double *coeff_b, int len_b, double *coeff_a, int len_a, double* zi, int len_zi, double* &filter_x, int& len_filter_x)
{

    len_filter_x = len_x;
    filter_x = new double[len_x];
    if (len_a == 1)
    {
        for (int m = 0; m<len_x; m++)
        {
            filter_x[m] = coeff_b[0] * x[m] + zi[0];
            for (int i = 1; i<len_b; i++)
            {
                zi[i - 1] = coeff_b[i] * x[m] + zi[i];//-coeff_a[i]*filter_x[m];
            }
        }
    }
    else
    {
        for (int m = 0; m<len_x; m++)
        {
            filter_x[m] = coeff_b[0] * x[m] + zi[0];
            for (int i = 1; i<len_b; i++)
            {
                zi[i - 1] = coeff_b[i] * x[m] + zi[i] - coeff_a[i] * filter_x[m];
            }
        }
    }
}

songwenshuai

我看不太懂 您能帮我下吗？

impulse

songwenshuai 发表于 2017-10-30 08:53
我看不太懂 您能帮我下吗？

22楼有函数使用代码啊，注意zi长度