Initialize filter so that filtered output begins with initial value of the input

90 vues (au cours des 30 derniers jours)
Peter Nagy
Peter Nagy le 12 Mar 2017
I would like to implement a recursive exponential filter with the following function: filter([alpha 0],[1 alpha-1],data); This is equivalent of the following equation: y(n)=y(n-1) * (1-alpha) + alpha * x(n) where x is the input and y is the output. With the above function I get alpha * x(1) for y(1) which means that y(0) is assumed to be zero. I expected that with the following function filter([alpha 0],[1 alpha-1],data,data(1)); I would get x(1) for y(1), but it is not the case. What is the problem? How can I make sure that y(1)=x(1)? Thanks in advance.
  1 commentaire
Christoph F.
Christoph F. le 25 Avr 2017
y(n) = y(n-1)*(1-a) + a*x(n)
And you want x(1) = y(1) so:
x(1) = y(0)*(1-a) + a*x(1)
(1-a)*x(1) = (1-a)*y(0)
y(0) = x(1)
You need to set the initial condition y(0) of your filter to x(1). Note that the exact values of the initial conditions depend on the filter structure, so if you're using something other than direct form 1, the calculation of the initial values for the delay line is somewhat different.
Example, step response to x=1, with transition
a = 0.25; y=0; x=1; for ii=1:40; y = y*(1-a) + a*x; yout(ii)=y; end;plot( yout);
Example with y(0) set to x(1):
a = 0.25; y=1; x=1; for ii=1:40; y = y*(1-a) + a*x; yout(ii)=y; end;plot( yout);
I don't think using filtfilt is the answer, since it effectively doubles the filter order and changes the filters amplitude response significantly. Also, it is not possible to use forward-backward-filtering if real-time output is required.

Connectez-vous pour commenter.

Réponses (3)

Peter Nagy
Peter Nagy le 13 Mar 2017
I also decided to use filtfilt which works. I just wanted to figure out why filter didn't work the way I originally expected. For the sake of the community I'm glad that I managed to figure it out. At the very bottom of the help of filter filter delay (z) is defined. For y(1) the equation is: y(1)=b1*x(1)+z1(0), for the recursive exponential filter I wanted to implement b(1)=alpha, therefore y(1)=alpha*x(1)+z1(0) If I want y(1) to be equal to x(1): x(1)=alpha*x(1)+z1(0) solving for z1(0) yields z1(0)=x1*(1-alpha)
If I use the following command: filter([alpha 0],[1 alpha-1],data,data(1)*(1-alpha));
I get the x(1) for y(1).
Peter

Star Strider
Star Strider le 12 Mar 2017
I always use the filtfilt (link) function for filtering because it avoids this problem and also has a maximally flat phase response, so that all filters, regardless of design, have the phase response that a hardware Bessel filter would have.
If you want to use the filter (link) function, see the section of the documentation on Filter Data in Sections (link) for details on setting the initial conditions.

Christoph F.
Christoph F. le 25 Avr 2017
Oops. That was supposed to be an answer, not a comment. So I'm posting it again as an answer.
y(n) = y(n-1)*(1-a) + a*x(n)
And you want x(1) = y(1) so:
x(1) = y(0)*(1-a) + a*x(1)
(1-a)*x(1) = (1-a)*y(0)
y(0) = x(1)
You need to set the initial condition y(0) of your filter to x(1). Note that the exact values of the initial conditions depend on the filter structure, so if you're using something other than direct form 1, the calculation of the initial values for the delay line is somewhat different.
Example, step response to x=1, with transition
a = 0.25; y=0; x=1; for ii=1:40; y = y*(1-a) + a*x; yout(ii)=y; end;plot( yout);
Example with y(0) set to x(1):
a = 0.25; y=1; x=1; for ii=1:40; y = y*(1-a) + a*x; yout(ii)=y; end;plot( yout);
I don't think using filtfilt is the answer, since it effectively doubles the filter order and changes the filters amplitude response significantly. Also, it is not possible to use forward-backward-filtering if real-time output is required.

Community Treasure Hunt

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

Start Hunting!

Translated by