signal.signaltools

Module Contents

Functions

_valfrommode(mode)
_bvalfromboundary(boundary)
_inputs_swap_needed(mode,shape1,shape2) If in ‘valid’ mode, returns whether or not the input arrays need to be
correlate(in1,in2,mode=”full”,method=”auto”) r
_centered(arr,newshape)
fftconvolve(in1,in2,mode=”full”) Convolve two N-dimensional arrays using FFT.
_numeric_arrays(arrays,kinds=”buifc”) See if a list of arrays are all numeric.
_prod(iterable) Product of a list of numbers.
_fftconv_faster(x,h,mode) See if using fftconvolve or _correlateND is faster. The boolean value
_reverse_and_conj(x) Reverse array x in all dimensions and perform the complex conjugate
_np_conv_ok(volume,kernel,mode) See if numpy supports convolution of volume and kernel (i.e. both are
_timeit_fast(stmt=”pass”,setup=”pass”,repeat=3) Returns the time the statement/function took, in seconds.
choose_conv_method(in1,in2,mode=”full”,measure=False) Find the fastest convolution/correlation method.
convolve(in1,in2,mode=”full”,method=”auto”) Convolve two N-dimensional arrays.
order_filter(a,domain,rank) Perform an order filter on an N-dimensional array.
medfilt(volume,kernel_size=None) Perform a median filter on an N-dimensional array.
wiener(im,mysize=None,noise=None) Perform a Wiener filter on an N-dimensional array.
convolve2d(in1,in2,mode=”full”,boundary=”fill”,fillvalue=0) Convolve two 2-dimensional arrays.
correlate2d(in1,in2,mode=”full”,boundary=”fill”,fillvalue=0) Cross-correlate two 2-dimensional arrays.
medfilt2d(input,kernel_size=3) Median filter a 2-dimensional array.
lfilter(b,a,x,axis=None,zi=None) Filter data along one-dimension with an IIR or FIR filter.
lfiltic(b,a,y,x=None) Construct initial conditions for lfilter given input and output vectors.
deconvolve(signal,divisor) Deconvolves divisor out of signal using inverse filtering.
hilbert(x,N=None,axis=None) Compute the analytic signal, using the Hilbert transform.
hilbert2(x,N=None) Compute the ‘2-D’ analytic signal of x
cmplx_sort(p) Sort roots based on magnitude.
unique_roots(p,tol=0.001,rtype=”min”) Determine unique roots and their multiplicities from a list of roots.
invres(r,p,k,tol=0.001,rtype=”avg”) Compute b(s) and a(s) from partial fraction expansion.
residue(b,a,tol=0.001,rtype=”avg”) Compute partial-fraction expansion of b(s) / a(s).
residuez(b,a,tol=0.001,rtype=”avg”) Compute partial-fraction expansion of b(z) / a(z).
invresz(r,p,k,tol=0.001,rtype=”avg”) Compute b(z) and a(z) from partial fraction expansion.
resample(x,num,t=None,axis=0,window=None) Resample x to num samples using Fourier method along the given axis.
resample_poly(x,up,down,axis=0,window=tuple) Resample x along the given axis using polyphase filtering.
vectorstrength(events,period) Determine the vector strength of the events corresponding to the given
detrend(data,axis=None,type=”linear”,bp=0) Remove linear trend along axis from data.
lfilter_zi(b,a) Construct initial conditions for lfilter for step response steady-state.
sosfilt_zi(sos) Construct initial conditions for sosfilt for step response steady-state.
_filtfilt_gust(b,a,x,axis=None,irlen=None) Forward-backward IIR filter that uses Gustafsson’s method.
filtfilt(b,a,x,axis=None,padtype=”odd”,padlen=None,method=”pad”,irlen=None) Apply a digital filter forward and backward to a signal.
_validate_pad(padtype,padlen,x,axis,ntaps) Helper to validate padding for filtfilt
sosfilt(sos,x,axis=None,zi=None) Filter data along one dimension using cascaded second-order sections.
sosfiltfilt(sos,x,axis=None,padtype=”odd”,padlen=None) A forward-backward digital filter using cascaded second-order sections.
decimate(x,q,n=None,ftype=”iir”,axis=None,zero_phase=True) Downsample the signal after applying an anti-aliasing filter.
_valfrommode(mode)
_bvalfromboundary(boundary)
_inputs_swap_needed(mode, shape1, shape2)

If in ‘valid’ mode, returns whether or not the input arrays need to be swapped depending on whether shape1 is at least as large as shape2 in every dimension.

This is important for some of the correlation and convolution implementations in this module, where the larger array input needs to come before the smaller array input when operating in this mode.

Note that if the mode provided is not ‘valid’, False is immediately returned.

correlate(in1, in2, mode="full", method="auto")

r Cross-correlate two N-dimensional arrays.

Cross-correlate in1 and in2, with the output size determined by the mode argument.

in1 : array_like
First input.
in2 : array_like
Second input. Should have the same number of dimensions as in1.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear cross-correlation of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
method : str {‘auto’, ‘direct’, ‘fft’}, optional

A string indicating which method to use to calculate the correlation.

direct
The correlation is determined directly from sums, the definition of correlation.
fft
The Fast Fourier Transform is used to perform the correlation more quickly (only available for numerical arrays.)
auto

Automatically chooses direct or Fourier method based on an estimate of which is faster (default). See convolve Notes for more detail.

New in version 0.19.0.

correlate : array
An N-dimensional array containing a subset of the discrete linear cross-correlation of in1 with in2.

choose_conv_method : contains more documentation on method.

The correlation z of two d-dimensional arrays x and y is defined as:

z[...,k,...] = sum[..., i_l, ...] x[..., i_l,...] * conj(y[..., i_l - k,...])

This way, if x and y are 1-D arrays and z = correlate(x, y, 'full') then

for

where is the length of x, , and is 0 when m is outside the range of y.

method='fft' only works for numerical arrays as it relies on fftconvolve. In certain cases (i.e., arrays of objects or when rounding integers can lose precision), method='direct' is always used.

Implement a matched filter using cross-correlation, to recover a signal that has passed through a noisy channel.

>>> from scipy import signal
>>> sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 128)
>>> sig_noise = sig + np.random.randn(len(sig))
>>> corr = signal.correlate(sig_noise, np.ones(128), mode='same') / 128
>>> import matplotlib.pyplot as plt
>>> clock = np.arange(64, len(sig), 128)
>>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True)
>>> ax_orig.plot(sig)
>>> ax_orig.plot(clock, sig[clock], 'ro')
>>> ax_orig.set_title('Original signal')
>>> ax_noise.plot(sig_noise)
>>> ax_noise.set_title('Signal with noise')
>>> ax_corr.plot(corr)
>>> ax_corr.plot(clock, corr[clock], 'ro')
>>> ax_corr.axhline(0.5, ls=':')
>>> ax_corr.set_title('Cross-correlated with rectangular pulse')
>>> ax_orig.margins(0, 0.1)
>>> fig.tight_layout()
>>> fig.show()
_centered(arr, newshape)
fftconvolve(in1, in2, mode="full")

Convolve two N-dimensional arrays using FFT.

Convolve in1 and in2 using the fast Fourier transform method, with the output size determined by the mode argument.

This is generally much faster than convolve for large arrays (n > ~500), but can be slower when only a few output values are needed, and can only output float arrays (int or object array inputs will be cast to float).

As of v0.19, convolve automatically chooses this method or the direct method based on an estimation of which is faster.

in1 : array_like
First input.
in2 : array_like
Second input. Should have the same number of dimensions as in1. If operating in ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
out : array
An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

Autocorrelation of white noise is an impulse.

>>> from scipy import signal
>>> sig = np.random.randn(1000)
>>> autocorr = signal.fftconvolve(sig, sig[::-1], mode='full')
>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1)
>>> ax_orig.plot(sig)
>>> ax_orig.set_title('White noise')
>>> ax_mag.plot(np.arange(-len(sig)+1,len(sig)), autocorr)
>>> ax_mag.set_title('Autocorrelation')
>>> fig.tight_layout()
>>> fig.show()

Gaussian blur implemented using FFT convolution. Notice the dark borders around the image, due to the zero-padding beyond its boundaries. The convolve2d function allows for other types of image boundaries, but is far slower.

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> kernel = np.outer(signal.gaussian(70, 8), signal.gaussian(70, 8))
>>> blurred = signal.fftconvolve(face, kernel, mode='same')
>>> fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1,
...                                                      figsize=(6, 15))
>>> ax_orig.imshow(face, cmap='gray')
>>> ax_orig.set_title('Original')
>>> ax_orig.set_axis_off()
>>> ax_kernel.imshow(kernel, cmap='gray')
>>> ax_kernel.set_title('Gaussian kernel')
>>> ax_kernel.set_axis_off()
>>> ax_blurred.imshow(blurred, cmap='gray')
>>> ax_blurred.set_title('Blurred')
>>> ax_blurred.set_axis_off()
>>> fig.show()
_numeric_arrays(arrays, kinds="buifc")

See if a list of arrays are all numeric.

ndarrays : array or list of arrays
arrays to check if numeric.
numeric_kinds : string-like
The dtypes of the arrays to be checked. If the dtype.kind of the ndarrays are not in this string the function returns False and otherwise returns True.
_prod(iterable)

Product of a list of numbers. Faster than np.prod for short lists like array shapes.

_fftconv_faster(x, h, mode)

See if using fftconvolve or _correlateND is faster. The boolean value returned depends on the sizes and shapes of the input values.

The big O ratios were found to hold across different machines, which makes sense as it’s the ratio that matters (the effective speed of the computer is found in both big O constants). Regardless, this had been tuned on an early 2015 MacBook Pro with 8GB RAM and an Intel i5 processor.

_reverse_and_conj(x)

Reverse array x in all dimensions and perform the complex conjugate

_np_conv_ok(volume, kernel, mode)

See if numpy supports convolution of volume and kernel (i.e. both are 1D ndarrays and of the appropriate shape). Numpy’s ‘same’ mode uses the size of the larger input, while Scipy’s uses the size of the first input.

_timeit_fast(stmt="pass", setup="pass", repeat=3)

Returns the time the statement/function took, in seconds.

Faster, less precise version of IPython’s timeit. stmt can be a statement written as a string or a callable.

Will do only 1 loop (like IPython’s timeit) with no repetitions (unlike IPython) for very slow functions. For fast functions, only does enough loops to take 5 ms, which seems to produce similar results (on Windows at least), and avoids doing an extraneous cycle that isn’t measured.

choose_conv_method(in1, in2, mode="full", measure=False)

Find the fastest convolution/correlation method.

This primarily exists to be called during the method='auto' option in convolve and correlate, but can also be used when performing many convolutions of the same input shapes and dtypes, determining which method to use for all of them, either to avoid the overhead of the ‘auto’ option or to use accurate real-world measurements.

in1 : array_like
The first argument passed into the convolution function.
in2 : array_like
The second argument passed into the convolution function.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
measure : bool, optional
If True, run and time the convolution of in1 and in2 with both methods and return the fastest. If False (default), predict the fastest method using precomputed values.
method : str
A string indicating which convolution method is fastest, either ‘direct’ or ‘fft’
times : dict, optional
A dictionary containing the times (in seconds) needed for each method. This value is only returned if measure=True.

convolve correlate

For large n, measure=False is accurate and can quickly determine the fastest method to perform the convolution. However, this is not as accurate for small n (when any dimension in the input or output is small).

In practice, we found that this function estimates the faster method up to a multiplicative factor of 5 (i.e., the estimated method is at most 5 times slower than the fastest method). The estimation values were tuned on an early 2015 MacBook Pro with 8GB RAM but we found that the prediction held fairly accurately across different machines.

If measure=True, time the convolutions. Because this function uses fftconvolve, an error will be thrown if it does not support the inputs. There are cases when fftconvolve supports the inputs but this function returns direct (e.g., to protect against floating point integer precision).

New in version 0.19.

Estimate the fastest method for a given input:

>>> from scipy import signal
>>> a = np.random.randn(1000)
>>> b = np.random.randn(1000000)
>>> method = signal.choose_conv_method(a, b, mode='same')
>>> method
'fft'

This can then be applied to other arrays of the same dtype and shape:

>>> c = np.random.randn(1000)
>>> d = np.random.randn(1000000)
>>> # `method` works with correlate and convolve
>>> corr1 = signal.correlate(a, b, mode='same', method=method)
>>> corr2 = signal.correlate(c, d, mode='same', method=method)
>>> conv1 = signal.convolve(a, b, mode='same', method=method)
>>> conv2 = signal.convolve(c, d, mode='same', method=method)
convolve(in1, in2, mode="full", method="auto")

Convolve two N-dimensional arrays.

Convolve in1 and in2, with the output size determined by the mode argument.

in1 : array_like
First input.
in2 : array_like
Second input. Should have the same number of dimensions as in1.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
method : str {‘auto’, ‘direct’, ‘fft’}, optional

A string indicating which method to use to calculate the convolution.

direct
The convolution is determined directly from sums, the definition of convolution.
fft
The Fourier Transform is used to perform the convolution by calling fftconvolve.
auto

Automatically chooses direct or Fourier method based on an estimate of which is faster (default). See Notes for more detail.

New in version 0.19.0.

convolve : array
An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.
numpy.polymul : performs polynomial multiplication (same operation, but
also accepts poly1d objects)

choose_conv_method : chooses the fastest appropriate convolution method fftconvolve

By default, convolve and correlate use method='auto', which calls choose_conv_method to choose the fastest method using pre-computed values (choose_conv_method can also measure real-world timing with a keyword argument). Because fftconvolve relies on floating point numbers, there are certain constraints that may force method=direct (more detail in choose_conv_method docstring).

Smooth a square pulse using a Hann window:

>>> from scipy import signal
>>> sig = np.repeat([0., 1., 0.], 100)
>>> win = signal.hann(50)
>>> filtered = signal.convolve(sig, win, mode='same') / sum(win)
>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True)
>>> ax_orig.plot(sig)
>>> ax_orig.set_title('Original pulse')
>>> ax_orig.margins(0, 0.1)
>>> ax_win.plot(win)
>>> ax_win.set_title('Filter impulse response')
>>> ax_win.margins(0, 0.1)
>>> ax_filt.plot(filtered)
>>> ax_filt.set_title('Filtered signal')
>>> ax_filt.margins(0, 0.1)
>>> fig.tight_layout()
>>> fig.show()
order_filter(a, domain, rank)

Perform an order filter on an N-dimensional array.

Perform an order filter on the array in. The domain argument acts as a mask centered over each pixel. The non-zero elements of domain are used to select elements surrounding each input pixel which are placed in a list. The list is sorted, and the output for that pixel is the element corresponding to rank in the sorted list.

a : ndarray
The N-dimensional input array.
domain : array_like
A mask array with the same number of dimensions as a. Each dimension should have an odd number of elements.
rank : int
A non-negative integer which selects the element from the sorted list (0 corresponds to the smallest element, 1 is the next smallest element, etc.).
out : ndarray
The results of the order filter in an array with the same shape as a.
>>> from scipy import signal
>>> x = np.arange(25).reshape(5, 5)
>>> domain = np.identity(3)
>>> x
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])
>>> signal.order_filter(x, domain, 0)
array([[  0.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   1.,   2.,   0.],
       [  0.,   5.,   6.,   7.,   0.],
       [  0.,  10.,  11.,  12.,   0.],
       [  0.,   0.,   0.,   0.,   0.]])
>>> signal.order_filter(x, domain, 2)
array([[  6.,   7.,   8.,   9.,   4.],
       [ 11.,  12.,  13.,  14.,   9.],
       [ 16.,  17.,  18.,  19.,  14.],
       [ 21.,  22.,  23.,  24.,  19.],
       [ 20.,  21.,  22.,  23.,  24.]])
medfilt(volume, kernel_size=None)

Perform a median filter on an N-dimensional array.

Apply a median filter to the input array using a local window-size given by kernel_size.

volume : array_like
An N-dimensional input array.
kernel_size : array_like, optional
A scalar or an N-length list giving the size of the median filter window in each dimension. Elements of kernel_size should be odd. If kernel_size is a scalar, then this scalar is used as the size in each dimension. Default size is 3 for each dimension.
out : ndarray
An array the same size as input containing the median filtered result.
wiener(im, mysize=None, noise=None)

Perform a Wiener filter on an N-dimensional array.

Apply a Wiener filter to the N-dimensional array im.

im : ndarray
An N-dimensional array.
mysize : int or array_like, optional
A scalar or an N-length list giving the size of the Wiener filter window in each dimension. Elements of mysize should be odd. If mysize is a scalar, then this scalar is used as the size in each dimension.
noise : float, optional
The noise-power to use. If None, then noise is estimated as the average of the local variance of the input.
out : ndarray
Wiener filtered result with the same shape as im.
convolve2d(in1, in2, mode="full", boundary="fill", fillvalue=0)

Convolve two 2-dimensional arrays.

Convolve in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue.

in1 : array_like
First input.
in2 : array_like
Second input. Should have the same number of dimensions as in1. If operating in ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
boundary : str {‘fill’, ‘wrap’, ‘symm’}, optional

A flag indicating how to handle boundaries:

fill
pad input arrays with fillvalue. (default)
wrap
circular boundary conditions.
symm
symmetrical boundary conditions.
fillvalue : scalar, optional
Value to fill pad input arrays with. Default is 0.
out : ndarray
A 2-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

Compute the gradient of an image by 2D convolution with a complex Scharr operator. (Horizontal operator is real, vertical is imaginary.) Use symmetric boundary condition to avoid creating edges at the image boundaries.

>>> from scipy import signal
>>> from scipy import misc
>>> ascent = misc.ascent()
>>> scharr = np.array([[ -3-3j, 0-10j,  +3 -3j],
...                    [-10+0j, 0+ 0j, +10 +0j],
...                    [ -3+3j, 0+10j,  +3 +3j]]) # Gx + j*Gy
>>> grad = signal.convolve2d(ascent, scharr, boundary='symm', mode='same')
>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_mag, ax_ang) = plt.subplots(3, 1, figsize=(6, 15))
>>> ax_orig.imshow(ascent, cmap='gray')
>>> ax_orig.set_title('Original')
>>> ax_orig.set_axis_off()
>>> ax_mag.imshow(np.absolute(grad), cmap='gray')
>>> ax_mag.set_title('Gradient magnitude')
>>> ax_mag.set_axis_off()
>>> ax_ang.imshow(np.angle(grad), cmap='hsv') # hsv is cyclic, like angles
>>> ax_ang.set_title('Gradient orientation')
>>> ax_ang.set_axis_off()
>>> fig.show()
correlate2d(in1, in2, mode="full", boundary="fill", fillvalue=0)

Cross-correlate two 2-dimensional arrays.

Cross correlate in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue.

in1 : array_like
First input.
in2 : array_like
Second input. Should have the same number of dimensions as in1. If operating in ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
mode : str {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full
The output is the full discrete linear cross-correlation of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
boundary : str {‘fill’, ‘wrap’, ‘symm’}, optional

A flag indicating how to handle boundaries:

fill
pad input arrays with fillvalue. (default)
wrap
circular boundary conditions.
symm
symmetrical boundary conditions.
fillvalue : scalar, optional
Value to fill pad input arrays with. Default is 0.
correlate2d : ndarray
A 2-dimensional array containing a subset of the discrete linear cross-correlation of in1 with in2.

Use 2D cross-correlation to find the location of a template in a noisy image:

>>> from scipy import signal
>>> from scipy import misc
>>> face = misc.face(gray=True) - misc.face(gray=True).mean()
>>> template = np.copy(face[300:365, 670:750])  # right eye
>>> template -= template.mean()
>>> face = face + np.random.randn(*face.shape) * 50  # add noise
>>> corr = signal.correlate2d(face, template, boundary='symm', mode='same')
>>> y, x = np.unravel_index(np.argmax(corr), corr.shape)  # find the match
>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_template, ax_corr) = plt.subplots(3, 1,
...                                                     figsize=(6, 15))
>>> ax_orig.imshow(face, cmap='gray')
>>> ax_orig.set_title('Original')
>>> ax_orig.set_axis_off()
>>> ax_template.imshow(template, cmap='gray')
>>> ax_template.set_title('Template')
>>> ax_template.set_axis_off()
>>> ax_corr.imshow(corr, cmap='gray')
>>> ax_corr.set_title('Cross-correlation')
>>> ax_corr.set_axis_off()
>>> ax_orig.plot(x, y, 'ro')
>>> fig.show()
medfilt2d(input, kernel_size=3)

Median filter a 2-dimensional array.

Apply a median filter to the input array using a local window-size given by kernel_size (must be odd).

input : array_like
A 2-dimensional input array.
kernel_size : array_like, optional
A scalar or a list of length 2, giving the size of the median filter window in each dimension. Elements of kernel_size should be odd. If kernel_size is a scalar, then this scalar is used as the size in each dimension. Default is a kernel of size (3, 3).
out : ndarray
An array the same size as input containing the median filtered result.
lfilter(b, a, x, axis=None, zi=None)

Filter data along one-dimension with an IIR or FIR filter.

Filter a data sequence, x, using a digital filter. This works for many fundamental data types (including Object type). The filter is a direct form II transposed implementation of the standard difference equation (see Notes).

b : array_like
The numerator coefficient vector in a 1-D sequence.
a : array_like
The denominator coefficient vector in a 1-D sequence. If a[0] is not 1, then both a and b are normalized by a[0].
x : array_like
An N-dimensional input array.
axis : int, optional
The axis of the input data array along which to apply the linear filter. The filter is applied to each subarray along this axis. Default is -1.
zi : array_like, optional
Initial conditions for the filter delays. It is a vector (or array of vectors for an N-dimensional input) of length max(len(a), len(b)) - 1. If zi is None or is not given then initial rest is assumed. See lfiltic for more information.
y : array
The output of the digital filter.
zf : array, optional
If zi is None, this is not returned, otherwise, zf holds the final filter delay values.

lfiltic : Construct initial conditions for lfilter. lfilter_zi : Compute initial state (steady state of step response) for

lfilter.

filtfilt : A forward-backward filter, to obtain a filter with linear phase. savgol_filter : A Savitzky-Golay filter. sosfilt: Filter data using cascaded second-order sections. sosfiltfilt: A forward-backward filter using second-order sections.

The filter function is implemented as a direct II transposed structure. This means that the filter implements:

a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M]
                      - a[1]*y[n-1] - ... - a[N]*y[n-N]

where M is the degree of the numerator, N is the degree of the denominator, and n is the sample number. It is implemented using the following difference equations (assuming M = N):

a[0]*y[n] = b[0] * x[n]               + d[0][n-1]
  d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n-1]
  d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n-1]
...
d[N-2][n] = b[N-1]*x[n] - a[N-1]*y[n] + d[N-1][n-1]
d[N-1][n] = b[N] * x[n] - a[N] * y[n]

where d are the state variables.

The rational transfer function describing this filter in the z-transform domain is:

                    -1              -M
        b[0] + b[1]z  + ... + b[M] z
Y(z) = -------------------------------- X(z)
                    -1              -N
        a[0] + a[1]z  + ... + a[N] z

Generate a noisy signal to be filtered:

>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t = np.linspace(-1, 1, 201)
>>> x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +
...      0.1*np.sin(2*np.pi*1.25*t + 1) +
...      0.18*np.cos(2*np.pi*3.85*t))
>>> xn = x + np.random.randn(len(t)) * 0.08

Create an order 3 lowpass butterworth filter:

>>> b, a = signal.butter(3, 0.05)

Apply the filter to xn. Use lfilter_zi to choose the initial condition of the filter:

>>> zi = signal.lfilter_zi(b, a)
>>> z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0])

Apply the filter again, to have a result filtered at an order the same as filtfilt:

>>> z2, _ = signal.lfilter(b, a, z, zi=zi*z[0])

Use filtfilt to apply the filter:

>>> y = signal.filtfilt(b, a, xn)

Plot the original signal and the various filtered versions:

>>> plt.figure
>>> plt.plot(t, xn, 'b', alpha=0.75)
>>> plt.plot(t, z, 'r--', t, z2, 'r', t, y, 'k')
>>> plt.legend(('noisy signal', 'lfilter, once', 'lfilter, twice',
...             'filtfilt'), loc='best')
>>> plt.grid(True)
>>> plt.show()
lfiltic(b, a, y, x=None)

Construct initial conditions for lfilter given input and output vectors.

Given a linear filter (b, a) and initial conditions on the output y and the input x, return the initial conditions on the state vector zi which is used by lfilter to generate the output given the input.

b : array_like
Linear filter term.
a : array_like
Linear filter term.
y : array_like

Initial conditions.

If N = len(a) - 1, then y = {y[-1], y[-2], ..., y[-N]}.

If y is too short, it is padded with zeros.

x : array_like, optional

Initial conditions.

If M = len(b) - 1, then x = {x[-1], x[-2], ..., x[-M]}.

If x is not given, its initial conditions are assumed zero.

If x is too short, it is padded with zeros.

zi : ndarray
The state vector zi = {z_0[-1], z_1[-1], ..., z_K-1[-1]}, where K = max(M, N).

lfilter, lfilter_zi

deconvolve(signal, divisor)

Deconvolves divisor out of signal using inverse filtering.

Returns the quotient and remainder such that signal = convolve(divisor, quotient) + remainder

signal : array_like
Signal data, typically a recorded signal
divisor : array_like
Divisor data, typically an impulse response or filter that was applied to the original signal
quotient : ndarray
Quotient, typically the recovered original signal
remainder : ndarray
Remainder

Deconvolve a signal that’s been filtered:

>>> from scipy import signal
>>> original = [0, 1, 0, 0, 1, 1, 0, 0]
>>> impulse_response = [2, 1]
>>> recorded = signal.convolve(impulse_response, original)
>>> recorded
array([0, 2, 1, 0, 2, 3, 1, 0, 0])
>>> recovered, remainder = signal.deconvolve(recorded, impulse_response)
>>> recovered
array([ 0.,  1.,  0.,  0.,  1.,  1.,  0.,  0.])
numpy.polydiv : performs polynomial division (same operation, but
also accepts poly1d objects)
hilbert(x, N=None, axis=None)

Compute the analytic signal, using the Hilbert transform.

The transformation is done along the last axis by default.

x : array_like
Signal data. Must be real.
N : int, optional
Number of Fourier components. Default: x.shape[axis]
axis : int, optional
Axis along which to do the transformation. Default: -1.
xa : ndarray
Analytic signal of x, of each 1-D array along axis

scipy.fftpack.hilbert : Return Hilbert transform of a periodic sequence x.

The analytic signal x_a(t) of signal x(t) is:

where F is the Fourier transform, U the unit step function, and y the Hilbert transform of x. [1]_

In other words, the negative half of the frequency spectrum is zeroed out, turning the real-valued signal into a complex signal. The Hilbert transformed signal can be obtained from np.imag(hilbert(x)), and the original signal from np.real(hilbert(x)).

In this example we use the Hilbert transform to determine the amplitude envelope and instantaneous frequency of an amplitude-modulated signal.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.signal import hilbert, chirp
>>> duration = 1.0
>>> fs = 400.0
>>> samples = int(fs*duration)
>>> t = np.arange(samples) / fs

We create a chirp of which the frequency increases from 20 Hz to 100 Hz and apply an amplitude modulation.

>>> signal = chirp(t, 20.0, t[-1], 100.0)
>>> signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) )

The amplitude envelope is given by magnitude of the analytic signal. The instantaneous frequency can be obtained by differentiating the instantaneous phase in respect to time. The instantaneous phase corresponds to the phase angle of the analytic signal.

>>> analytic_signal = hilbert(signal)
>>> amplitude_envelope = np.abs(analytic_signal)
>>> instantaneous_phase = np.unwrap(np.angle(analytic_signal))
>>> instantaneous_frequency = (np.diff(instantaneous_phase) /
...                            (2.0*np.pi) * fs)
>>> fig = plt.figure()
>>> ax0 = fig.add_subplot(211)
>>> ax0.plot(t, signal, label='signal')
>>> ax0.plot(t, amplitude_envelope, label='envelope')
>>> ax0.set_xlabel("time in seconds")
>>> ax0.legend()
>>> ax1 = fig.add_subplot(212)
>>> ax1.plot(t[1:], instantaneous_frequency)
>>> ax1.set_xlabel("time in seconds")
>>> ax1.set_ylim(0.0, 120.0)
[1]Wikipedia, “Analytic signal”. http://en.wikipedia.org/wiki/Analytic_signal
[2]Leon Cohen, “Time-Frequency Analysis”, 1995. Chapter 2.
[3]Alan V. Oppenheim, Ronald W. Schafer. Discrete-Time Signal Processing, Third Edition, 2009. Chapter 12. ISBN 13: 978-1292-02572-8
hilbert2(x, N=None)

Compute the ‘2-D’ analytic signal of x

x : array_like
2-D signal data.
N : int or tuple of two ints, optional
Number of Fourier components. Default is x.shape
xa : ndarray
Analytic signal of x taken along axes (0,1).
[1]Wikipedia, “Analytic signal”, http://en.wikipedia.org/wiki/Analytic_signal
cmplx_sort(p)

Sort roots based on magnitude.

p : array_like
The roots to sort, as a 1-D array.
p_sorted : ndarray
Sorted roots.
indx : ndarray
Array of indices needed to sort the input p.
unique_roots(p, tol=0.001, rtype="min")

Determine unique roots and their multiplicities from a list of roots.

p : array_like
The list of roots.
tol : float, optional
The tolerance for two roots to be considered equal. Default is 1e-3.
rtype : {‘max’, ‘min, ‘avg’}, optional

How to determine the returned root if multiple roots are within tol of each other.

  • ‘max’: pick the maximum of those roots.
  • ‘min’: pick the minimum of those roots.
  • ‘avg’: take the average of those roots.
pout : ndarray
The list of unique roots, sorted from low to high.
mult : ndarray
The multiplicity of each root.

This utility function is not specific to roots but can be used for any sequence of values for which uniqueness and multiplicity has to be determined. For a more general routine, see numpy.unique.

>>> from scipy import signal
>>> vals = [0, 1.3, 1.31, 2.8, 1.25, 2.2, 10.3]
>>> uniq, mult = signal.unique_roots(vals, tol=2e-2, rtype='avg')

Check which roots have multiplicity larger than 1:

>>> uniq[mult > 1]
array([ 1.305])
invres(r, p, k, tol=0.001, rtype="avg")

Compute b(s) and a(s) from partial fraction expansion.

If M is the degree of numerator b and N the degree of denominator a:

        b(s)     b[0] s**(M) + b[1] s**(M-1) + ... + b[M]
H(s) = ------ = ------------------------------------------
        a(s)     a[0] s**(N) + a[1] s**(N-1) + ... + a[N]

then the partial-fraction expansion H(s) is defined as:

    r[0]       r[1]             r[-1]
= -------- + -------- + ... + --------- + k(s)
  (s-p[0])   (s-p[1])         (s-p[-1])

If there are any repeated roots (closer together than tol), then H(s) has terms like:

  r[i]      r[i+1]              r[i+n-1]
-------- + ----------- + ... + -----------
(s-p[i])  (s-p[i])**2          (s-p[i])**n

This function is used for polynomials in positive powers of s or z, such as analog filters or digital filters in controls engineering. For negative powers of z (typical for digital filters in DSP), use invresz.

r : array_like
Residues.
p : array_like
Poles.
k : array_like
Coefficients of the direct polynomial term.
tol : float, optional
The tolerance for two roots to be considered equal. Default is 1e-3.
rtype : {‘max’, ‘min, ‘avg’}, optional

How to determine the returned root if multiple roots are within tol of each other.

  • ‘max’: pick the maximum of those roots.
  • ‘min’: pick the minimum of those roots.
  • ‘avg’: take the average of those roots.
b : ndarray
Numerator polynomial coefficients.
a : ndarray
Denominator polynomial coefficients.

residue, invresz, unique_roots

residue(b, a, tol=0.001, rtype="avg")

Compute partial-fraction expansion of b(s) / a(s).

If M is the degree of numerator b and N the degree of denominator a:

        b(s)     b[0] s**(M) + b[1] s**(M-1) + ... + b[M]
H(s) = ------ = ------------------------------------------
        a(s)     a[0] s**(N) + a[1] s**(N-1) + ... + a[N]

then the partial-fraction expansion H(s) is defined as:

    r[0]       r[1]             r[-1]
= -------- + -------- + ... + --------- + k(s)
  (s-p[0])   (s-p[1])         (s-p[-1])

If there are any repeated roots (closer together than tol), then H(s) has terms like:

  r[i]      r[i+1]              r[i+n-1]
-------- + ----------- + ... + -----------
(s-p[i])  (s-p[i])**2          (s-p[i])**n

This function is used for polynomials in positive powers of s or z, such as analog filters or digital filters in controls engineering. For negative powers of z (typical for digital filters in DSP), use residuez.

b : array_like
Numerator polynomial coefficients.
a : array_like
Denominator polynomial coefficients.
r : ndarray
Residues.
p : ndarray
Poles.
k : ndarray
Coefficients of the direct polynomial term.

invres, residuez, numpy.poly, unique_roots

residuez(b, a, tol=0.001, rtype="avg")

Compute partial-fraction expansion of b(z) / a(z).

If M is the degree of numerator b and N the degree of denominator a:

        b(z)     b[0] + b[1] z**(-1) + ... + b[M] z**(-M)
H(z) = ------ = ------------------------------------------
        a(z)     a[0] + a[1] z**(-1) + ... + a[N] z**(-N)

then the partial-fraction expansion H(z) is defined as:

        r[0]                   r[-1]
= --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ...
  (1-p[0]z**(-1))         (1-p[-1]z**(-1))

If there are any repeated roots (closer than tol), then the partial fraction expansion has terms like:

     r[i]              r[i+1]                    r[i+n-1]
-------------- + ------------------ + ... + ------------------
(1-p[i]z**(-1))  (1-p[i]z**(-1))**2         (1-p[i]z**(-1))**n

This function is used for polynomials in negative powers of z, such as digital filters in DSP. For positive powers, use residue.

b : array_like
Numerator polynomial coefficients.
a : array_like
Denominator polynomial coefficients.
r : ndarray
Residues.
p : ndarray
Poles.
k : ndarray
Coefficients of the direct polynomial term.

invresz, residue, unique_roots

invresz(r, p, k, tol=0.001, rtype="avg")

Compute b(z) and a(z) from partial fraction expansion.

If M is the degree of numerator b and N the degree of denominator a:

        b(z)     b[0] + b[1] z**(-1) + ... + b[M] z**(-M)
H(z) = ------ = ------------------------------------------
        a(z)     a[0] + a[1] z**(-1) + ... + a[N] z**(-N)

then the partial-fraction expansion H(z) is defined as:

        r[0]                   r[-1]
= --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ...
  (1-p[0]z**(-1))         (1-p[-1]z**(-1))

If there are any repeated roots (closer than tol), then the partial fraction expansion has terms like:

     r[i]              r[i+1]                    r[i+n-1]
-------------- + ------------------ + ... + ------------------
(1-p[i]z**(-1))  (1-p[i]z**(-1))**2         (1-p[i]z**(-1))**n

This function is used for polynomials in negative powers of z, such as digital filters in DSP. For positive powers, use invres.

r : array_like
Residues.
p : array_like
Poles.
k : array_like
Coefficients of the direct polynomial term.
tol : float, optional
The tolerance for two roots to be considered equal. Default is 1e-3.
rtype : {‘max’, ‘min, ‘avg’}, optional

How to determine the returned root if multiple roots are within tol of each other.

  • ‘max’: pick the maximum of those roots.
  • ‘min’: pick the minimum of those roots.
  • ‘avg’: take the average of those roots.
b : ndarray
Numerator polynomial coefficients.
a : ndarray
Denominator polynomial coefficients.

residuez, unique_roots, invres

resample(x, num, t=None, axis=0, window=None)

Resample x to num samples using Fourier method along the given axis.

The resampled signal starts at the same value as x but is sampled with a spacing of len(x) / num * (spacing of x). Because a Fourier method is used, the signal is assumed to be periodic.

x : array_like
The data to be resampled.
num : int
The number of samples in the resampled signal.
t : array_like, optional
If t is given, it is assumed to be the sample positions associated with the signal data in x.
axis : int, optional
The axis of x that is resampled. Default is 0.
window : array_like, callable, string, float, or tuple, optional
Specifies the window applied to the signal in the Fourier domain. See below for details.
resampled_x or (resampled_x, resampled_t)
Either the resampled array, or, if t was given, a tuple containing the resampled array and the corresponding resampled positions.

decimate : Downsample the signal after applying an FIR or IIR filter. resample_poly : Resample using polyphase filtering and an FIR filter.

The argument window controls a Fourier-domain window that tapers the Fourier spectrum before zero-padding to alleviate ringing in the resampled values for sampled signals you didn’t intend to be interpreted as band-limited.

If window is a function, then it is called with a vector of inputs indicating the frequency bins (i.e. fftfreq(x.shape[axis]) ).

If window is an array of the same length as x.shape[axis] it is assumed to be the window to be applied directly in the Fourier domain (with dc and low-frequency first).

For any other type of window, the function scipy.signal.get_window is called to generate the window.

The first sample of the returned vector is the same as the first sample of the input vector. The spacing between samples is changed from dx to dx * len(x) / num.

If t is not None, then it represents the old sample positions, and the new sample positions will be returned as well as the new samples.

As noted, resample uses FFT transformations, which can be very slow if the number of input or output samples is large and prime; see scipy.fftpack.fft.

Note that the end of the resampled data rises to meet the first sample of the next cycle:

>>> from scipy import signal
>>> x = np.linspace(0, 10, 20, endpoint=False)
>>> y = np.cos(-x**2/6.0)
>>> f = signal.resample(y, 100)
>>> xnew = np.linspace(0, 10, 100, endpoint=False)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'go-', xnew, f, '.-', 10, y[0], 'ro')
>>> plt.legend(['data', 'resampled'], loc='best')
>>> plt.show()
resample_poly(x, up, down, axis=0, window=tuple)

Resample x along the given axis using polyphase filtering.

The signal x is upsampled by the factor up, a zero-phase low-pass FIR filter is applied, and then it is downsampled by the factor down. The resulting sample rate is up / down times the original sample rate. Values beyond the boundary of the signal are assumed to be zero during the filtering step.

x : array_like
The data to be resampled.
up : int
The upsampling factor.
down : int
The downsampling factor.
axis : int, optional
The axis of x that is resampled. Default is 0.
window : string, tuple, or array_like, optional
Desired window to use to design the low-pass filter, or the FIR filter coefficients to employ. See below for details.
resampled_x : array
The resampled array.

decimate : Downsample the signal after applying an FIR or IIR filter. resample : Resample up or down using the FFT method.

This polyphase method will likely be faster than the Fourier method in scipy.signal.resample when the number of samples is large and prime, or when the number of samples is large and up and down share a large greatest common denominator. The length of the FIR filter used will depend on max(up, down) // gcd(up, down), and the number of operations during polyphase filtering will depend on the filter length and down (see scipy.signal.upfirdn for details).

The argument window specifies the FIR low-pass filter design.

If window is an array_like it is assumed to be the FIR filter coefficients. Note that the FIR filter is applied after the upsampling step, so it should be designed to operate on a signal at a sampling frequency higher than the original by a factor of up//gcd(up, down). This function’s output will be centered with respect to this array, so it is best to pass a symmetric filter with an odd number of samples if, as is usually the case, a zero-phase filter is desired.

For any other type of window, the functions scipy.signal.get_window and scipy.signal.firwin are called to generate the appropriate filter coefficients.

The first sample of the returned vector is the same as the first sample of the input vector. The spacing between samples is changed from dx to dx * down / float(up).

Note that the end of the resampled data rises to meet the first sample of the next cycle for the FFT method, and gets closer to zero for the polyphase method:

>>> from scipy import signal
>>> x = np.linspace(0, 10, 20, endpoint=False)
>>> y = np.cos(-x**2/6.0)
>>> f_fft = signal.resample(y, 100)
>>> f_poly = signal.resample_poly(y, 100, 20)
>>> xnew = np.linspace(0, 10, 100, endpoint=False)
>>> import matplotlib.pyplot as plt
>>> plt.plot(xnew, f_fft, 'b.-', xnew, f_poly, 'r.-')
>>> plt.plot(x, y, 'ko-')
>>> plt.plot(10, y[0], 'bo', 10, 0., 'ro')  # boundaries
>>> plt.legend(['resample', 'resamp_poly', 'data'], loc='best')
>>> plt.show()
vectorstrength(events, period)

Determine the vector strength of the events corresponding to the given period.

The vector strength is a measure of phase synchrony, how well the timing of the events is synchronized to a single period of a periodic signal.

If multiple periods are used, calculate the vector strength of each. This is called the “resonating vector strength”.

events : 1D array_like
An array of time points containing the timing of the events.
period : float or array_like
The period of the signal that the events should synchronize to. The period is in the same units as events. It can also be an array of periods, in which case the outputs are arrays of the same length.
strength : float or 1D array
The strength of the synchronization. 1.0 is perfect synchronization and 0.0 is no synchronization. If period is an array, this is also an array with each element containing the vector strength at the corresponding period.
phase : float or array
The phase that the events are most strongly synchronized to in radians. If period is an array, this is also an array with each element containing the phase for the corresponding period.
van Hemmen, JL, Longtin, A, and Vollmayr, AN. Testing resonating vector
strength: Auditory system, electric fish, and noise. Chaos 21, 047508 (2011); :doi:`10.1063/1.3670512`.
van Hemmen, JL. Vector strength after Goldberg, Brown, and von Mises:
biological and mathematical perspectives. Biol Cybern. 2013 Aug;107(4):385-96. :doi:`10.1007/s00422-013-0561-7`.
van Hemmen, JL and Vollmayr, AN. Resonating vector strength: what happens
when we vary the “probing” frequency while keeping the spike times fixed. Biol Cybern. 2013 Aug;107(4):491-94. :doi:`10.1007/s00422-013-0560-8`.
detrend(data, axis=None, type="linear", bp=0)

Remove linear trend along axis from data.

data : array_like
The input data.
axis : int, optional
The axis along which to detrend the data. By default this is the last axis (-1).
type : {‘linear’, ‘constant’}, optional
The type of detrending. If type == 'linear' (default), the result of a linear least-squares fit to data is subtracted from data. If type == 'constant', only the mean of data is subtracted.
bp : array_like of ints, optional
A sequence of break points. If given, an individual linear fit is performed for each part of data between two break points. Break points are specified as indices into data.
ret : ndarray
The detrended input data.
>>> from scipy import signal
>>> randgen = np.random.RandomState(9)
>>> npoints = 1000
>>> noise = randgen.randn(npoints)
>>> x = 3 + 2*np.linspace(0, 1, npoints) + noise
>>> (signal.detrend(x) - noise).max() < 0.01
True
lfilter_zi(b, a)

Construct initial conditions for lfilter for step response steady-state.

Compute an initial state zi for the lfilter function that corresponds to the steady state of the step response.

A typical use of this function is to set the initial state so that the output of the filter starts at the same value as the first element of the signal to be filtered.

b, a : array_like (1-D)
The IIR filter coefficients. See lfilter for more information.
zi : 1-D ndarray
The initial state for the filter.

lfilter, lfiltic, filtfilt

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) = A*z(n) + B*x(n)
y(n)   = C*z(n) + D*x(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.

Given the filter coefficients a and b, the state space matrices for the transposed direct form II implementation of the linear filter, which is the implementation used by scipy.signal.lfilter, are:

A = scipy.linalg.companion(a).T
B = b[1:] - a[1:]*b[0]

assuming a[0] is 1.0; if a[0] is not 1, a and b are first divided by a[0].

The following code creates a lowpass Butterworth filter. Then it applies that filter to an array whose values are all 1.0; the output is also all 1.0, as expected for a lowpass filter. If the zi argument of lfilter had not been given, the output would have shown the transient signal.

>>> from numpy import array, ones
>>> from scipy.signal import lfilter, lfilter_zi, butter
>>> b, a = butter(5, 0.25)
>>> zi = lfilter_zi(b, a)
>>> y, zo = lfilter(b, a, ones(10), zi=zi)
>>> y
array([1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])

Another example:

>>> x = array([0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0])
>>> y, zf = lfilter(b, a, x, zi=zi*x[0])
>>> y
array([ 0.5       ,  0.5       ,  0.5       ,  0.49836039,  0.48610528,
    0.44399389,  0.35505241])

Note that the zi argument to lfilter was computed using lfilter_zi and scaled by x[0]. Then the output y has no transient until the input drops from 0.5 to 0.0.

sosfilt_zi(sos)

Construct initial conditions for sosfilt for step response steady-state.

Compute an initial state zi for the sosfilt function that corresponds to the steady state of the step response.

A typical use of this function is to set the initial state so that the output of the filter starts at the same value as the first element of the signal to be filtered.

sos : array_like
Array of second-order filter coefficients, must have shape (n_sections, 6). See sosfilt for the SOS filter format specification.
zi : ndarray
Initial conditions suitable for use with sosfilt, shape (n_sections, 2).

sosfilt, zpk2sos

New in version 0.16.0.

Filter a rectangular pulse that begins at time 0, with and without the use of the zi argument of scipy.signal.sosfilt.

>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> sos = signal.butter(9, 0.125, output='sos')
>>> zi = signal.sosfilt_zi(sos)
>>> x = (np.arange(250) < 100).astype(int)
>>> f1 = signal.sosfilt(sos, x)
>>> f2, zo = signal.sosfilt(sos, x, zi=zi)
>>> plt.plot(x, 'k--', label='x')
>>> plt.plot(f1, 'b', alpha=0.5, linewidth=2, label='filtered')
>>> plt.plot(f2, 'g', alpha=0.25, linewidth=4, label='filtered with zi')
>>> plt.legend(loc='best')
>>> plt.show()
_filtfilt_gust(b, a, x, axis=None, irlen=None)

Forward-backward IIR filter that uses Gustafsson’s method.

Apply the IIR filter defined by (b,a) to x twice, first forward then backward, using Gustafsson’s initial conditions [1]_.

Let y_fb be the result of filtering first forward and then backward, and let y_bf be the result of filtering first backward then forward. Gustafsson’s method is to compute initial conditions for the forward pass and the backward pass such that y_fb == y_bf.

b : scalar or 1-D ndarray
Numerator coefficients of the filter.
a : scalar or 1-D ndarray
Denominator coefficients of the filter.
x : ndarray
Data to be filtered.
axis : int, optional
Axis of x to be filtered. Default is -1.
irlen : int or None, optional
The length of the nonnegligible part of the impulse response. If irlen is None, or if the length of the signal is less than 2 * irlen, then no part of the impulse response is ignored.
y : ndarray
The filtered data.
x0 : ndarray
Initial condition for the forward filter.
x1 : ndarray
Initial condition for the backward filter.

Typically the return values x0 and x1 are not needed by the caller. The intended use of these return values is in unit tests.

[1]F. Gustaffson. Determining the initial states in forward-backward filtering. Transactions on Signal Processing, 46(4):988-992, 1996.
filtfilt(b, a, x, axis=None, padtype="odd", padlen=None, method="pad", irlen=None)

Apply a digital filter forward and backward to a signal.

This function applies a linear digital filter twice, once forward and once backwards. The combined filter has zero phase and a filter order twice that of the original.

The function provides options for handling the edges of the signal.

b : (N,) array_like
The numerator coefficient vector of the filter.
a : (N,) array_like
The denominator coefficient vector of the filter. If a[0] is not 1, then both a and b are normalized by a[0].
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)).
method : str, optional
Determines the method for handling the edges of the signal, either “pad” or “gust”. When method is “pad”, the signal is padded; the type of padding is determined by padtype and padlen, and irlen is ignored. When method is “gust”, Gustafsson’s method is used, and padtype and padlen are ignored.
irlen : int or None, optional
When method is “gust”, irlen specifies the length of the impulse response of the filter. If irlen is None, no part of the impulse response is ignored. For a long signal, specifying irlen can significantly improve the performance of the filter.
y : ndarray
The filtered output with the same shape as x.

sosfiltfilt, lfilter_zi, lfilter, lfiltic, savgol_filter, sosfilt

When method is “pad”, the function pads 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 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.

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

The examples will use several functions from scipy.signal.

>>> 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 rate, 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.
>>> np.random.seed(123456)

sig is a random input signal to be filtered.

>>> n = 60
>>> sig = np.random.randn(n)**3 + 3*np.random.randn(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 = np.random.randn(5000)
>>> 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)))
1.80056858312e-10
_validate_pad(padtype, padlen, x, axis, ntaps)

Helper to validate padding for filtfilt

sosfilt(sos, x, axis=None, zi=None)

Filter data along one dimension using cascaded second-order sections.

Filter a data sequence, x, using a digital IIR filter defined by sos. This is implemented by performing lfilter for each second-order section. See lfilter for details.

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
An N-dimensional input array.
axis : int, optional
The axis of the input data array along which to apply the linear filter. The filter is applied to each subarray along this axis. Default is -1.
zi : array_like, optional
Initial conditions for the cascaded filter delays. It is a (at least 2D) vector of shape (n_sections, ..., 2, ...), where ..., 2, ... denotes the shape of x, but with x.shape[axis] replaced by 2. If zi is None or is not given then initial rest (i.e. all zeros) is assumed. Note that these initial conditions are not the same as the initial conditions given by lfiltic or lfilter_zi.
y : ndarray
The output of the digital filter.
zf : ndarray, optional
If zi is None, this is not returned, otherwise, zf holds the final filter delay values.

zpk2sos, sos2zpk, sosfilt_zi, sosfiltfilt, sosfreqz

The filter function is implemented as a series of second-order filters with direct-form II transposed structure. It is designed to minimize numerical precision errors for high-order filters.

New in version 0.16.0.

Plot a 13th-order filter’s impulse response using both lfilter and sosfilt, showing the instability that results from trying to do a 13th-order filter in a single stage (the numerical error pushes some poles outside of the unit circle):

>>> import matplotlib.pyplot as plt
>>> from scipy import signal
>>> b, a = signal.ellip(13, 0.009, 80, 0.05, output='ba')
>>> sos = signal.ellip(13, 0.009, 80, 0.05, output='sos')
>>> x = signal.unit_impulse(700)
>>> y_tf = signal.lfilter(b, a, x)
>>> y_sos = signal.sosfilt(sos, x)
>>> plt.plot(y_tf, 'r', label='TF')
>>> plt.plot(y_sos, 'k', label='SOS')
>>> plt.legend(loc='best')
>>> plt.show()
sosfiltfilt(sos, x, axis=None, padtype="odd", padlen=None)

A forward-backward digital filter using cascaded second-order sections.

See filtfilt for more complete information about this method.

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 * (2 * len(sos) + 1 - min((sos[:, 2] == 0).sum(),
                            (sos[:, 5] == 0).sum()))

The extra subtraction at the end attempts to compensate for poles and zeros at the origin (e.g. for odd-order filters) to yield equivalent estimates of padlen to those of filtfilt for second-order section filters built with scipy.signal functions.

y : ndarray
The filtered output with the same shape as x.

filtfilt, sosfilt, sosfilt_zi, sosfreqz

New in version 0.18.0.

>>> from scipy.signal import sosfiltfilt, butter
>>> import matplotlib.pyplot as plt

Create an interesting signal to filter.

>>> n = 201
>>> t = np.linspace(0, 1, n)
>>> np.random.seed(123)
>>> x = 1 + (t < 0.5) - 0.25*t**2 + 0.05*np.random.randn(n)

Create a lowpass Butterworth filter, and use it to filter x.

>>> sos = butter(4, 0.125, output='sos')
>>> y = sosfiltfilt(sos, x)

For comparison, apply an 8th order filter using sosfilt. The filter is initialized using the mean of the first four values of x.

>>> from scipy.signal import sosfilt, sosfilt_zi
>>> sos8 = butter(8, 0.125, output='sos')
>>> zi = x[:4].mean() * sosfilt_zi(sos8)
>>> y2, zo = sosfilt(sos8, x, zi=zi)

Plot the results. Note that the phase of y matches the input, while y2 has a significant phase delay.

>>> plt.plot(t, x, alpha=0.5, label='x(t)')
>>> plt.plot(t, y, label='y(t)')
>>> plt.plot(t, y2, label='y2(t)')
>>> plt.legend(framealpha=1, shadow=True)
>>> plt.grid(alpha=0.25)
>>> plt.xlabel('t')
>>> plt.show()
decimate(x, q, n=None, ftype="iir", axis=None, zero_phase=True)

Downsample the signal after applying an anti-aliasing filter.

By default, an order 8 Chebyshev type I filter is used. A 30 point FIR filter with Hamming window is used if ftype is ‘fir’.

x : array_like
The signal to be downsampled, as an N-dimensional array.
q : int
The downsampling factor. When using IIR downsampling, it is recommended to call decimate multiple times for downsampling factors higher than 13.
n : int, optional
The order of the filter (1 less than the length for ‘fir’). Defaults to 8 for ‘iir’ and 20 times the downsampling factor for ‘fir’.
ftype : str {‘iir’, ‘fir’} or dlti instance, optional
If ‘iir’ or ‘fir’, specifies the type of lowpass filter. If an instance of an dlti object, uses that object to filter before downsampling.
axis : int, optional
The axis along which to decimate.
zero_phase : bool, optional

Prevent phase shift by filtering with filtfilt instead of lfilter when using an IIR filter, and shifting the outputs back by the filter’s group delay when using an FIR filter. The default value of True is recommended, since a phase shift is generally not desired.

New in version 0.18.0.

y : ndarray
The down-sampled signal.

resample : Resample up or down using the FFT method. resample_poly : Resample using polyphase filtering and an FIR filter.

The zero_phase keyword was added in 0.18.0. The possibility to use instances of dlti as ftype was added in 0.18.0.