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

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.

[Christoph F.]

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.

Answer 1

[Peter Nagy]

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

Answer 2

[Star Strider]

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.

Answer 3

[Christoph F.]

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.