the data with the values at the end points. On
both the forward and backward passes, the initial condition of the
filter is found by using `lfilter_zi` and scaling it by the end point of
the extended data.

When `method` is "gust", Gustafsson's method [1]_ is used.  Initial
conditions are chosen for the forward and backward passes so that the
forward-backward filter gives the same result as the backward-forward
filter.

The option to use Gustaffson's method was added in scipy version 0.16.0.

References
----------
.. [1] F. Gustaffson, "Determining the initial states in forward-backward
       filtering", Transactions on Signal Processing, Vol. 46, pp. 988-992,
       1996.

Examples
--------
The examples will use several functions from `scipy.signal`.

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt

First we create a one second signal that is the sum of two pure sine
waves, with frequencies 5 Hz and 250 Hz, sampled at 2000 Hz.

>>> t = np.linspace(0, 1.0, 2001)
>>> xlow = np.sin(2 * np.pi * 5 * t)
>>> xhigh = np.sin(2 * np.pi * 250 * t)
>>> x = xlow + xhigh

Now create a lowpass Butterworth filter with a cutoff of 0.125 times
the Nyquist frequency, or 125 Hz, and apply it to ``x`` with `filtfilt`.
The result should be approximately ``xlow``, with no phase shift.

>>> b, a = signal.butter(8, 0.125)
>>> y = signal.filtfilt(b, a, x, padlen=150)
>>> np.abs(y - xlow).max()
9.1086182074789912e-06

We get a fairly clean result for this artificial example because
the odd extension is exact, and with the moderately long padding,
the filter's transients have dissipated by the time the actual data
is reached.  In general, transient effects at the edges are
unavoidable.

The following example demonstrates the option ``method="gust"``.

First, create a filter.

>>> b, a = signal.ellip(4, 0.01, 120, 0.125)  # Filter to be applied.

`sig` is a random input signal to be filtered.

>>> rng = np.random.default_rng()
>>> n = 60
>>> sig = rng.standard_normal(n)**3 + 3*rng.standard_normal(n).cumsum()

Apply `filtfilt` to `sig`, once using the Gustafsson method, and
once using padding, and plot the results for comparison.

>>> fgust = signal.filtfilt(b, a, sig, method="gust")
>>> fpad = signal.filtfilt(b, a, sig, padlen=50)
>>> plt.plot(sig, 'k-', label='input')
>>> plt.plot(fgust, 'b-', linewidth=4, label='gust')
>>> plt.plot(fpad, 'c-', linewidth=1.5, label='pad')
>>> plt.legend(loc='best')
>>> plt.show()

The `irlen` argument can be used to improve the performance
of Gustafsson's method.

Estimate the impulse response length of the filter.

>>> z, p, k = signal.tf2zpk(b, a)
>>> eps = 1e-9
>>> r = np.max(np.abs(p))
>>> approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r)))
>>> approx_impulse_len
137

Apply the filter to a longer signal, with and without the `irlen`
argument.  The difference between `y1` and `y2` is small.  For long
signals, using `irlen` gives a significant performance improvement.

>>> x = rng.standard_normal(4000)
>>> y1 = signal.filtfilt(b, a, x, method='gust')
>>> y2 = signal.filtfilt(b, a, x, method='gust', irlen=approx_impulse_len)
>>> print(np.max(np.abs(y1 - y2)))
2.875334415008979e-10

)