Implementation of the matlab function filtfilt in C language

Question

afef 2017-8-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/353172-implementation-of-the-matlab-function-filtfilt-in-c-language

编辑： Oybek Amonov 2020-11-5

I used in Matlab the function filtfilt for my Low pass fir filter :

d1=designfilt('lowpassfir','PassbandFrequency',0.45,'StopbandFrequency',0.5,'Pass bandRipple',3,'StopbandAttenuation',60,'DesignMethod','equiripple'); a = filtfilt(d1,deriv1);

And I tried to implement this function in C language I don't get the same output that I get with Matlab and I don't understand the problem. Can anyone help me please?

This is the code:

    #define ORDER 65
      #define NP 2560
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];
      }
      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];
      }
}
main()
{
float x[NP]= { -84.786,0.000,0.000,0.000,0.000,-0.781....};
float y[NP],a[ORDER+1],b[ORDER+1];
int i,j;
for (i=0;i<NP;i++)
   printf("%f\n",x[i]);
/* test coef from
[b,a]=butter(3,30/500); in MATLAB
*/
b[0]=-0.0056;
b[1]=-0.0202;
b[2]=-0.0341;
b[3]=-0.0303;
b[4]=-0.0058;
b[5]= 0.0157;
b[6]=0.0117;
b[7]=-0.0081;
b[8]=-0.0126;
b[9]=0.0038;
b[10]=0.0133;
b[11]=-0.000;
b[12]=-0.0139;
b[13]=-0.0027;
b[14]=0.0147;
b[15]=0.0065;
b[16]=-0.0150;
b[17]=-0.0110;
b[18]=0.0147;
b[19]=0.0163;
b[20]=-0.0134;
b[21]=-0.0228;
b[22]=0.0107;
b[23]=0.0308;
b[24]=-0.0057;
b[25]=-0.0412;
b[26]=-0.0030;
b[27]=0.0561;
b[28]=0.0197;
b[29]=-0.0833;
b[30]=-0.0617;
b[31]=0.1726;
b[32]=0.4244;
b[33]=0.4244;
b[34]=0.1726;
b[35]=-0.0617;
b[36]=-0.0833;
b[37]=0.0197;
b[38]=0.0561;
b[39]=-0.0030;
b[40]=-0.0412;
b[41]=-0.0057;
b[42]=0.0308;
b[43]=0.0107;
b[44]=-0.0228;
b[45]=-0.0134;
b[46]=0.0163;
b[47]=0.0147;
b[48]=-0.0110;
b[49]=-0.0150;
b[50]=0.0065;
b[51]=0.0147;
b[52]=-0.0027;
b[53]=-0.0139;
b[54]=-0.0005;
b[55]=0.0133;
b[56]=0.0038;
b[57]=-0.0126;
b[58]=-0.0081;
b[59]=0.0117;
b[60]=0.0157;
b[61]=-0.0058;
b[62]=-0.0303;
b[63]=-0.0341;
b[64]=-0.0202;
b[65]=-0.0056;
a[0]=0;
a[1]=0;
a[2]=0;
a[3]=0;
a[4]=0;
a[5]=0;
a[6]=0;
a[7]=0;
a[8]=0;
a[9]=0;
a[10]=0;
a[11]=0;
a[12]=0;
a[13]=0;
a[14]=0;
a[15]=0;
a[16]=0;
a[17]=0;
a[18]=0;
a[19]=0;
a[20]=0;
a[21]=0;
a[22]=0;
a[23]=0;
a[24]=0;
a[25]=0;
a[26]=0;
a[27]=0;
a[28]=0;
a[29]=0;
a[30]=0;
a[31]=0;
a[32]=0;
a[33]=0;
a[34]=0;
a[35]=0;
a[36]=0;
a[37]=0;
a[38]=0;
a[39]=0;
a[40]=0;
a[41]=0;
a[42]=0;
a[43]=0;
a[44]=0;
a[45]=0;
a[46]=0;
a[47]=0;
a[48]=0;
a[49]=0;
a[50]=0;
a[51]=0;
a[52]=0;
a[53]=0;
a[54]=0;
a[55]=0;
a[56]=0;
a[57]=0;
a[58]=0;
a[59]=0;
a[60]=0;
a[61]=0;
a[62]=0;
a[63]=0;
a[64]=0;
a[65]=0;
filter(ORDER,a,b,NP,x,y);
for (i=0;i<NP;i++)
{ x[i]=y[NP-i-1];}
filter(ORDER,a,b,NP,x,y);
for (i=0;i<NP;i++)
    { x[i]=y[NP-i-1];}
    for (i=0;i<NP;i++)
    { y[i]=x[i];}
     for (i=0;i<NP;i++)
    printf("%f\n",y[i]);
    }

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Michele Giugliano 2018-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/353172-implementation-of-the-matlab-function-filtfilt-in-c-language#answer_310993

Not sure it (still) helps you: the MATLAB filtfilt.m implementation has an internal method to reduce edge effects, while initialising the internal forward and backward filtering. Have a look at filtfilt.m and at the paper cited therein:

https://pdfs.semanticscholar.org/f98b/09d9ef5b8cedfc7d1de0138f9d13b9b87b58.pdf

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Jan 2018-3-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/353172-implementation-of-the-matlab-function-filtfilt-in-c-language#answer_310998

编辑：Jan 2018-3-21

在 MATLAB Online 中打开

See https://www.mathworks.com/matlabcentral/fileexchange/32261-filterm for an implementation of FILTFILT. Here the treatment of the initial and final sequence is process in Matlab and only the filtering in implemented in C.

This is not correct:

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

a and b must be applied simultaneously and you need to store the intermediate states. See this M-version of the code:

function [Y, z] = emulateFILTER(b, a, X, z)
b    = b ./ a(1);
a    = a ./ a(1);
n    = length(a);
z(n) = 0;
Y    = zeros(size(X));
for m = 1:length(Y)
   Xm   = X(m);
   Y(m) = b(1) * Xm + z(1);
   Ym   = Y(m);
   for i = 2:n
      z(i - 1) = b(i) * Xm + z(i) - a(i) * Ym;
   end
end
z = z(1:n - 1);

Your hard-coded filter parameters have 3 significant digits only. This will cause inaccurate results in addition.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Oybek Amonov 2020-11-5

编辑：Oybek Amonov 2020-11-5

Is z the actual result here? I do not understand what z is here. As I understand it, Y(m) is the result of forward filtering, so I though z would be the actual zero-phase filtering result.

请先登录，再进行评论。