Created
October 6, 2014 15:48
-
-
Save jfosorio/3a4ff3eb150d665fe73b to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
def sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None): | |
""" | |
A forward-backward filter. | |
This function applies a linear filter twice, once forward | |
and once backwards. The combined filter has linear phase. | |
Before applying the filter, the function can pad the data along the | |
given axis in one of three ways: odd, even or constant. The odd | |
and even extensions have the corresponding symmetry about the end point | |
of the data. The constant extension extends the data with the values | |
at end points. On both the forward and backwards passes, the | |
initial condition of the filter is found by using `lfilter_zi` and | |
scaling it by the end point of the extended data. | |
Parameters | |
---------- | |
sos : array_like | |
Array of second-order filter coefficients, must have shape | |
``(n_sections, 6)``. Each row corresponds to a second-order | |
section, with the first three columns providing the numerator | |
coefficients and the last three providing the denominator | |
coefficients. | |
x : array_like | |
The array of data to be filtered. | |
axis : int, optional | |
The axis of `x` to which the filter is applied. | |
Default is -1. | |
padtype : str or None, optional | |
Must be 'odd', 'even', 'constant', or None. This determines the | |
type of extension to use for the padded signal to which the filter | |
is applied. If `padtype` is None, no padding is used. The default | |
is 'odd'. | |
padlen : int or None, optional | |
The number of elements by which to extend `x` at both ends of | |
`axis` before applying the filter. This value must be less than | |
`x.shape[axis]-1`. `padlen=0` implies no padding. | |
The default value is 3*max(len(a),len(b)). | |
Returns | |
------- | |
y : ndarray | |
The filtered output, an array of type numpy.float64 with the same | |
shape as `x`. | |
See Also | |
-------- | |
soslfilter_zi, soslfilter | |
Examples | |
-------- | |
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 rate, or 125 Hz, and apply it to x with filtfilt. The | |
result should be approximately xlow, with no phase shift. | |
>>> from scipy import signal | |
>>> sos = signal.butter(8, 0.125, output='sos') | |
>>> y = signal.sosfiltfilt(sos, 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. | |
""" | |
if padtype not in ['even', 'odd', 'constant', None]: | |
raise ValueError(("Unknown value '%s' given to padtype. padtype must " | |
"be 'even', 'odd', 'constant', or None.") % | |
padtype) | |
# b = np.atleast_1d(b) | |
# a = np.atleast_1d(a) | |
x = np.asarray(x) | |
ntaps = sos.shape[0]*2 | |
if padtype is None: | |
padlen = 0 | |
if padlen is None: | |
# Original padding; preserved for backwards compatibility. | |
edge = ntaps * 3 | |
else: | |
edge = padlen | |
# x's 'axis' dimension must be bigger than edge. | |
if x.shape[axis] <= edge: | |
raise ValueError("The length of the input vector x must be at least " | |
"padlen, which is %d." % edge) | |
if padtype is not None and edge > 0: | |
# Make an extension of length `edge` at each | |
# end of the input array. | |
if padtype == 'even': | |
ext = even_ext(x, edge, axis=axis) | |
elif padtype == 'odd': | |
ext = odd_ext(x, edge, axis=axis) | |
else: | |
ext = const_ext(x, edge, axis=axis) | |
else: | |
ext = x | |
# Get the steady state of the filter's step response. | |
zi = sosfilt_zi(sos) | |
# Reshape zi and create x0 so that zi*x0 broadcasts | |
# to the correct value for the 'zi' keyword argument | |
# to lfilter. | |
zi_shape = [1] * x.ndim | |
zi_shape[axis] = zi.size | |
zi = np.reshape(zi, zi_shape) | |
x0 = axis_slice(ext, stop=1, axis=axis) | |
zix0 = zi * x0 | |
zix0 = reshape(zix0, (sos.shape[0], 2 )) | |
# Forward filter. | |
(y, zf) = sosfilt(sos, ext, axis=axis, zi=zix0) | |
# Backward filter. | |
# Create y0 so zi*y0 broadcasts appropriately. | |
y0 = axis_slice(y, start=-1, axis=axis) | |
ziy0 = zi * y0 | |
ziy0 = reshape(zix0, (sos.shape[0], 2 )) | |
(y, zf) = sosfilt(sos, axis_reverse(y, axis=axis), axis=axis, zi=ziy0) | |
# Reverse y. | |
y = axis_reverse(y, axis=axis) | |
if edge > 0: | |
# Slice the actual signal from the extended signal. | |
y = axis_slice(y, start=edge, stop=-edge, axis=axis) | |
return y |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment