9
$\begingroup$

Is anybody familiar with Gustafson's algorithm for minimizing transients in forward backward filtering [1]? I'm trying to implement it and my first guess was to check Matlab's filtfilt.m, since they are referencing the paper. In the Matlab function also a linear equation system is solved in order to find initial conditions zi that minimize the startup transients, but the relationship between reference and code is not obvious to me. The only lines of code regarding the minimization are (nfilt is the length of the coefficient vectors):

zi = ( eye(nfilt-1) - [-a(2:nfilt), [eye(nfilt-2); zeros(1,nfilt-2)]] ) \...
 ( b(2:nfilt) - b(1)*a(2:nfilt) );

Can anybody point me in the right direction on how those lines relate to the algorithm described in Gustafson's article?

[1] Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992.

Improve this question
$\endgroup$
8
  • $\begingroup$ initial states of any IIR filter should be zero at the beginning of the forward-filtering pass and should be zero at the beginning of the backward-filtering pass. in both passes, the signal file (or buffer) being filtered will get longer by the apparent length of the IIR (how long it takes for the output to decay closely enough to zero that you can choose to cut off the rest of the decay). $\endgroup$
    – robert bristow-johnson
    Feb 24, 2017 at 22:09
  • $\begingroup$ In the paper, the author claims that the forward and backward filters are different in their state space representation. Can you explain why? $\endgroup$
    – Maxtron
    May 14, 2018 at 18:08
  • $\begingroup$ well, as i understand the use of filtfilt() i cannot see why. i have not read the Gustafson paper (i'm not IEEE and can't get it for free, anyone who has a copy is welcome to email me a .pdf of it). in using the concept of filtfilt, one can do it to an entire file of samples (for me it would be an audio or sound file, like a .wav) first forward filter the sound with it zero-padded on the end by as long as you expect the impulse response of the forward filter to be. that lengthens the file, but the output gets virtually to zero. then run the resulting file through the filter backwards. $\endgroup$
    – robert bristow-johnson
    May 14, 2018 at 19:36
  • $\begingroup$ there's another use, based on a paper by Powell and Chau that does filtfilt in real time by breaking the input into blocks of samples, zero-padding each block, filtering the blocks backward but keeping the "tails" flipping it back around to the forward direction and overlap-adding. Powell-Chau did not do this, but i think this is a good application of Truncated IIR filters, so you know when the decaying block output ends. $\endgroup$
    – robert bristow-johnson
    May 14, 2018 at 19:43
  • 1
    $\begingroup$ @robertbristow-johnson: I ran into this copy of the Gustafsson paper. $\endgroup$
    – djvg
    Oct 10, 2018 at 6:32

2 Answers 2

Reset to default
5
$\begingroup$

For anyone who is interested, i coincidentally found a paper describing the method implemented in matlab's filtfilt.m. A link to the paper is attached. At least to my understanding matlab's filtfilt.m doesn't implement the Gustafson algorithm.

Sadovsky, P.; Bartusek, K: Optimisation of the Transient Response of a Digital Filter, Radioengineering Vol. 9, No. 2, 2000

Improve this answer
$\endgroup$
1
  • $\begingroup$ Also see the scipy documentation for lfilter_zi, which is the default in scipy.signal.filtfilt for determining the initial conditions,as you can see in the source. In this case 'odd' padding is used by default, but it can use Gustafsson's method as an option (have a look at the definition for _filtfilt_gust in the source). $\endgroup$
    – djvg
    Oct 10, 2018 at 8:49
0
$\begingroup$

The zi = (...)\(...) line in the OP's question determines the initial state of the filter. I believe the exact same approach is used by Python's scipy. According to the scipy docs:

A linear filter with order m has a state space representation (A, B, C, D), for which the output y of the filter can be expressed as:

z(n+1) = Az(n) + Bx(n)

y(n) = Cz(n) + Dx(n)

where z(n) is a vector of length m, A has shape (m, m), B has shape (m, 1), C has shape (1, m) and D has shape (1, 1) (assuming x(n) is a scalar). lfilter_zi solves:

zi = A*zi + B

In other words, it finds the initial condition for which the response to an input of all ones is a constant.

(my emphasis)

Note that filtfilt calculates the initial state zi, scales it by the first sample value, then passes it on to filter, which actually applies it (docs).

A basic example of how this initial state zi can be applied in a filter, using the state space representation, is provided here.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.