Implementation of the matlab function filtfilt in C language

41 visualizaciones (últimos 30 días)
afef
afef el 18 de Ag. de 2017
Editada: Oybek Amonov el 5 de Nov. de 2020
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]);
}

Respuestas (2)

Michele Giugliano
Michele Giugliano el 20 de Mzo. de 2018
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

Jan
Jan el 21 de Mzo. de 2018
Editada: Jan el 21 de Mzo. de 2018
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 comentario
Oybek Amonov
Oybek Amonov el 5 de Nov. de 2020
Editada: Oybek Amonov el 5 de Nov. de 2020
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.

Iniciar sesión para comentar.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by