cuSignal API Reference#

Convolution#

Convolve#

cusignal.convolution.convolve.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.

Parameters:

in1array_like

The first argument passed into the convolution function.

in2array_like

The second argument passed into the convolution function.

modestr {‘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.

measurebool, 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.

Returns:

methodstr: A string indicating which convolution method is fastest, either ‘direct’ or ‘fft’
timesdict, optional: A dictionary containing the times (in seconds) needed for each method. This value is only returned if measure=True.

See also

convolve
correlate

Examples

Estimate the fastest method for a given input:

>>> import cusignal
>>> import cupy as cp
>>> a = cp.random.randn(1000)
>>> b = cp.random.randn(1000000)
>>> method = cusignal.choose_conv_method(a, b, mode='same')
>>> method
'fft'

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

>>> c = cp.random.randn(1000)
>>> d = cp.random.randn(1000000)
>>> # `method` works with correlate and convolve
>>> corr1 = cusignal.correlate(a, b, mode='same', method=method)
>>> corr2 = cusignal.correlate(c, d, mode='same', method=method)
>>> conv1 = cusignal.convolve(a, b, mode='same', method=method)
>>> conv2 = cusignal.convolve(c, d, mode='same', method=method)

cusignal.convolution.convolve.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.

Parameters:

in1array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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.

methodstr {‘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).

Returns:

convolvearray: An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

See also

choose_conv_method: chooses the fastest appropriate convolution method
fftconvolve

Notes

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).

Examples

Smooth a square pulse using a Hann window:

>>> import cusignal
>>> import cupy as cp
>>> sig = cp.repeat(cp.asarray([0., 1., 0.]), 100)
>>> win = cusignal.hann(50)
>>> filtered = cusignal.convolve(sig, win, mode='same') / cp.sum(win)

>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True)
>>> ax_orig.plot(cp.asnumpy(sig))
>>> ax_orig.set_title('Original pulse')
>>> ax_orig.margins(0, 0.1)
>>> ax_win.plot(cp.asnumpy(win))
>>> ax_win.set_title('Filter impulse response')
>>> ax_win.margins(0, 0.1)
>>> ax_filt.plot(cp.asnumpy(filtered))
>>> ax_filt.set_title('Filtered signal')
>>> ax_filt.margins(0, 0.1)
>>> fig.tight_layout()
>>> fig.show()

cusignal.convolution.convolve.convolve1d2o(in1, in2, mode='valid', method='direct')#

Convolve a 1-dimensional arrays with a 2nd order filter. This results in a second order convolution.

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

Parameters:

in1array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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.

methodstr {‘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).

Returns:

outndarray: A 1-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

See also

convolve
convolve1d2o
convolve1d3o

Examples

Convolution of a 2nd order filter on a 1d signal

>>> import cusignal as cs
>>> import numpy as np
>>> d = 50
>>> a = np.random.uniform(-1,1,(200))
>>> b = np.random.uniform(-1,1,(d,d))
>>> c = cs.convolve1d2o(a,b)

cusignal.convolution.convolve.convolve1d3o(in1, in2, mode='valid', method='direct')#

Convolve a 1-dimensional array with a 3rd order filter. This results in a second order convolution.

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

Parameters:

in1array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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.

methodstr {‘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).

Returns:

outndarray: A 1-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

See also

convolve
convolve1d2o
convolve1d3o

Examples

Convolution of a 3rd order filter on a 1d signal

>>> import cusignal as cs
>>> import numpy as np
>>> d = 50
>>> a = np.random.uniform(-1,1,(200))
>>> b = np.random.uniform(-1,1,(d,d,d))
>>> c = cs.convolve1d3o(a,b)

cusignal.convolution.convolve.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. Parameters ———- in1 : array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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.

boundarystr {‘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.

fillvaluescalar, optional

Value to fill pad input arrays with. Default is 0.

Returns:

outndarray: A 2-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

Examples

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.

>>> import cusignal
>>> import cupy as cp
>>> from scipy import misc
>>> ascent = cp.asarray(misc.ascent())
>>> scharr = cp.array([[ -3-3j, 0-10j,  +3 -3j],
...                    [-10+0j, 0+ 0j, +10 +0j],
...                    [ -3+3j, 0+10j,  +3 +3j]]) # Gx + j*Gy
>>> grad = cusignal.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(cp.asnumpy(ascent), cmap='gray')
>>> ax_orig.set_title('Original')
>>> ax_orig.set_axis_off()
>>> ax_mag.imshow(cp.asnumpy(cp.absolute(grad)), cmap='gray')
>>> ax_mag.set_title('Gradient magnitude')
>>> ax_mag.set_axis_off()
>>> ax_ang.imshow(cp.asnumpy(cp.angle(grad)), cmap='hsv')
>>> ax_ang.set_title('Gradient orientation')
>>> ax_ang.set_axis_off()
>>> fig.show()

cusignal.convolution.convolve.fftconvolve(in1, in2, mode='full', axes=None)#

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.

Parameters:

in1array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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. axis : tuple, optional

axesint or array_like of ints or None, optional

Axes over which to compute the convolution. The default is over all axes.

Returns:

outarray: An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.

Examples

Autocorrelation of white noise is an impulse.

>>> import cusignal
>>> import cupy as cp
>>> import numpy as np
>>> sig = cp.random.randn(1000)
>>> autocorr = cusignal.fftconvolve(sig, sig[::-1], mode='full')

>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1)
>>> ax_orig.plot(cp.asnumpy(sig))
>>> ax_orig.set_title('White noise')
>>> ax_mag.plot(cp.asnumpy(cp.arange(-len(sig)+1,len(sig))),
        cp.asnumpy(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 = cp.outer(cusignal.gaussian(70, 8), cusignal.gaussian(70, 8))
>>> blurred = cusignal.convolve(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(cp.asnumpy(kernel), cmap='gray')
>>> ax_kernel.set_title('Gaussian kernel')
>>> ax_kernel.set_axis_off()
>>> ax_blurred.imshow(cp.asnumpy(blurred), cmap='gray')
>>> ax_blurred.set_title('Blurred')
>>> ax_blurred.set_axis_off()
>>> fig.show()

Correlate#

cusignal.convolution.correlate.correlate(in1, in2, mode='full', method='auto')#

Cross-correlate two N-dimensional arrays.

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

Parameters:

in1array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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.

methodstr {‘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.

Returns:

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

See also

choose_conv_method: contains more documentation on method.

Notes

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

\[z[k] = (x * y)(k - N + 1) = \sum_{l=0}^{||x||-1}x_l y_{l-k+N-1}^{*}\]

for \(k = 0, 1, ..., ||x|| + ||y|| - 2\)

where \(||x||\) is the length of x, \(N = \max(||x||,||y||)\), and \(y_m\) 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.

Examples

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

>>> import cusignal
>>> import cupy as cp
>>> sig = cp.repeat(cp.array([0., 1., 1., 0., 1., 0., 0., 1.]), 128)
>>> sig_noise = sig + cp.random.randn(len(sig))
>>> corr = cusignal.correlate(sig_noise, cp.ones(128), mode='same') / 128

>>> import matplotlib.pyplot as plt
>>> clock = cp.arange(64, len(sig), 128)
>>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True)
>>> ax_orig.plot(cp.asnumpy(sig))
>>> ax_orig.plot(cp.asnumpy(clock), cp.asnumpy(sig[clock]), 'ro')
>>> ax_orig.set_title('Original signal')
>>> ax_noise.plot(cp.asnumpy(sig_noise))
>>> ax_noise.set_title('Signal with noise')
>>> ax_corr.plot(cp.asnumpy(corr))
>>> ax_corr.plot(cp.asnumpy(clock), cp.asnumpy(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()

cusignal.convolution.correlate.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.

Parameters:

in1array_like

First input.

in2array_like

Second input. Should have the same number of dimensions as in1.

modestr {‘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.

boundarystr {‘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.

fillvaluescalar, optional

Value to fill pad input arrays with. Default is 0.

Returns:

correlate2dndarray: A 2-dimensional array containing a subset of the discrete linear cross-correlation of in1 with in2.

Examples

Use 2D cross-correlation to find the location of a template in a noisy image: >>> import cusignal >>> import cupy as cp >>> from scipy import misc >>> face = cp.asarray(misc.face(gray=True) - misc.face(gray=True).mean()) >>> template = cp.copy(face[300:365, 670:750]) # right eye >>> template -= template.mean() >>> face = face + cp.random.randn(*face.shape) * 50 # add noise >>> corr = cusignal.correlate2d(face, template, boundary=’symm’, mode=’same’) >>> y, x = cp.unravel_index(cp.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(cp.asnumpy(face), cmap=’gray’) >>> ax_orig.set_title(‘Original’) >>> ax_orig.set_axis_off() >>> ax_template.imshow(cp.asnumpy(template), cmap=’gray’) >>> ax_template.set_title(‘Template’) >>> ax_template.set_axis_off() >>> ax_corr.imshow(cp.asnumpy(corr), cmap=’gray’) >>> ax_corr.set_title(‘Cross-correlation’) >>> ax_corr.set_axis_off() >>> ax_orig.plot(cp.asnumpy(x), cp.asnumpy(y), ‘ro’) >>> fig.show()

cusignal.convolution.correlate.correlation_lags(in1_len, in2_len, mode='full')#

Calculates the lag / displacement indices array for 1D cross-correlation. Parameters ———- in1_size : int

First input size.

in2_sizeint: Second input size.
modestr {‘full’, ‘valid’, ‘same’}, optional: A string indicating the size of the output. See the documentation correlate for more information.

Notes#

Cross-correlation for continuous functions \(f\) and \(g\) is defined as: .. math:

\left ( f\star g \right )\left ( \tau \right )
\triangleq \int_{t_0}^{t_0 +T}
\overline{f\left ( t \right )}g\left ( t+\tau \right )dt

Where \(\tau\) is defined as the displacement, also known as the lag. Cross correlation for discrete functions \(f\) and \(g\) is defined as: .. math:

\left ( f\star g \right )\left [ n \right ]
\triangleq \sum_{-\infty}^{\infty}
\overline{f\left [ m \right ]}g\left [ m+n \right ]

Where \(n\) is the lag. Examples ——– Cross-correlation of a signal with its time-delayed self. >>> import cusignal >>> import cupy as cp >>> from cupy.random import default_rng >>> rng = default_rng() >>> x = rng.standard_normal(1000) >>> y = cp.concatenate([rng.standard_normal(100), x]) >>> correlation = cusignal.correlate(x, y, mode=”full”) >>> lags = cusignal.correlation_lags(x.size, y.size, mode=”full”) >>> lag = lags[cp.argmax(correlation)]

Modulation/Demodulation#

Demodulation#

cusignal.demod.demod.fm_demod(x, axis=-1)#

Demodulate Frequency Modulated Signal

Parameters:

xndarray: Received complex valued signal or batch of signals

Returns:

yndarray: The demodulated output with the same shape as x.

Estimation#

Kalman Filter#

class cusignal.estimation.filters.KalmanFilter(dim_x, dim_z, dim_u=0, points=1, dtype=<class 'numpy.float32'>)#

This is a multi-point Kalman Filter implementation of rlabbe/filterpy, with a subset of functionality.

All Kalman Filter matrices are stack on the X axis. This is to allow for optimal global accesses on the GPU.

Parameters:

dim_xint: Number of state variables for the Kalman filter. For example, if you are tracking the position and velocity of an object in two dimensions, dim_x would be 4. This is used to set the default size of P, Q, and u
dim_zint: Number of of measurement inputs. For example, if the sensor provides you with position in (x,y), dim_z would be 2.
dim_uint (optional): Size of the control input, if it is being used. Default value of 0 indicates it is not used.
pointsint (optional): Number of Kalman Filter points to track.
dtypedtype (optional): Data type of compute.

References

[1]

Dan Simon. “Optimal State Estimation.” John Wiley & Sons. p. 208-212. (2006)

[2]

Roger Labbe. “Kalman and Bayesian Filters in Python” rlabbe/Kalman-and-Bayesian-Filters-in-Python

Examples

Here is a filter that tracks position and velocity using a sensor that only reads position.

First construct the object with the required dimensionality, number of points, and data type.

import cupy as cp
import numpy as np

from cusignal import KalmanFilter

points = 1024
kf = KalmanFilter(dim_x=4, dim_z=2, points=points, dtype=cp.float64)

Assign the initial value for the state (position and velocity) for all Kalman Filter points.

initial_location = np.array(
    [[10.0, 10.0, 0.0, 0.0]], dtype=dt
).T  # x, y, v_x, v_y
kf.x = cp.repeat(
    cp.asarray(initial_location[cp.newaxis, :, :]), points, axis=0
)

Define the state transition matrix for all Kalman Filter points:

F = np.array(
    [
        [1.0, 0.0, 1.0, 0.0],  # x = x0 + v_x*dt
        [0.0, 1.0, 0.0, 1.0],  # y = y0 + v_y*dt
        [0.0, 0.0, 1.0, 0.0],  # dx = v_x
        [1.0, 0.0, 0.0, 1.0],
    ],  # dy = v_y
    dtype=dt,
)
kf.F = cp.repeat(cp.asarray(F[cp.newaxis, :, :]), points, axis=0)

Define the measurement function for all Kalman Filter points:

H = np.array(
    [[1.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 1.0]],
    dtype=dt,  # x_0  # y_0
)
kf.H = cp.repeat(cp.asarray(H[cp.newaxis, :, :]), points, axis=0)

Define the covariance matrix for all Kalman Filter points:

initial_estimate_error = np.eye(dim_x, dtype=dt) * np.array(
    [1.0, 1.0, 2.0, 2.0], dtype=dt
)
kf.P = cp.repeat(
    cp.asarray(initial_estimate_error[cp.newaxis, :, :]),
    points,
    axis=0,
)

Define the measurement noise for all Kalman Filter points:

measurement_noise = np.eye(dim_z, dtype=dt) * 0.01
kf.R = cp.repeat(
    cp.asarray(measurement_noise[cp.newaxis, :, :]), points, axis=0
)

Define the process noise for all Kalman Filter points:

Now just perform the standard predict/update loop: Note: This example just uses the same sensor reading for all points

kf.predict()
z = get_sensor_reading() (dim_z, 1)
kf.z = cp.repeat(z[cp.newaxis, :, :], points, axis=0)
kf.update()

Results are in:

Attributes

x	(array(points, dim_x, 1)) Current state estimate. Any call to update() or predict() updates this variable.
P	(array(points, dim_x, dim_x)) Current state covariance matrix. Any call to update() or predict() updates this variable.
z	(array(points, dim_z, 1)) Last measurement used in update(). Read only.
R	(array(points, dim_z, dim_z)) Measurement noise matrix
Q	(array(points, dim_x, dim_x)) Process noise matrix
F	(array(points, dim_x, dim_x)) State Transition matrix
H	(array(points, dim_z, dim_x)) Measurement function
_alpha_sq	(float (points, 1, 1)) Fading memory setting. 1.0 gives the normal Kalman filter, and values slightly larger than 1.0 (such as 1.02) give a fading memory effect - previous measurements have less influence on the filter’s estimates. This formulation of the Fading memory filter (there are many) is due to Dan Simon [1].

Methods

`predict`([u, B, F, Q])	Predict next state (prior) using the Kalman filter state propagation equations.
`update`(z[, R, H])	Add a new measurement (z) to the Kalman filter.

predict(u=None, B=None, F=None, Q=None)#

Predict next state (prior) using the Kalman filter state propagation equations.

Parameters:

unarray, default 0: Optional control vector.
Barray(points, dim_x, dim_u), or None: Optional control transition matrix; a value of None will cause the filter to use self.B.
Farray(points, dim_x, dim_x), or None: Optional state transition matrix; a value of None will cause the filter to use self.F.
Qarray(points, dim_x, dim_x), scalar, or None: Optional process noise matrix; a value of None will cause the filter to use self.Q.

update(z, R=None, H=None)#

Add a new measurement (z) to the Kalman filter.

Parameters:

zarray(points, dim_z, 1): measurement for this update. z can be a scalar if dim_z is 1, otherwise it must be convertible to a column vector. If you pass in a value of H, z must be a column vector the of the correct size.
Rarray(points, dim_z, dim_z), scalar, or None: Optionally provide R to override the measurement noise for this one call, otherwise self.R will be used.
Harray(points, dim_z, dim_x), or None: Optionally provide H to override the measurement function for this one call, otherwise self.H will be used.

Filtering#

Resample#

cusignal.filtering.resample.decimate(x, q, n=None, axis=-1, zero_phase=True, gpupath=True)#

Downsample the signal after applying an anti-aliasing filter.

Parameters:

xarray_like: The signal to be downsampled, as an N-dimensional array.
qint: The downsampling factor.
nint or array_like, optional: The order of the filter (1 less than the length for FIR) to calculate, or the FIR filter coefficients to employ. Defaults to calculating the coefficients for 20 times the downsampling factor.
axisint, optional: The axis along which to decimate.
zero_phasebool, optional: Prevent 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.
gpupathbool, Optional: Optional path for filter design. gpupath == False may be desirable if filter sizes are small.

Returns:

yndarray: The down-sampled signal.

See also

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

Notes

Only FIR filter types are currently supported in cuSignal.

cusignal.filtering.resample.resample(x, num, t=None, axis=0, window=None, domain='time')#

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.

Parameters:

xarray_like

The data to be resampled.

numint

The number of samples in the resampled signal.

tarray_like, optional

If t is given, it is assumed to be the sample positions associated with the signal data in x.

axisint, optional

The axis of x that is resampled. Default is 0.

windowarray_like, callable, string, float, or tuple, optional

Specifies the window applied to the signal in the Fourier domain. See below for details.

domainstring, optional

A string indicating the domain of the input x:

time: Consider the input x as time-domain. (Default)
freq: Consider the input x as frequency-domain.

Returns:

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.

See also

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

Notes

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 cusignal.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.

Examples

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

>>> import cusignal
>>> import cupy as cp

>>> x = cp.linspace(0, 10, 20, endpoint=False)
>>> y = cp.cos(-x**2/6.0)
>>> f = cusignal.resample(y, 100)
>>> xnew = cp.linspace(0, 10, 100, endpoint=False)

>>> import matplotlib.pyplot as plt
>>> plt.plot(cp.asnumpy(x), cp.asnumpy(y), 'go-', cp.asnumpy(xnew),                 cp.asnumpy(f), '.-', 10, cp.asnumpy(y[0]), 'ro')
>>> plt.legend(['data', 'resampled'], loc='best')
>>> plt.show()

cusignal.filtering.resample.resample_poly(x, up, down, axis=0, window=('kaiser', 5.0), gpupath=True)#

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.

Parameters:

xarray_like: The data to be resampled.
upint: The upsampling factor.
downint: The downsampling factor.
axisint, optional: The axis of x that is resampled. Default is 0.
windowstring, 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.
gpupathbool, Optional: Optional path for filter design. gpupath == False may be desirable if filter sizes are small.

Returns:

resampled_xarray: The resampled array.

See also

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

Notes

This polyphase method will likely be faster than the Fourier method in cusignal.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 cusignal.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 cusignal.get_window and cusignal.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).

Examples

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:

>>> import cusignal
>>> import cupy as cp

>>> x = cp.linspace(0, 10, 20, endpoint=False)
>>> y = cp.cos(-x**2/6.0)
>>> f_fft = cusignal.resample(y, 100)
>>> f_poly = cusignal.resample_poly(y, 100, 20)
>>> xnew = cp.linspace(0, 10, 100, endpoint=False)

>>> import matplotlib.pyplot as plt
>>> plt.plot(cp.asnumpy(xnew), cp.asnumpy(f_fft), 'b.-',                  cp.asnumpy(xnew), cp.asnumpy(f_poly), 'r.-')
>>> plt.plot(cp.asnumpy(x), cp.asnumpy(y), 'ko-')
>>> plt.plot(10, cp.asnumpy(y[0]), 'bo', 10, 0., 'ro')  # boundaries
>>> plt.legend(['resample', 'resamp_poly', 'data'], loc='best')
>>> plt.show()

cusignal.filtering.resample.upfirdn(h, x, up=1, down=1, axis=-1)#

Upsample, FIR filter, and downsample.

Parameters:

harray_like: 1-dimensional FIR (finite-impulse response) filter coefficients.
xarray_like: Input signal array.
upint, optional: Upsampling rate. Default is 1.
downint, optional: Downsampling rate. Default is 1.
axisint, 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.

Returns:

yndarray: The output signal array. Dimensions will be the same as x except for along axis, which will change size according to the h, up, and down parameters.

Notes

The algorithm is an implementation of the block diagram shown on page 129 of the Vaidyanathan text [1] (Figure 4.3-8d).

The direct approach of upsampling by factor of P with zero insertion, FIR filtering of length N, and downsampling by factor of Q is O(N*Q) per output sample. The polyphase implementation used here is O(N/P).

References

[1]

P. P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice Hall, 1993.

Examples

Simple operations: >>> from cusignal import upfirdn >>> upfirdn([1, 1, 1], [1, 1, 1]) # FIR filter array([ 1., 2., 3., 2., 1.]) >>> upfirdn([1], [1, 2, 3], 3) # upsampling with zeros insertion array([ 1., 0., 0., 2., 0., 0., 3., 0., 0.]) >>> upfirdn([1, 1, 1], [1, 2, 3], 3) # upsampling with sample-and-hold array([ 1., 1., 1., 2., 2., 2., 3., 3., 3.]) >>> upfirdn([.5, 1, .5], [1, 1, 1], 2) # linear interpolation array([ 0.5, 1. , 1. , 1. , 1. , 1. , 0.5, 0. ]) >>> upfirdn([1], cp.arange(10), 1, 3) # decimation by 3 array([ 0., 3., 6., 9.]) >>> upfirdn([.5, 1, .5], cp.arange(10), 2, 3) # linear interp, rate 2/3 array([ 0. , 1. , 2.5, 4. , 5.5, 7. , 8.5, 0. ]) Apply a single filter to multiple signals: >>> x = cp.reshape(cp.arange(8), (4, 2)) >>> x array([[0, 1],

[2, 3], [4, 5], [6, 7]])

Apply along the last dimension of x: >>> h = [1, 1] >>> upfirdn(h, x, 2) array([[ 0., 0., 1., 1.],

[ 2., 2., 3., 3.], [ 4., 4., 5., 5.], [ 6., 6., 7., 7.]])

Apply along the 0th dimension of x: >>> upfirdn(h, x, 2, axis=0) array([[ 0., 1.],

[ 0., 1.], [ 2., 3.], [ 2., 3.], [ 4., 5.], [ 4., 5.], [ 6., 7.], [ 6., 7.]])

FIR Filters#

Filter data along one-dimension with an FIR filter.

Filter a data sequence, x, using a digital filter. This works for many fundamental data types (including Object type). Please note, cuSignal doesn’t support IIR filters presently, and this implementation is optimized for large filtering operations (and inherently depends on fftconvolve)

Parameters#

barray_like: The numerator coefficient vector in a 1-D sequence.
xarray_like: An N-dimensional input array.
axisint, 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.
ziarray_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.

Returns#

yarray: The output of the digital filter.
zfarray, optional: If zi is None, this is not returned, otherwise, zf holds the final filter delay values.

cusignal.filter_design.fir_filter_design.cmplx_sort(p)#

Sort roots based on magnitude.

Parameters:

parray_like: The roots to sort, as a 1-D array.

Returns:

p_sortedndarray: Sorted roots.
indxndarray: Array of indices needed to sort the input p.

Examples

>>> import cusignal
>>> vals = [1, 4, 1+1.j, 3]
>>> p_sorted, indx = cusignal.cmplx_sort(vals)
>>> p_sorted
array([1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j])
>>> indx
array([0, 2, 3, 1])

cusignal.filter_design.fir_filter_design.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None, gpupath=True)#

FIR filter design using the window method.

This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if numtaps is odd and Type II if numtaps is even.

Type II filters always have zero response at the Nyquist frequency, so a ValueError exception is raised if firwin is called with numtaps even and having a passband whose right end is at the Nyquist frequency.

Parameters:

numtapsint

Length of the filter (number of coefficients, i.e. the filter order + 1). numtaps must be odd if a passband includes the Nyquist frequency.

cutofffloat or 1D array_like

Cutoff frequency of filter (expressed in the same units as fs) OR an array of cutoff frequencies (that is, band edges). In the latter case, the frequencies in cutoff should be positive and monotonically increasing between 0 and fs/2. The values 0 and fs/2 must not be included in cutoff.

widthfloat or None, optional

If width is not None, then assume it is the approximate width of the transition region (expressed in the same units as fs) for use in Kaiser FIR filter design. In this case, the window argument is ignored.

windowstring or tuple of string and parameter values, optional

Desired window to use. See cusignal.get_window for a list of windows and required parameters.

pass_zero{True, False, ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’},

optional If True, the gain at the frequency 0 (i.e. the “DC gain”) is 1. If False, the DC gain is 0. Can also be a string argument for the desired filter type (equivalent to btype in IIR design functions).

New in version 1.3.0: Support for string arguments.

scalebool, optional

Set to True to scale the coefficients so that the frequency response is exactly unity at a certain frequency. That frequency is either:

0 (DC) if the first passband starts at 0 (i.e. pass_zero is True)
fs/2 (the Nyquist frequency) if the first passband ends at fs/2 (i.e the filter is a single band highpass filter); center of first passband otherwise

nyqfloat, optional

Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in cutoff must be between 0 and nyq. Default is 1.

fsfloat, optional

The sampling frequency of the signal. Each frequency in cutoff must be between 0 and fs/2. Default is 2.

gpupathbool, Optional

Optional path for filter design. gpupath == False may be desirable if filter sizes are small.

Returns:

h(numtaps,) ndarray: Coefficients of length numtaps FIR filter.

Raises:

ValueError: If any value in cutoff is less than or equal to 0 or greater than or equal to fs/2, if the values in cutoff are not strictly monotonically increasing, or if numtaps is even but a passband includes the Nyquist frequency.

See also

firwin2
firls
minimum_phase
remez

Examples

Low-pass from 0 to f:

>>> import cusignal
>>> numtaps = 3
>>> f = 0.1
>>> cusignal.firwin(numtaps, f)
array([ 0.06799017,  0.86401967,  0.06799017])

Use a specific window function:

>>> cusignal.firwin(numtaps, f, window='nuttall')
array([  3.56607041e-04,   9.99286786e-01,   3.56607041e-04])

High-pass (‘stop’ from 0 to f):

>>> cusignal.firwin(numtaps, f, pass_zero=False)
array([-0.00859313,  0.98281375, -0.00859313])

Band-pass:

>>> f1, f2 = 0.1, 0.2
>>> cusignal.firwin(numtaps, [f1, f2], pass_zero=False)
array([ 0.06301614,  0.88770441,  0.06301614])

Band-stop:

>>> cusignal.firwin(numtaps, [f1, f2])
array([-0.00801395,  1.0160279 , -0.00801395])

Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1]):

>>> f3, f4 = 0.3, 0.4
>>> cusignal.firwin(numtaps, [f1, f2, f3, f4])
array([-0.01376344,  1.02752689, -0.01376344])

Multi-band (passbands are [f1, f2] and [f3,f4]):

>>> cusignal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
array([ 0.04890915,  0.91284326,  0.04890915])

cusignal.filter_design.fir_filter_design.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None, antisymmetric=False, fs=None, gpupath=True)#

FIR filter design using the window method. From the given frequencies freq and corresponding gains gain, this function constructs an FIR filter with linear phase and (approximately) the given frequency response. Parameters ———- numtaps : int

The number of taps in the FIR filter. numtaps must be less than nfreqs.

freqarray_like, 1-D: The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being Nyquist. The Nyquist frequency is half fs. The values in freq must be nondecreasing. A value can be repeated once to implement a discontinuity. The first value in freq must be 0, and the last value must be fs/2. Values 0 and fs/2 must not be repeated.
gainarray_like: The filter gains at the frequency sampling points. Certain constraints to gain values, depending on the filter type, are applied, see Notes for details.
nfreqsint, optional: The size of the interpolation mesh used to construct the filter. For most efficient behavior, this should be a power of 2 plus 1 (e.g, 129, 257, etc). The default is one more than the smallest power of 2 that is not less than numtaps. nfreqs must be greater than numtaps.
windowstring or (string, float) or float, or None, optional: Window function to use. Default is “hamming”. See scipy.signal.get_window for the complete list of possible values. If None, no window function is applied.
nyqfloat, optional: Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in freq must be between 0 and nyq. Default is 1.
antisymmetricbool, optional: Whether resulting impulse response is symmetric/antisymmetric. See Notes for more details.
fsfloat, optional: The sampling frequency of the signal. Each frequency in cutoff must be between 0 and fs/2. Default is 2.

Returns#

tapsndarray: The filter coefficients of the FIR filter, as a 1-D array of length numtaps.

References#

[1]

Oppenheim, A. V. and Schafer, R. W., “Discrete-Time Signal Processing”, Prentice-Hall, Englewood Cliffs, New Jersey (1989). (See, for example, Section 7.4.)

[2]

Smith, Steven W., “The Scientist and Engineer’s Guide to Digital Signal Processing”, Ch. 17. http://www.dspguide.com/ch17/1.htm

cusignal.filter_design.fir_filter_design.kaiser_atten(numtaps, width)#

Compute the attenuation of a Kaiser FIR filter. Given the number of taps N and the transition width width, compute the attenuation a in dB, given by Kaiser’s formula:

a = 2.285 * (N - 1) * pi * width + 7.95

Parameters#

numtapsint: The number of taps in the FIR filter.
widthfloat: The desired width of the transition region between passband and stopband (or, in general, at any discontinuity) for the filter.

Returns#

afloat: The attenuation of the ripple, in dB.

cusignal.filter_design.fir_filter_design.kaiser_beta(a)#

Compute the Kaiser parameter beta, given the attenuation a. Parameters ———- a : float

The desired attenuation in the stopband and maximum ripple in the passband, in dB. This should be a positive number.

Returns#

betafloat: The beta parameter to be used in the formula for a Kaiser window.

References#

Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.

Channelizer#

Polyphase channelize signal into n channels

Parameters#

xarray_like: The input data to be channelized
harray_like: The 1-D input filter; will be split into n channels of int number of taps
n_chansint: Number of channels for channelizer

Returns#

yy : channelized output matrix

Notes#

Currently only supports simple channelizer where channel spacing is equivalent to the number of channels used (zero overlap). Number of filter taps (len of filter / n_chans) must be <=32.

Filter Design#

Resample#

cusignal.filter_design.fir_filter_design.cmplx_sort(p)#

Sort roots based on magnitude.

Parameters:

parray_like: The roots to sort, as a 1-D array.

Returns:

p_sortedndarray: Sorted roots.
indxndarray: Array of indices needed to sort the input p.

Examples

>>> import cusignal
>>> vals = [1, 4, 1+1.j, 3]
>>> p_sorted, indx = cusignal.cmplx_sort(vals)
>>> p_sorted
array([1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j])
>>> indx
array([0, 2, 3, 1])

cusignal.filter_design.fir_filter_design.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None, gpupath=True)#

FIR filter design using the window method.

This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if numtaps is odd and Type II if numtaps is even.

Parameters:

numtapsint

Length of the filter (number of coefficients, i.e. the filter order + 1). numtaps must be odd if a passband includes the Nyquist frequency.

cutofffloat or 1D array_like

widthfloat or None, optional

windowstring or tuple of string and parameter values, optional

Desired window to use. See cusignal.get_window for a list of windows and required parameters.

pass_zero{True, False, ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’},

New in version 1.3.0: Support for string arguments.

scalebool, optional

Set to True to scale the coefficients so that the frequency response is exactly unity at a certain frequency. That frequency is either:

0 (DC) if the first passband starts at 0 (i.e. pass_zero is True)
fs/2 (the Nyquist frequency) if the first passband ends at fs/2 (i.e the filter is a single band highpass filter); center of first passband otherwise

nyqfloat, optional

Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in cutoff must be between 0 and nyq. Default is 1.

fsfloat, optional

The sampling frequency of the signal. Each frequency in cutoff must be between 0 and fs/2. Default is 2.

gpupathbool, Optional

Optional path for filter design. gpupath == False may be desirable if filter sizes are small.

Returns:

h(numtaps,) ndarray: Coefficients of length numtaps FIR filter.

Raises:

ValueError: If any value in cutoff is less than or equal to 0 or greater than or equal to fs/2, if the values in cutoff are not strictly monotonically increasing, or if numtaps is even but a passband includes the Nyquist frequency.

See also

firwin2
firls
minimum_phase
remez

Examples

Low-pass from 0 to f:

>>> import cusignal
>>> numtaps = 3
>>> f = 0.1
>>> cusignal.firwin(numtaps, f)
array([ 0.06799017,  0.86401967,  0.06799017])

Use a specific window function:

>>> cusignal.firwin(numtaps, f, window='nuttall')
array([  3.56607041e-04,   9.99286786e-01,   3.56607041e-04])

High-pass (‘stop’ from 0 to f):

>>> cusignal.firwin(numtaps, f, pass_zero=False)
array([-0.00859313,  0.98281375, -0.00859313])

Band-pass:

>>> f1, f2 = 0.1, 0.2
>>> cusignal.firwin(numtaps, [f1, f2], pass_zero=False)
array([ 0.06301614,  0.88770441,  0.06301614])

Band-stop:

>>> cusignal.firwin(numtaps, [f1, f2])
array([-0.00801395,  1.0160279 , -0.00801395])

Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1]):

>>> f3, f4 = 0.3, 0.4
>>> cusignal.firwin(numtaps, [f1, f2, f3, f4])
array([-0.01376344,  1.02752689, -0.01376344])

Multi-band (passbands are [f1, f2] and [f3,f4]):

>>> cusignal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
array([ 0.04890915,  0.91284326,  0.04890915])

cusignal.filter_design.fir_filter_design.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None, antisymmetric=False, fs=None, gpupath=True)#

The number of taps in the FIR filter. numtaps must be less than nfreqs.

freqarray_like, 1-D: The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being Nyquist. The Nyquist frequency is half fs. The values in freq must be nondecreasing. A value can be repeated once to implement a discontinuity. The first value in freq must be 0, and the last value must be fs/2. Values 0 and fs/2 must not be repeated.
gainarray_like: The filter gains at the frequency sampling points. Certain constraints to gain values, depending on the filter type, are applied, see Notes for details.
nfreqsint, optional: The size of the interpolation mesh used to construct the filter. For most efficient behavior, this should be a power of 2 plus 1 (e.g, 129, 257, etc). The default is one more than the smallest power of 2 that is not less than numtaps. nfreqs must be greater than numtaps.
windowstring or (string, float) or float, or None, optional: Window function to use. Default is “hamming”. See scipy.signal.get_window for the complete list of possible values. If None, no window function is applied.
nyqfloat, optional: Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in freq must be between 0 and nyq. Default is 1.
antisymmetricbool, optional: Whether resulting impulse response is symmetric/antisymmetric. See Notes for more details.
fsfloat, optional: The sampling frequency of the signal. Each frequency in cutoff must be between 0 and fs/2. Default is 2.

Returns#

tapsndarray: The filter coefficients of the FIR filter, as a 1-D array of length numtaps.

References#

[1]

Oppenheim, A. V. and Schafer, R. W., “Discrete-Time Signal Processing”, Prentice-Hall, Englewood Cliffs, New Jersey (1989). (See, for example, Section 7.4.)

[2]

Smith, Steven W., “The Scientist and Engineer’s Guide to Digital Signal Processing”, Ch. 17. http://www.dspguide.com/ch17/1.htm

cusignal.filter_design.fir_filter_design.kaiser_atten(numtaps, width)#

Compute the attenuation of a Kaiser FIR filter. Given the number of taps N and the transition width width, compute the attenuation a in dB, given by Kaiser’s formula:

a = 2.285 * (N - 1) * pi * width + 7.95

Parameters#

numtapsint: The number of taps in the FIR filter.
widthfloat: The desired width of the transition region between passband and stopband (or, in general, at any discontinuity) for the filter.

Returns#

afloat: The attenuation of the ripple, in dB.

cusignal.filter_design.fir_filter_design.kaiser_beta(a)#

Compute the Kaiser parameter beta, given the attenuation a. Parameters ———- a : float

The desired attenuation in the stopband and maximum ripple in the passband, in dB. This should be a positive number.

Returns#

betafloat: The beta parameter to be used in the formula for a Kaiser window.

References#

Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.

Peak Finding#

cusignal.peak_finding.peak_finding.argrelextrema(data, comparator, axis=0, order=1, mode='clip')#

Calculate the relative extrema of data.

Parameters:

datandarray: Array in which to find the relative extrema.
comparatorcallable: Function to use to compare two data points. Should take two arrays as arguments.
axisint, optional: Axis over which to select from data. Default is 0.
orderint, optional: How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.

Returns:

extrematuple of ndarrays: Indices of the maxima in arrays of integers. extrema[k] is the array of indices of axis k of data. Note that the return value is a tuple even when data is one-dimensional.

See also

argrelmin, argrelmax

Examples

>>> from cusignal import argrelextrema
>>> import cupy as cp
>>> x = cp.array([2, 1, 2, 3, 2, 0, 1, 0])
>>> argrelextrema(x, cp.greater)
(array([3, 6]),)
>>> y = cp.array([[1, 2, 1, 2],
...               [2, 2, 0, 0],
...               [5, 3, 4, 4]])
...
>>> argrelextrema(y, cp.less, axis=1)
(array([0, 2]), array([2, 1]))

cusignal.peak_finding.peak_finding.argrelmax(data, axis=0, order=1, mode='clip')#

Calculate the relative maxima of data.

Parameters:

datandarray: Array in which to find the relative maxima.
axisint, optional: Axis over which to select from data. Default is 0.
orderint, optional: How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.

Returns:

extrematuple of ndarrays: Indices of the maxima in arrays of integers. extrema[k] is the array of indices of axis k of data. Note that the return value is a tuple even when data is one-dimensional.

See also

argrelextrema, argrelmin, find_peaks

Notes

This function uses argrelextrema with cp.greater as comparator. Therefore it requires a strict inequality on both sides of a value to consider it a maximum. This means flat maxima (more than one sample wide) are not detected. In case of one-dimensional data find_peaks can be used to detect all local maxima, including flat ones.

Examples

>>> from cusignal import argrelmax
>>> import cupy as cp
>>> x = cp.array([2, 1, 2, 3, 2, 0, 1, 0])
>>> argrelmax(x)
(array([3, 6]),)
>>> y = cp.array([[1, 2, 1, 2],
...               [2, 2, 0, 0],
...               [5, 3, 4, 4]])
...
>>> argrelmax(y, axis=1)
(array([0]), array([1]))

cusignal.peak_finding.peak_finding.argrelmin(data, axis=0, order=1, mode='clip')#

Calculate the relative minima of data.

Parameters:

datandarray: Array in which to find the relative minima.
axisint, optional: Axis over which to select from data. Default is 0.
orderint, optional: How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.

Returns:

extrematuple of ndarrays: Indices of the minima in arrays of integers. extrema[k] is the array of indices of axis k of data. Note that the return value is a tuple even when data is one-dimensional.

See also

argrelextrema, argrelmax, find_peaks

Notes

This function uses argrelextrema with cp.less as comparator. Therefore it requires a strict inequality on both sides of a value to consider it a minimum. This means flat minima (more than one sample wide) are not detected. In case of one-dimensional data find_peaks can be used to detect all local minima, including flat ones, by calling it with negated data.

Examples

>>> from cusignal import argrelmin
>>> import cupy as cp
>>> x = cp.array([2, 1, 2, 3, 2, 0, 1, 0])
>>> argrelmin(x)
(array([1, 5]),)
>>> y = cp.array([[1, 2, 1, 2],
...               [2, 2, 0, 0],
...               [5, 3, 4, 4]])
...
>>> argrelmin(y, axis=1)
(array([0, 2]), array([2, 1]))

Window Functions#

Windows#

cusignal.windows.windows.barthann(M, sym=True)#

Return a modified Bartlett-Hann window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.barthann(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Bartlett-Hann window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Bartlett-Hann window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.bartlett(M, sym=True)#

Return a Bartlett window.

The Bartlett window is very similar to a triangular window, except that the end points are at zero. It is often used in signal processing for tapering a signal, without generating too much ripple in the frequency domain.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The triangular window, with the first and last samples equal to zero and the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

See also

triang: A triangular window that does not touch zero at the ends

Notes

The Bartlett window is defined as

\[w(n) = \frac{2}{M-1} \left( \frac{M-1}{2} - \left|n - \frac{M-1}{2}\right| \right)\]

Most references to the Bartlett window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. Note that convolution with this window produces linear interpolation. It is also known as an apodization (which means”removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function. The Fourier transform of the Bartlett is the product of two sinc functions. Note the excellent discussion in Kanasewich. [2]

References

[1]

M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika 37, 1-16, 1950.

[2]

E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 109-110.

[3]

A.V. Oppenheim and R.W. Schafer, “Discrete-Time Signal Processing”, Prentice-Hall, 1999, pp. 468-471.

[4]

Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function

[5]

W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, “Numerical Recipes”, Cambridge University Press, 1986, page 429.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.bartlett(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Bartlett window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Bartlett window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.blackman(M, sym=True)#

Return a Blackman window.

The Blackman window is a taper formed by using the first three terms of a summation of cosines. It was designed to have close to the minimal leakage possible. It is close to optimal, only slightly worse than a Kaiser window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The Blackman window is defined as

\[w(n) = 0.42 - 0.5 \cos(2\pi n/M) + 0.08 \cos(4\pi n/M)\]

The “exact Blackman” window was designed to null out the third and fourth sidelobes, but has discontinuities at the boundaries, resulting in a 6 dB/oct fall-off. This window is an approximation of the “exact” window, which does not null the sidelobes as well, but is smooth at the edges, improving the fall-off rate to 18 dB/oct. [3]

Most references to the Blackman window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function. It is known as a “near optimal” tapering function, almost as good (by some measures) as the Kaiser window.

References

[1]

Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.

[2]

Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.

[3]

Harris, Fredric J. (Jan 1978). “On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform”. Proceedings of the IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.blackman(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Blackman window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(cp.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Blackman window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.blackmanharris(M, sym=True)#

Return a minimum 4-term Blackman-Harris window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.blackmanharris(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Blackman-Harris window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Blackman-Harris window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.bohman(M, sym=True)#

Return a Bohman window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.bohman(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Bohman window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Bohman window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.boxcar(M, sym=True)#

Return a boxcar or rectangular window.

Also known as a rectangular window or Dirichlet window, this is equivalent to no window at all.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: Whether the window is symmetric. (Has no effect for boxcar.)

Returns:

wndarray: The window, with the maximum value normalized to 1.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.boxcar(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Boxcar window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the boxcar window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.chebwin(M, at, sym=True)#

Return a Dolph-Chebyshev window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
atfloat: Attenuation (in dB).
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value always normalized to 1

Notes

This window optimizes for the narrowest main lobe width for a given order M and sidelobe equiripple attenuation at, using Chebyshev polynomials. It was originally developed by Dolph to optimize the directionality of radio antenna arrays.

Unlike most windows, the Dolph-Chebyshev is defined in terms of its frequency response:

\[W(k) = \frac {\cos\{M \cos^{-1}[\beta \cos(\frac{\pi k}{M})]\}} {\cosh[M \cosh^{-1}(\beta)]}\]

where

\[\beta = \cosh \left [\frac{1}{M} \cosh^{-1}(10^\frac{A}{20}) \right ]\]

and 0 <= abs(k) <= M-1. A is the attenuation in decibels (at).

The time domain window is then generated using the IFFT, so power-of-two M are the fastest to generate, and prime number M are the slowest.

The equiripple condition in the frequency domain creates impulses in the time domain, which appear at the ends of the window.

References

[1]

C. Dolph, “A current distribution for broadside arrays which optimizes the relationship between beam width and side-lobe level”, Proceedings of the IEEE, Vol. 34, Issue 6

[2]

Peter Lynch, “The Dolph-Chebyshev Window: A Simple Optimal Filter”, American Meteorological Society (April 1997) http://mathsci.ucd.ie/~plynch/Publications/Dolph.pdf

[3]

F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transforms”, Proceedings of the IEEE, Vol. 66, No. 1, January 1978

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.chebwin(51, at=100)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Dolph-Chebyshev window (100 dB)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Dolph-Chebyshev window (100 dB)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.cosine(M, sym=True)#

Return a window with a simple cosine shape.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

New in version 0.13.0.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.cosine(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Cosine window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the cosine window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
>>> plt.show()

cusignal.windows.windows.exponential(M, center=None, tau=1.0, sym=True)#

Return an exponential (or Poisson) window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
centerfloat, optional: Parameter defining the center location of the window function. The default value if not given is center = (M-1) / 2. This parameter must take its default value for symmetric windows.
taufloat, optional: Parameter defining the decay. For center = 0 use tau = -(M-1) / ln(x) if x is the fraction of the window remaining at the end.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The Exponential window is defined as

\[w(n) = e^{-|n-center| / \tau}\]

References

S. Gade and H. Herlufsen, “Windows to FFT analysis (Part I)”, Technical Review 3, Bruel & Kjaer, 1987.

Examples

Plot the symmetric window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> M = 51
>>> tau = 3.0
>>> window = cusignal.exponential(M, tau=tau)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Exponential Window (tau=3.0)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -35, 0])
>>> plt.title("Frequency response of the Exponential window (tau=3.0)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

This function can also generate non-symmetric windows:

>>> tau2 = -(M-1) / np.log(0.01)
>>> window2 = cusignal.exponential(M, 0, tau2, False)
>>> plt.figure()
>>> plt.plot(cp.asnumpy(window2))
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

cusignal.windows.windows.flattop(M, sym=True)#

Return a flat top window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

Flat top windows are used for taking accurate measurements of signal amplitude in the frequency domain, with minimal scalloping error from the center of a frequency bin to its edges, compared to others. This is a 5th-order cosine window, with the 5 terms optimized to make the main lobe maximally flat. [1]

References

[1]

D’Antona, Gabriele, and A. Ferrero, “Digital Signal Processing for Measurement Systems”, Springer Media, 2006, p. 70 :doi:`10.1007/0-387-28666-7`.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.flattop(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Flat top window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the flat top window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.gaussian(M, std, sym=True)#

Return a Gaussian window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
stdfloat: The standard deviation, sigma.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The Gaussian window is defined as

\[w(n) = e^{ -\frac{1}{2}\left(\frac{n}{\sigma}\right)^2 }\]

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.gaussian(51, std=7)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title(r"Gaussian window ($\sigma$=7)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title(r"Frequency response of the Gaussian window ($\sigma$=7)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.general_cosine(M, a, sym=True)#

Generic weighted sum of cosine terms window

Parameters:

Mint: Number of points in the output window
aarray_like: Sequence of weighting coefficients. This uses the convention of being centered on the origin, so these will typically all be positive numbers, not alternating sign.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

References

[1]

A. Nuttall, “Some windows with very good sidelobe behavior,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.

[2]

Heinzel G. et al., “Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows”, February 15, 2002 https://holometer.fnal.gov/GH_FFT.pdf

Examples

Heinzel describes a flat-top window named “HFT90D” with formula: [2]

\[w_j = 1 - 1.942604 \cos(z) + 1.340318 \cos(2z) - 0.440811 \cos(3z) + 0.043097 \cos(4z)\]

where

\[z = \frac{2 \pi j}{N}, j = 0...N - 1\]

Since this uses the convention of starting at the origin, to reproduce the window, we need to convert every other coefficient to a positive number:

>>> HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097]

The paper states that the highest sidelobe is at -90.2 dB. Reproduce Figure 42 by plotting the window and its frequency response, and confirm the sidelobe level in red:

>>> from cusignal.windows import general_cosine
>>> from cupy.fft import fft, fftshift
>>> import cupy as cp
>>> import matplotlib.pyplot as plt

>>> window = general_cosine(1000, HFT90D, sym=False)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("HFT90D window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 10000) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(cp.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-50/1000, 50/1000, -140, 0])
>>> plt.title("Frequency response of the HFT90D window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
>>> plt.axhline(-90.2, color='red')
>>> plt.show()

cusignal.windows.windows.general_gaussian(M, p, sig, sym=True)#

Return a window with a generalized Gaussian shape.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
pfloat: Shape parameter. p = 1 is identical to gaussian, p = 0.5 is the same shape as the Laplace distribution.
sigfloat: The standard deviation, sigma.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The generalized Gaussian window is defined as

\[w(n) = e^{ -\frac{1}{2}\left|\frac{n}{\sigma}\right|^{2p} }\]

the half-power point is at

\[(2 \log(2))^{1/(2 p)} \sigma\]

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.general_gaussian(51, p=1.5, sig=7)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title(r"Generalized Gaussian window (p=1.5, $\sigma$=7)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title(r"Freq. resp. of the gen. Gaussian "
...           r"window (p=1.5, $\sigma$=7)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.general_hamming(M, alpha, sym=True)#

Return a generalized Hamming window.

The generalized Hamming window is constructed by multiplying a rectangular window by one period of a cosine function [1].

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
alphafloat: The window coefficient, \(\alpha\)
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

See also

hamming, hann

Notes

The generalized Hamming window is defined as

\[w(n) = \alpha - \left(1 - \alpha\right) \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1\]

Both the common Hamming window and Hann window are special cases of the generalized Hamming window with \(\alpha\) = 0.54 and \(\alpha\) = 0.5, respectively [2].

References

[1]

DSPRelated, “Generalized Hamming Window Family”, https://www.dsprelated.com/freebooks/sasp/Generalized_Hamming_Window_Family.html

[2]

Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function

[3]

Riccardo Piantanida ESA, “Sentinel-1 Level 1 Detailed Algorithm Definition”, https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Level-1-Detailed-Algorithm-Definition

[4]

Matthieu Bourbigot ESA, “Sentinel-1 Product Definition”, https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Product-Definition

Examples

The Sentinel-1A/B Instrument Processing Facility uses generalized Hamming windows in the processing of spaceborne Synthetic Aperture Radar (SAR) data [3]. The facility uses various values for the \(\alpha\) parameter based on operating mode of the SAR instrument. Some common \(\alpha\) values include 0.75, 0.7 and 0.52 [4]. As an example, we plot these different windows.

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> fig1, spatial_plot = plt.subplots()
>>> spatial_plot.set_title("Generalized Hamming Windows")
>>> spatial_plot.set_ylabel("Amplitude")
>>> spatial_plot.set_xlabel("Sample")

>>> fig2, freq_plot = plt.subplots()
>>> freq_plot.set_title("Frequency Responses")
>>> freq_plot.set_ylabel("Normalized magnitude [dB]")
>>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]")

>>> for alpha in [0.75, 0.7, 0.52]:
...     window = cusignal.general_hamming(41, alpha)
...     spatial_plot.plot(cp.asnumpy(window), label="{:.2f}".format(alpha))
...     A = fft(window, 2048) / (len(window)/2.0)
...     freq = cp.linspace(-0.5, 0.5, len(A))
...     response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
...     freq_plot.plot(
...         cp.asnumpy(freq), cp.asnumpy(response),
...         label="{:.2f}".format(alpha)
...     )
>>> freq_plot.legend(loc="upper right")
>>> spatial_plot.legend(loc="upper right")

cusignal.windows.windows.get_window(window, Nx, fftbins=True)#

Return a window of a given length and type.

Parameters:

windowstring, float, or tuple: The type of window to create. See below for more details.
Nxint: The number of samples in the window.
fftbinsbool, optional: If True (default), create a “periodic” window, ready to use with ifftshift and be multiplied by the result of an FFT (see also fftpack.fftfreq). If False, create a “symmetric” window, for use in filter design.

Returns:

get_windowndarray: Returns a window of length Nx and type window

Notes

Window types:

~cusignal.windows.windows.boxcar
~cusignal.windows.windows.triang
~cusignal.windows.windows.blackman
~cusignal.windows.windows.hamming
~cusignal.windows.windows.hann
~cusignal.windows.windows.bartlett
~cusignal.windows.windows.flattop
~cusignal.windows.windows.parzen
~cusignal.windows.windows.bohman
~cusignal.windows.windows.blackmanharris
~cusignal.windows.windows.nuttall
~cusignal.windows.windows.barthann
~cusignal.windows.windows.kaiser (needs beta)
~cusignal.windows.windows.gaussian (needs standard deviation)
~cusignal.windows.windows.general_gaussian
(needs power, width)
~cusignal.windows.windows.slepian (needs width)
~cusignal.windows.windows.dpss
(needs normalized half-bandwidth)
~cusignal.windows.windows.chebwin (needs attenuation)
~cusignal.windows.windows.exponential (needs decay scale)
~cusignal.windows.windows.tukey (needs taper fraction)

If the window requires no parameters, then window can be a string.

If the window requires parameters, then window must be a tuple with the first argument the string name of the window, and the next arguments the needed parameters.

If window is a floating point number, it is interpreted as the beta parameter of the ~cusignal.windows.windows.kaiser window.

Each of the window types listed above is also the name of a function that can be called directly to create a window of that type.

Examples

>>> import cusignal
>>> cusignal.get_window('triang', 7)
array([ 0.125,  0.375,  0.625,  0.875,  0.875,  0.625,  0.375])
>>> cusignal.get_window(('kaiser', 4.0), 9)
array([0.08848053, 0.32578323, 0.63343178, 0.89640418, 1.,
       0.89640418, 0.63343178, 0.32578323, 0.08848053])
>>> cusignal.get_window(4.0, 9)
array([0.08848053, 0.32578323, 0.63343178, 0.89640418, 1.,
       0.89640418, 0.63343178, 0.32578323, 0.08848053])

cusignal.windows.windows.hamming(M, sym=True)#

Return a Hamming window.

The Hamming window is a taper formed by using a raised cosine with non-zero endpoints, optimized to minimize the nearest side lobe.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The Hamming window is defined as

\[w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1\]

The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and is described in Blackman and Tukey. It was recommended for smoothing the truncated autocovariance function in the time domain. Most references to the Hamming window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function.

References

[1]

Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.

[2]

E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 109-110.

[3]

Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function

[4]

W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, “Numerical Recipes”, Cambridge University Press, 1986, page 425.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.hamming(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Hamming window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Hamming window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.hann(M, sym=True)#

Return a Hann window.

The Hann window is a taper formed by using a raised cosine or sine-squared with ends that touch zero.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The Hann window is defined as

\[w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1\]

The window was named for Julius von Hann, an Austrian meteorologist. It is also known as the Cosine Bell. It is sometimes erroneously referred to as the “Hanning” window, from the use of “hann” as a verb in the original paper and confusion with the very similar Hamming window.

Most references to the Hann window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function.

References

[1]

Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.

[2]

E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 106-108.

[3]

Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function

[4]

W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, “Numerical Recipes”, Cambridge University Press, 1986, page 425.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.hann(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Hann window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(np.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Hann window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.kaiser(M, beta, sym=True)#

Return a Kaiser window.

The Kaiser window is a taper formed by using a Bessel function.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
betafloat: Shape parameter, determines trade-off between main-lobe width and side lobe level. As beta gets large, the window narrows.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Notes

The Kaiser window is defined as

\[w(n) = I_0\left( \beta \sqrt{1-\frac{4n^2}{(M-1)^2}} \right)/I_0(\beta)\]

with

\[\quad -\frac{M-1}{2} \leq n \leq \frac{M-1}{2},\]

where \(I_0\) is the modified zeroth-order Bessel function.

The Kaiser was named for Jim Kaiser, who discovered a simple approximation to the DPSS window based on Bessel functions. The Kaiser window is a very good approximation to the Digital Prolate Spheroidal Sequence, or Slepian window, which is the transform which maximizes the energy in the main lobe of the window relative to total energy.

The Kaiser can approximate other windows by varying the beta parameter. (Some literature uses alpha = beta/pi.) [4]

beta	Window shape
0	Rectangular
5	Similar to a Hamming
6	Similar to a Hann
8.6	Similar to a Blackman

A beta value of 14 is probably a good starting point. Note that as beta gets large, the window narrows, and so the number of samples needs to be large enough to sample the increasingly narrow spike, otherwise NaNs will be returned.

Most references to the Kaiser window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function.

References

[1]

J. F. Kaiser, “Digital Filters” - Ch 7 in “Systems analysis by digital computer”, Editors: F.F. Kuo and J.F. Kaiser, p 218-285. John Wiley and Sons, New York, (1966).

[2]

E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 177-178.

[3]

Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function

[4]

F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978. :doi:`10.1109/PROC.1978.10837`.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.kaiser(51, beta=14)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title(r"Kaiser window ($\beta$=14)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title(r"Frequency response of the Kaiser window ($\beta$=14)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.nuttall(M, sym=True)#

Return a minimum 4-term Blackman-Harris window according to Nuttall.

This variation is called “Nuttall4c” by Heinzel. [2]

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

References

[1]

A. Nuttall, “Some windows with very good sidelobe behavior,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.

[2]

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.nuttall(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Nuttall window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Nuttall window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.parzen(M, sym=True)#

cusignal.windows.windows.taylor(M, nbar=4, sll=30, norm=True, sym=True)#

Return a Taylor window. The Taylor window taper function approximates the Dolph-Chebyshev window’s constant sidelobe level for a parameterized number of near-in sidelobes, but then allows a taper beyond [2]. The SAR (synthetic aperature radar) community commonly uses Taylor weighting for image formation processing because it provides strong, selectable sidelobe suppression with minimum broadening of the mainlobe [Rc5e098a1a171-1]. Parameters ———- M : int

Number of points in the output window. If zero or less, an empty array is returned.

nbarint, optional: Number of nearly constant level sidelobes adjacent to the mainlobe.
sllfloat, optional: Desired suppression of sidelobe level in decibels (dB) relative to the DC gain of the mainlobe. This should be a positive number.
normbool, optional: When True (default), divides the window by the largest (middle) value for odd-length windows or the value that would occur between the two repeated middle values for even-length windows such that all values are less than or equal to 1. When False the DC gain will remain at 1 (0 dB) and the sidelobes will be sll dB down.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns#

outarray: The window. When norm is True (default), the maximum value is normalized to 1 (though the value 1 does not appear if M is even and sym is True).

Examples#

Plot the window and its frequency response: >>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt >>> window = signal.windows.taylor(51, nbar=20, sll=100, norm=False) >>> plt.plot(window) >>> plt.title(“Taylor window (100 dB)”) >>> plt.ylabel(“Amplitude”) >>> plt.xlabel(“Sample”) >>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title(“Frequency response of the Taylor window (100 dB)”) >>> plt.ylabel(“Normalized magnitude [dB]”) >>> plt.xlabel(“Normalized frequency [cycles per sample]”)

cusignal.windows.windows.triang(M, sym=True)#

Return a triangular window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

See also

bartlett: A triangular window that touches zero

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.triang(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Triangular window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(cp.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the triangular window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

cusignal.windows.windows.tukey(M, alpha=0.5, sym=True)#

Return a Tukey window, also known as a tapered cosine window.

Parameters:

Mint: Number of points in the output window. If zero or less, an empty array is returned.
alphafloat, optional: Shape parameter of the Tukey window, representing the fraction of the window inside the cosine tapered region. If zero, the Tukey window is equivalent to a rectangular window. If one, the Tukey window is equivalent to a Hann window.
symbool, optional: When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns:

wndarray: The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).

References

[1]

Harris, Fredric J. (Jan 1978). “On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform”. Proceedings of the IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`

[2]

Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function#Tukey_window

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt

>>> window = cusignal.tukey(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Tukey window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.ylim([0, 1.1])

>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Tukey window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

Waveforms#

Waveform Generation#

cusignal.waveforms.waveforms.chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True, type='real')#

Frequency-swept cosine generator.

In the following, ‘Hz’ should be interpreted as ‘cycles per unit’; there is no requirement here that the unit is one second. The important distinction is that the units of rotation are cycles, not radians. Likewise, t could be a measurement of space instead of time.

Parameters:

tarray_like: Times at which to evaluate the waveform.
f0float: Frequency (e.g. Hz) at time t=0.
t1float: Time at which f1 is specified.
f1float: Frequency (e.g. Hz) of the waveform at time t1.
method{‘linear’, ‘quadratic’, ‘logarithmic’, ‘hyperbolic’}, optional: Kind of frequency sweep. If not given, linear is assumed. See Notes below for more details.
phifloat, optional: Phase offset, in degrees. Default is 0.
vertex_zerobool, optional: This parameter is only used when method is ‘quadratic’. It determines whether the vertex of the parabola that is the graph of the frequency is at t=0 or t=t1.
type{‘real’, ‘complex’}, optional: Specify output chirp type, only applicable when method is ‘linear’.

Returns:

yndarray: A numpy array containing the signal evaluated at t with the requested time-varying frequency. More precisely, the function returns cos(phase + (pi/180)*phi) where phase is the integral (from 0 to t) of 2*pi*f(t). f(t) is defined below.

Examples

The following will be used in the examples:

>>> from cusignal import chirp, spectrogram
>>> import matplotlib.pyplot as plt
>>> import cupy as cp

For the first example, we’ll plot the waveform for a linear chirp from 6 Hz to 1 Hz over 10 seconds:

>>> t = cp.linspace(0, 10, 5001)
>>> w = chirp(t, f0=6, f1=1, t1=10, method='linear')
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(w))
>>> plt.title("Linear Chirp, f(0)=6, f(10)=1")
>>> plt.xlabel('t (sec)')
>>> plt.show()

For the remaining examples, we’ll use higher frequency ranges, and demonstrate the result using cusignal.spectrogram. We’ll use a 10 second interval sampled at 8000 Hz.

>>> fs = 8000
>>> T = 10
>>> t = cp.linspace(0, T, T*fs, endpoint=False)

Quadratic chirp from 1500 Hz to 250 Hz over 10 seconds (vertex of the parabolic curve of the frequency is at t=0):

>>> w = chirp(t, f0=1500, f1=250, t1=10, method='quadratic')
>>> ff, tt, Sxx = spectrogram(w, fs=fs, noverlap=256, nperseg=512,
...                           nfft=2048)
>>> plt.pcolormesh(cp.asnumpy(tt), cp.asnumpy(ff[:513]),
                   cp.asnumpy(Sxx[:513]), cmap='gray_r')
>>> plt.title('Quadratic Chirp, f(0)=1500, f(10)=250')
>>> plt.xlabel('t (sec)')
>>> plt.ylabel('Frequency (Hz)')
>>> plt.grid()
>>> plt.show()

cusignal.waveforms.waveforms.gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False, retenv=False)#

Return a Gaussian modulated sinusoid:

exp(-a t^2) exp(1j*2*pi*fc*t).

If retquad is True, then return the real and imaginary parts (in-phase and quadrature). If retenv is True, then return the envelope (unmodulated signal). Otherwise, return the real part of the modulated sinusoid.

Parameters:

tndarray or the string ‘cutoff’: Input array.
fcint, optional: Center frequency (e.g. Hz). Default is 1000.
bwfloat, optional: Fractional bandwidth in frequency domain of pulse (e.g. Hz). Default is 0.5.
bwrfloat, optional: Reference level at which fractional bandwidth is calculated (dB). Default is -6.
tprfloat, optional: If t is ‘cutoff’, then the function returns the cutoff time for when the pulse amplitude falls below tpr (in dB). Default is -60.
retquadbool, optional: If True, return the quadrature (imaginary) as well as the real part of the signal. Default is False.
retenvbool, optional: If True, return the envelope of the signal. Default is False.

Returns:

yIndarray: Real part of signal. Always returned.
yQndarray: Imaginary part of signal. Only returned if retquad is True.
yenvndarray: Envelope of signal. Only returned if retenv is True.

See also

cusignal.morlet

Examples

Plot real component, imaginary component, and envelope for a 5 Hz pulse, sampled at 100 Hz for 2 seconds:

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
>>> t = cp.linspace(-1, 1, 2 * 100, endpoint=False)
>>> i, q, e = cusignal.gausspulse(t, fc=5, retquad=True, retenv=True)
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(i), cp.asnumpy(t), cp.asnumpy(q),
             cp.asnumpy(t), cp.asnumpy(e), '--')

cusignal.waveforms.waveforms.sawtooth(t, width=1.0)#

Return a periodic sawtooth or triangle waveform. The sawtooth waveform has a period 2*pi, rises from -1 to 1 on the interval 0 to width*2*pi, then drops from 1 to -1 on the interval width*2*pi to 2*pi. width must be in the interval [0, 1]. Note that this is not band-limited. It produces an infinite number of harmonics, which are aliased back and forth across the frequency spectrum. Parameters ———- t : array_like

Time.

widtharray_like, optional: Width of the rising ramp as a proportion of the total cycle. Default is 1, producing a rising ramp, while 0 produces a falling ramp. width = 0.5 produces a triangle wave. If an array, causes wave shape to change over time, and must be the same length as t.

Returns#

yndarray: Output array containing the sawtooth waveform.

Examples#

A 5 Hz waveform sampled at 500 Hz for 1 second: >>> from scipy import signal >>> import matplotlib.pyplot as plt >>> t = np.linspace(0, 1, 500) >>> plt.plot(t, signal.sawtooth(2 * np.pi * 5 * t))

cusignal.waveforms.waveforms.square(t, duty=0.5)#

Return a periodic square-wave waveform.

The square wave has a period 2*pi, has value +1 from 0 to 2*pi*duty and -1 from 2*pi*duty to 2*pi. duty must be in the interval [0,1].

Note that this is not band-limited. It produces an infinite number of harmonics, which are aliased back and forth across the frequency spectrum.

Parameters:

tarray_like: The input time array.
dutyarray_like, optional: Duty cycle. Default is 0.5 (50% duty cycle). If an array, causes wave shape to change over time, and must be the same length as t.

Returns:

yndarray: Output array containing the square waveform.

Examples

A 5 Hz waveform sampled at 500 Hz for 1 second:

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
>>> t = cp.linspace(0, 1, 500, endpoint=False)
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(cusignal.square(2 * cp.pi * 5 * t)))
>>> plt.ylim(-2, 2)

A pulse-width modulated sine wave:

>>> plt.figure()
>>> sig = cp.sin(2 * cp.pi * t)
>>> pwm = cusignal.square(2 * cp.pi * 30 * t, duty=(sig + 1)/2)
>>> plt.subplot(2, 1, 1)
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(sig))
>>> plt.subplot(2, 1, 2)
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(pwm))
>>> plt.ylim(-1.5, 1.5)

cusignal.waveforms.waveforms.unit_impulse(shape, idx=None, dtype=<class 'float'>)#

Unit impulse signal (discrete delta function) or unit basis vector.

Parameters:

shapeint or tuple of int: Number of samples in the output (1-D), or a tuple that represents the shape of the output (N-D).
idxNone or int or tuple of int or ‘mid’, optional: Index at which the value is 1. If None, defaults to the 0th element. If idx='mid', the impulse will be centered at shape // 2 in all dimensions. If an int, the impulse will be at idx in all dimensions.
dtypedata-type, optional: The desired data-type for the array, e.g., numpy.int8. Default is numpy.float64.

Returns:

yndarray: Output array containing an impulse signal.

Notes

The 1D case is also known as the Kronecker delta.

Examples

An impulse at the 0th element (\(\delta[n]\)):

>>> import cusignal
>>> import cupy as cp
>>> cusignal.unit_impulse(8)
array([ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])

Impulse offset by 2 samples (\(\delta[n-2]\)):

>>> cusignal.unit_impulse(7, 2)
array([ 0.,  0.,  1.,  0.,  0.,  0.,  0.])

2-dimensional impulse, centered:

>>> cusignal.unit_impulse((3, 3), 'mid')
array([[ 0.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  0.]])

Impulse at (2, 2), using broadcasting:

>>> cusignal.unit_impulse((4, 4), 2)
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.]])

Spectrum Analysis#

Spectral#

cusignal.spectral_analysis.spectral.coherence(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', axis=-1)#

Estimate the magnitude squared coherence estimate, Cxy, of discrete-time signals X and Y using Welch’s method.

Cxy = abs(Pxy)**2/(Pxx*Pyy), where Pxx and Pyy are power spectral density estimates of X and Y, and Pxy is the cross spectral density estimate of X and Y.

Parameters:

xarray_like: Time series of measurement values
yarray_like: Time series of measurement values
fsfloat, optional: Sampling frequency of the x and y time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
npersegint, optional: Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
noverlap: int, optional: Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None.
nfftint, optional: Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
detrendstr or function or False, optional: Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
axisint, optional: Axis along which the coherence is computed for both inputs; the default is over the last axis (i.e. axis=-1).

Returns:

fndarray: Array of sample frequencies.
Cxyndarray: Magnitude squared coherence of x and y.

See also

periodogram: Simple, optionally modified periodogram
lombscargle: Lomb-Scargle periodogram for unevenly sampled data
welch: Power spectral density by Welch’s method.
csd: Cross spectral density by Welch’s method.

Notes

An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.

References

[1]

P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967.

[2]

Stoica, Petre, and Randolph Moses, “Spectral Analysis of Signals” Prentice Hall, 2005

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt

Generate two test signals with some common features.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 20
>>> freq = 1234.0
>>> noise_power = 0.001 * fs / 2
>>> time = cp.arange(N) / fs
>>> b, a = cusignal.butter(2, 0.25, 'low')
>>> x = cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)
>>> # lfilter not implemented in cuSignal
>>> y = cusignal.lfilter(b, a, x)
>>> x += amp*cp.sin(2*cp.pi*freq*time)
>>> y += cp.random.normal(scale=0.1*cp.sqrt(noise_power), size=time.shape)

Compute and plot the coherence.

>>> f, Cxy = cusignal.coherence(x, y, fs, nperseg=1024)
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Cxy))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Coherence')
>>> plt.show()

cusignal.spectral_analysis.spectral.csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')#

Estimate the cross power spectral density, Pxy, using Welch’s method.

Parameters:

xarray_like: Time series of measurement values
yarray_like: Time series of measurement values
fsfloat, optional: Sampling frequency of the x and y time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
npersegint, optional: Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
noverlap: int, optional: Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None.
nfftint, optional: Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
detrendstr or function or False, optional: Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
return_onesidedbool, optional: If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
scaling{ ‘density’, ‘spectrum’ }, optional: Selects between computing the cross spectral density (‘density’) where Pxy has units of V**2/Hz and computing the cross spectrum (‘spectrum’) where Pxy has units of V**2, if x and y are measured in V and fs is measured in Hz. Defaults to ‘density’
axisint, optional: Axis along which the CSD is computed for both inputs; the default is over the last axis (i.e. axis=-1).
average{ ‘mean’, ‘median’ }, optional: Method to use when averaging periodograms. Defaults to ‘mean’.

Returns:

fndarray: Array of sample frequencies.
Pxyndarray: Cross spectral density or cross power spectrum of x,y.

See also

periodogram: Simple, optionally modified periodogram
lombscargle: Lomb-Scargle periodogram for unevenly sampled data
welch: Power spectral density by Welch’s method. [Equivalent to csd(x,x)]
coherence: Magnitude squared coherence by Welch’s method.

Notes

By convention, Pxy is computed with the conjugate FFT of X multiplied by the FFT of Y.

If the input series differ in length, the shorter series will be zero-padded to match.

References

[1]

[2]

Rabiner, Lawrence R., and B. Gold. “Theory and Application of Digital Signal Processing” Prentice-Hall, pp. 414-419, 1975

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt

Generate two test signals with some common features.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 20
>>> freq = 1234.0
>>> noise_power = 0.001 * fs / 2
>>> time = cp.arange(N) / fs
>>> b, a = cusignal.butter(2, 0.25, 'low')
>>> x = cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)
>>> # lfilter not currently implemented in cuSignal
>>> y = signal.lfilter(b, a, x)
>>> x += amp*cp.sin(2*cp.pi*freq*time)
>>> y += cp.random.normal(scale=0.1*cp.sqrt(noise_power), size=time.shape)

Compute and plot the magnitude of the cross spectral density.

>>> f, Pxy = cusignal.csd(x, y, fs, nperseg=1024)
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.abs(Pxy)))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('CSD [V**2/Hz]')
>>> plt.show()

cusignal.spectral_analysis.spectral.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2)#

Perform the inverse Short Time Fourier transform (iSTFT). Parameters ———- Zxx : array_like

STFT of the signal to be reconstructed. If a purely real array is passed, it will be cast to a complex data type.

fsfloat, optional: Sampling frequency of the time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window. Must match the window used to generate the STFT for faithful inversion.
npersegint, optional: Number of data points corresponding to each STFT segment. This parameter must be specified if the number of data points per segment is odd, or if the STFT was padded via nfft > nperseg. If None, the value depends on the shape of Zxx and input_onesided. If input_onesided is True, nperseg=2*(Zxx.shape[freq_axis] - 1). Otherwise, nperseg=Zxx.shape[freq_axis]. Defaults to None.
noverlapint, optional: Number of points to overlap between segments. If None, half of the segment length. Defaults to None. When specified, the COLA constraint must be met (see Notes below), and should match the parameter used to generate the STFT. Defaults to None.
nfftint, optional: Number of FFT points corresponding to each STFT segment. This parameter must be specified if the STFT was padded via nfft > nperseg. If None, the default values are the same as for nperseg, detailed above, with one exception: if input_onesided is True and nperseg==2*Zxx.shape[freq_axis] - 1, nfft also takes on that value. This case allows the proper inversion of an odd-length unpadded STFT using nfft=None. Defaults to None.
input_onesidedbool, optional: If True, interpret the input array as one-sided FFTs, such as is returned by stft with return_onesided=True and numpy.fft.rfft. If False, interpret the input as a a two-sided FFT. Defaults to True.
boundarybool, optional: Specifies whether the input signal was extended at its boundaries by supplying a non-None boundary argument to stft. Defaults to True.
time_axisint, optional: Where the time segments of the STFT is located; the default is the last axis (i.e. axis=-1).
freq_axisint, optional: Where the frequency axis of the STFT is located; the default is the penultimate axis (i.e. axis=-2).

Returns#

tndarray: Array of output data times.
xndarray: iSTFT of Zxx.

Examples#

>>> from scipy import signal
>>> import matplotlib.pyplot as plt
Generate a test signal, a 2 Vrms sine wave at 50Hz corrupted by
0.001 V**2/Hz of white noise sampled at 1024 Hz.
>>> fs = 1024
>>> N = 10*fs
>>> nperseg = 512
>>> amp = 2 * np.sqrt(2)
>>> noise_power = 0.001 * fs / 2
>>> time = cp.arange(N) / float(fs)
>>> carrier = amp * cp.sin(2*cp.pi*50*time)
>>> noise = cp.random.normal(scale=cp.sqrt(noise_power),
...                          size=time.shape)
>>> x = carrier + noise
Compute the STFT, and plot its magnitude
>>> f, t, Zxx = cusignal.stft(x, fs=fs, nperseg=nperseg)
>>> f = cp.asnumpy(f)
>>> t = cp.asnumpy(t)
>>> Zxx = cp.asnumpy(Zxx)
>>> plt.figure()
>>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud')
>>> plt.ylim([f[1], f[-1]])
>>> plt.title('STFT Magnitude')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.yscale('log')
>>> plt.show()
Zero the components that are 10% or less of the carrier magnitude,
then convert back to a time series via inverse STFT
>>> Zxx = cp.where(cp.abs(Zxx) >= amp/10, Zxx, 0)
>>> _, xrec = cusignal.istft(Zxx, fs)
>>> xrec = cp.asnumpy(xrec)
>>> x = cp.asnumpy(x)
>>> time = cp.asnumpy(time)
>>> carrier = cp.asnumpy(carrier)
Compare the cleaned signal with the original and true carrier signals.
>>> plt.figure()
>>> plt.plot(time, x, time, xrec, time, carrier)
>>> plt.xlim([2, 2.1])*+
>>> plt.xlabel('Time [sec]')
>>> plt.ylabel('Signal')
>>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
>>> plt.show()
Note that the cleaned signal does not start as abruptly as the original,
since some of the coefficients of the transient were also removed:
>>> plt.figure()
>>> plt.plot(time, x, time, xrec, time, carrier)
>>> plt.xlim([0, 0.1])
>>> plt.xlabel('Time [sec]')
>>> plt.ylabel('Signal')
>>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
>>> plt.show()

cusignal.spectral_analysis.spectral.lombscargle(x, y, freqs)#

Computes the Lomb-Scargle periodogram. The Lomb-Scargle periodogram was developed by Lomb [1] and further extended by Scargle [2] to find, and test the significance of weak periodic signals with uneven temporal sampling. When normalize is False (default) the computed periodogram is unnormalized, it takes the value (A**2) * N/4 for a harmonic signal with amplitude A for sufficiently large N. When normalize is True the computed periodogram is normalized by the residuals of the data around a constant reference model (at zero). Input arrays should be one-dimensional and will be cast to float64.

Parameters:

xarray_like: Sample times.
yarray_like: Measurement values.
freqsarray_like: Angular frequencies for output periodogram.
precenterbool, optional: Pre-center amplitudes by subtracting the mean.
normalizebool, optional: Compute normalized periodogram.

Returns:

pgramarray_like: Lomb-Scargle periodogram.

Raises:

ValueError: If the input arrays x and y do not have the same shape.

See also

istft: Inverse Short Time Fourier Transform
check_COLA: Check whether the Constant OverLap Add (COLA) constraint is met
welch: Power spectral density by Welch’s method
spectrogram: Spectrogram by Welch’s method
csd: Cross spectral density by Welch’s method

Notes

This subroutine calculates the periodogram using a slightly modified algorithm due to Townsend [3] which allows the periodogram to be calculated using only a single pass through the input arrays for each frequency. The algorithm running time scales roughly as O(x * freqs) or O(N^2) for a large number of samples and frequencies.

References

[1]

N.R. Lomb “Least-squares frequency analysis of unequally spaced data”, Astrophysics and Space Science, vol 39, pp. 447-462, 1976

[2]

J.D. Scargle “Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data”, The Astrophysical Journal, vol 263, pp. 835-853, 1982

[3]

R.H.D. Townsend, “Fast calculation of the Lomb-Scargle periodogram using graphics processing units.”, The Astrophysical Journal Supplement Series, vol 191, pp. 247-253, 2010

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
First define some input parameters for the signal:
>>> A = 2.
>>> w = 1.
>>> phi = 0.5 * cp.pi
>>> nin = 1000
>>> nout = 100000
>>> frac_points = 0.9 # Fraction of points to select
Randomly select a fraction of an array with timesteps:
>>> r = cp.random.rand(nin)
>>> x = cp.linspace(0.01, 10*cp.pi, nin)
>>> x = x[r >= frac_points]
Plot a sine wave for the selected times:
>>> y = A * cp.sin(w*x+phi)
Define the array of frequencies for which to compute the periodogram:
>>> f = cp.linspace(0.01, 10, nout)
Calculate Lomb-Scargle periodogram:
>>> pgram = cusignal.lombscargle(x, y, f, normalize=True)
Now make a plot of the input data:
>>> plt.subplot(2, 1, 1)
>>> plt.plot(cp.asnumpy(x), cp.asnumpy(y), 'b+')
Then plot the normalized periodogram:
>>> plt.subplot(2, 1, 2)
>>> plt.plot(cp.asnumpy(f), cp.asnumpy(pgram))
>>> plt.show()

cusignal.spectral_analysis.spectral.periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1)#

Estimate power spectral density using a periodogram.

Parameters:

xarray_like: Time series of measurement values
fsfloat, optional: Sampling frequency of the x time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to ‘boxcar’.
nfftint, optional: Length of the FFT used. If None the length of x will be used.
detrendstr or function or False, optional: Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
return_onesidedbool, optional: If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
scaling{ ‘density’, ‘spectrum’ }, optional: Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’
axisint, optional: Axis along which the periodogram is computed; the default is over the last axis (i.e. axis=-1).

Returns:

fndarray: Array of sample frequencies.
Pxxndarray: Power spectral density or power spectrum of x.

See also

welch: Estimate power spectral density using Welch’s method
lombscargle: Lomb-Scargle periodogram for unevenly sampled data

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
>>> cp.random.seed(1234)

Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*cp.sqrt(2)
>>> freq = 1234.0
>>> noise_power = 0.001 * fs / 2
>>> time = cp.arange(N) / fs
>>> x = amp*cp.sin(2*cp.pi*freq*time)
>>> x += cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)

Compute and plot the power spectral density.

>>> f, Pxx_den = cusignal.periodogram(x, fs)
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Pxx_den))
>>> plt.ylim([1e-7, 1e2])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.show()

If we average the last half of the spectral density, to exclude the peak, we can recover the noise power on the signal.

>>> cp.mean(Pxx_den[25000:])
0.00099728892368242854

Now compute and plot the power spectrum.

>>> f, Pxx_spec = cusignal.periodogram(x, fs, 'flattop',             scaling='spectrum')
>>> plt.figure()
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.sqrt(Pxx_spec)))
>>> plt.ylim([1e-4, 1e1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Linear spectrum [V RMS]')
>>> plt.show()

The peak height in the power spectrum is an estimate of the RMS amplitude.

>>> cp.sqrt(Pxx_spec.max())
2.0077340678640727

cusignal.spectral_analysis.spectral.spectrogram(x, fs=1.0, window=('tukey', 0.25), nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, mode='psd')#

Compute a spectrogram with consecutive Fourier transforms.

Spectrograms can be used as a way of visualizing the change of a nonstationary signal’s frequency content over time.

Parameters:

xarray_like: Time series of measurement values
fsfloat, optional: Sampling frequency of the x time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Tukey window with shape parameter of 0.25.
npersegint, optional: Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
noverlapint, optional: Number of points to overlap between segments. If None, noverlap = nperseg // 8. Defaults to None.
nfftint, optional: Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
detrendstr or function or False, optional: Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
return_onesidedbool, optional: If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
scaling{ ‘density’, ‘spectrum’ }, optional: Selects between computing the power spectral density (‘density’) where Sxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Sxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’.
axisint, optional: Axis along which the spectrogram is computed; the default is over the last axis (i.e. axis=-1).
modestr, optional: Defines what kind of return values are expected. Options are [‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’]. ‘complex’ is equivalent to the output of stft with no padding or boundary extension. ‘magnitude’ returns the absolute magnitude of the STFT. ‘angle’ and ‘phase’ return the complex angle of the STFT, with and without unwrapping, respectively.

Returns:

fndarray: Array of sample frequencies.
tndarray: Array of segment times.
Sxxndarray: Spectrogram of x. By default, the last axis of Sxx corresponds to the segment times.

See also

periodogram: Simple, optionally modified periodogram
lombscargle: Lomb-Scargle periodogram for unevenly sampled data
welch: Power spectral density by Welch’s method.
csd: Cross spectral density by Welch’s method.

Notes

An appropriate amount of overlap will depend on the choice of window and on your requirements. In contrast to welch’s method, where the entire data stream is averaged over, one may wish to use a smaller overlap (or perhaps none at all) when computing a spectrogram, to maintain some statistical independence between individual segments. It is for this reason that the default window is a Tukey window with 1/8th of a window’s length overlap at each end.

New in version 0.16.0.

References

[1]

Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”, Prentice Hall, 1999.

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt

Generate a test signal, a 2 Vrms sine wave whose frequency is slowly modulated around 3kHz, corrupted by white noise of exponentially decreasing magnitude sampled at 10 kHz.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2 * cp.sqrt(2)
>>> noise_power = 0.01 * fs / 2
>>> time = cp.arange(N) / float(fs)
>>> mod = 500*cp.cos(2*cp.pi*0.25*time)
>>> carrier = amp * cp.sin(2*cp.pi*3e3*time + mod)
>>> noise = cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)
>>> noise *= cp.exp(-time/5)
>>> x = carrier + noise

Compute and plot the spectrogram.

>>> f, t, Sxx = cusignal.spectrogram(x, fs)
>>> plt.pcolormesh(cp.asnumpy(t), cp.asnumpy(f), cp.asnumpy(Sxx))
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()

Note, if using output that is not one sided, then use the following:

>>> f, t, Sxx = cusignal.spectrogram(x, fs, return_onesided=False)
>>> plt.pcolormesh(cp.asnumpy(t), cp.fft.fftshift(f),         cp.fft.fftshift(Sxx, axes=0))
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()

cusignal.spectral_analysis.spectral.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend=False, return_onesided=True, boundary='zeros', padded=True, axis=-1)#

Compute the Short Time Fourier Transform (STFT).

STFTs can be used as a way of quantifying the change of a nonstationary signal’s frequency and phase content over time.

Parameters:

xarray_like: Time series of measurement values
fsfloat, optional: Sampling frequency of the x time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
npersegint, optional: Length of each segment. Defaults to 256.
noverlapint, optional: Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None. When specified, the COLA constraint must be met (see Notes below).
nfftint, optional: Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
detrendstr or function or False, optional: Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to False.
return_onesidedbool, optional: If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
boundarystr or None, optional: Specifies whether the input signal is extended at both ends, and how to generate the new values, in order to center the first windowed segment on the first input point. This has the benefit of enabling reconstruction of the first input point when the employed window function starts at zero. Valid options are ['even', 'odd', 'constant', 'zeros', None]. Defaults to ‘zeros’, for zero padding extension. I.e. [1, 2, 3, 4] is extended to [0, 1, 2, 3, 4, 0] for nperseg=3.
paddedbool, optional: Specifies whether the input signal is zero-padded at the end to make the signal fit exactly into an integer number of window segments, so that all of the signal is included in the output. Defaults to True. Padding occurs after boundary extension, if boundary is not None, and padded is True, as is the default.
axisint, optional: Axis along which the STFT is computed; the default is over the last axis (i.e. axis=-1).

Returns:

fndarray: Array of sample frequencies.
tndarray: Array of segment times.
Zxxndarray: STFT of x. By default, the last axis of Zxx corresponds to the segment times.

See also

welch: Power spectral density by Welch’s method.
spectrogram: Spectrogram by Welch’s method.
csd: Cross spectral density by Welch’s method.
lombscargle: Lomb-Scargle periodogram for unevenly sampled data

Notes

In order to enable inversion of an STFT via the inverse STFT in istft, the signal windowing must obey the constraint of “Nonzero OverLap Add” (NOLA), and the input signal must have complete windowing coverage (i.e. (x.shape[axis] - nperseg) % (nperseg-noverlap) == 0). The padded argument may be used to accomplish this.

Given a time-domain signal \(x[n]\), a window \(w[n]\), and a hop size \(H\) = nperseg - noverlap, the windowed frame at time index \(t\) is given by

\[x_{t}[n]=x[n]w[n-tH]\]

The overlap-add (OLA) reconstruction equation is given by

\[x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]}\]

The NOLA constraint ensures that every normalization term that appears in the denomimator of the OLA reconstruction equation is nonzero. Whether a choice of window, nperseg, and noverlap satisfy this constraint can be tested with check_NOLA.

References

[1]

Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”, Prentice Hall, 1999.

[2]

Daniel W. Griffin, Jae S. Lim “Signal Estimation from Modified Short-Time Fourier Transform”, IEEE 1984, 10.1109/TASSP.1984.1164317

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt

Generate a test signal, a 2 Vrms sine wave whose frequency is slowly modulated around 3kHz, corrupted by white noise of exponentially decreasing magnitude sampled at 10 kHz.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2 * cp.sqrt(2)
>>> noise_power = 0.01 * fs / 2
>>> time = cp.arange(N) / float(fs)
>>> mod = 500*cp.cos(2*cp.pi*0.25*time)
>>> carrier = amp * cp.sin(2*cp.pi*3e3*time + mod)
>>> noise = cp.random.normal(scale=cp.sqrt(noise_power),
...                          size=time.shape)
>>> noise *= cp.exp(-time/5)
>>> x = carrier + noise

Compute and plot the STFT’s magnitude.

>>> f, t, Zxx = cusignal.stft(x, fs, nperseg=1000)
>>> plt.pcolormesh(cp.asnumpy(t), cp.asnumpy(f), cp.asnumpy(cp.abs(Zxx)), \
    vmin=0, vmax=amp)
>>> plt.title('STFT Magnitude')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()

cusignal.spectral_analysis.spectral.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”.

Parameters:

events1D array_like: An array of time points containing the timing of the events.
periodfloat 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.

Returns:

strengthfloat 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.
phasefloat 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.

References

[1]

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`.

[2]

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`.

[3]

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`.

cusignal.spectral_analysis.spectral.welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')#

Estimate power spectral density using Welch’s method.

Welch’s method [1] computes an estimate of the power spectral density by dividing the data into overlapping segments, computing a modified periodogram for each segment and averaging the periodograms.

Parameters:

xarray_like: Time series of measurement values
fsfloat, optional: Sampling frequency of the x time series. Defaults to 1.0.
windowstr or tuple or array_like, optional: Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
npersegint, optional: Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
noverlapint, optional: Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None.
nfftint, optional: Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
detrendstr or function or False, optional: Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
return_onesidedbool, optional: If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
scaling{ ‘density’, ‘spectrum’ }, optional: Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’
axisint, optional: Axis along which the periodogram is computed; the default is over the last axis (i.e. axis=-1).
average{ ‘mean’, ‘median’ }, optional: Method to use when averaging periodograms. Defaults to ‘mean’.

Returns:

fndarray: Array of sample frequencies.
Pxxndarray: Power spectral density or power spectrum of x.

See also

periodogram: Simple, optionally modified periodogram
lombscargle: Lomb-Scargle periodogram for unevenly sampled data

Notes

If noverlap is 0, this method is equivalent to Bartlett’s method [2].

New in version 0.12.0.

References

[1]

[2]

M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika, vol. 37, pp. 1-16, 1950.

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
>>> cp.random.seed(1234)

Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*cp.sqrt(2)
>>> freq = 1234.0
>>> noise_power = 0.001 * fs / 2
>>> time = cp.arange(N) / fs
>>> x = amp*cp.sin(2*cp.pi*freq*time)
>>> x += cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)

Compute and plot the power spectral density.

>>> f, Pxx_den = cusignal.welch(x, fs, nperseg=1024)
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Pxx_den))
>>> plt.ylim([0.5e-3, 1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.show()

If we average the last half of the spectral density, to exclude the peak, we can recover the noise power on the signal.

>>> cp.mean(Pxx_den[256:])
0.0009924865443739191

Now compute and plot the power spectrum.

>>> f, Pxx_spec = cusignal.welch(x, fs, 'flattop', 1024, \
    scaling='spectrum')
>>> plt.figure()
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.sqrt(Pxx_spec)))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Linear spectrum [V RMS]')
>>> plt.show()

The peak height in the power spectrum is an estimate of the RMS amplitude.

>>> cp.sqrt(Pxx_spec.max())
2.0077340678640727

If we now introduce a discontinuity in the signal, by increasing the amplitude of a small portion of the signal by 50, we can see the corruption of the mean average power spectral density, but using a median average better estimates the normal behaviour.

>>> x[int(N//2):int(N//2)+10] *= 50.
>>> f, Pxx_den = cusignal.welch(x, fs, nperseg=1024)
>>> f_med, Pxx_den_med = cusignal.welch(x, fs, nperseg=1024,
                                      average='median')
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Pxx_den), label='mean')
>>> plt.semilogy(cp.asnumpy(f_med), cp.asnumpy(Pxx_den_med), \
    label='median')
>>> plt.ylim([0.5e-3, 1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.legend()
>>> plt.show()

Acoustics#

cusignal.acoustics.cepstrum.complex_cepstrum(x, n=None, axis=-1)#

Calculates the complex cepstrum of a real valued input sequence x where the cepstrum is defined as the inverse Fourier transform of the log magnitude DFT (spectrum) of a signal. It’s primarily used for source/speaker separation in speech signal processing.

The input is altered to have zero-phase at pi radians (180 degrees) Parameters ———- x : ndarray

Input sequence, if x is a matrix, return cepstrum in direction of axis

nint: Size of Fourier Transform; If none, will use length of input array
axis: int: Direction for cepstrum calculation

Returns#

cepsndarray: Complex cepstrum result

cusignal.acoustics.cepstrum.inverse_complex_cepstrum(ceps, ndelay)#

Compute the inverse complex cepstrum of a real sequence. ceps : ndarray

Real sequence to compute inverse complex cepstrum of.

ndelay: int: The amount of samples of circular delay added to x.

Returns#

xndarray: The inverse complex cepstrum of the real sequence ceps.

The inverse complex cepstrum is given by .. math:: x[n] = F^{-1}left{exp(F(c[n]))right} where \(c_[n]\) is the input signal and \(F\) and :math:`F_{-1} are respectively the forward and backward Fourier transform.

cusignal.acoustics.cepstrum.minimum_phase(x, n=None)#

Compute the minimum phase reconstruction of a real sequence. x : ndarray

Real sequence to compute the minimum phase reconstruction of.

n{None, int}, optional: Length of the Fourier transform.

Compute the minimum phase reconstruction of a real sequence using the real cepstrum. Returns ——- m : ndarray

The minimum phase reconstruction of the real sequence x.

cusignal.acoustics.cepstrum.real_cepstrum(x, n=None, axis=-1)#

Calculates the real cepstrum of an input sequence x where the cepstrum is defined as the inverse Fourier transform of the log magnitude DFT (spectrum) of a signal. It’s primarily used for source/speaker separation in speech signal processing

Parameters:

xndarray: Input sequence, if x is a matrix, return cepstrum in direction of axis
nint: Size of Fourier Transform; If none, will use length of input array
axis: int: Direction for cepstrum calculation

Returns:

cepsndarray: Complex cepstrum result

Wavelets#

cusignal.wavelets.wavelets.cwt(data, wavelet, widths)#

Continuous wavelet transform.

Performs a continuous wavelet transform on data, using the wavelet function. A CWT performs a convolution with data using the wavelet function, which is characterized by a width parameter and length parameter.

Parameters:

data(N,) ndarray: data on which to perform the transform.
waveletfunction: Wavelet function, which should take 2 arguments. The first argument is the number of points that the returned vector will have (len(wavelet(length,width)) == length). The second is a width parameter, defining the size of the wavelet (e.g. standard deviation of a gaussian). See ricker, which satisfies these requirements.
widths(M,) sequence: Widths to use for transform.

Returns:

cwt: (M, N) ndarray: Will have shape of (len(widths), len(data)).

Notes

length = min(10 * width[ii], len(data))
cwt[ii,:] = cusignal.convolve(data, wavelet(length,
                            width[ii]), mode='same')

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
>>> t = cp.linspace(-1, 1, 200, endpoint=False)
>>> sig  = cp.cos(2 * cp.pi * 7 * t) + cusignal.gausspulse(t - 0.4, fc=2)
>>> widths = cp.arange(1, 31)
>>> cwtmatr = cusignal.cwt(sig, cusignal.ricker, widths)
>>> plt.imshow(abs(cp.asnumpy(cwtmatr)), extent=[-1, 1, 31, 1],
               cmap='PRGn', aspect='auto', vmax=abs(cwtmatr).max(),
               vmin=-abs(cwtmatr).max())
>>> plt.show()

cusignal.wavelets.wavelets.morlet(M, w=5.0, s=1.0, complete=True)#

Complex Morlet wavelet.

Parameters:

Mint: Length of the wavelet.
wfloat, optional: Omega0. Default is 5
sfloat, optional: Scaling factor, windowed from -s*2*pi to +s*2*pi. Default is 1.
completebool, optional: Whether to use the complete or the standard version.

Returns:

morlet(M,) ndarray

See also

cusignal.gausspulse

Notes

The standard version:

pi**-0.25 * exp(1j*w*x) * exp(-0.5*(x**2))

This commonly used wavelet is often referred to simply as the Morlet wavelet. Note that this simplified version can cause admissibility problems at low values of w.

The complete version:

pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2))

This version has a correction term to improve admissibility. For w greater than 5, the correction term is negligible.

Note that the energy of the return wavelet is not normalised according to s.

The fundamental frequency of this wavelet in Hz is given by f = 2*s*w*r / M where r is the sampling rate.

Note: This function was created before cwt and is not compatible with it.

cusignal.wavelets.wavelets.morlet2(M, s, w=5)#

Complex Morlet wavelet, designed to work with cwt. Returns the complete version of morlet wavelet, normalised according to s:

exp(1j*w*x/s) * exp(-0.5*(x/s)**2) * pi**(-0.25) * sqrt(1/s)

Parameters#

Mint: Length of the wavelet.
sfloat: Width parameter of the wavelet.
wfloat, optional: Omega0. Default is 5

Returns#

morlet : (M,) ndarray See Also ——– morlet : Implementation of Morlet wavelet, incompatible with cwt Notes —– .. versionadded:: 1.4.0 This function was designed to work with cwt. Because morlet2 returns an array of complex numbers, the dtype argument of cwt should be set to complex128 for best results. Note the difference in implementation with morlet. The fundamental frequency of this wavelet in Hz is given by:

f = w*fs / (2*s*np.pi)

where fs is the sampling rate and s is the wavelet width parameter. Similarly we can get the wavelet width parameter at f:

s = w*fs / (2*f*np.pi)

Examples#

>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> M = 100
>>> s = 4.0
>>> w = 2.0
>>> wavelet = signal.morlet2(M, s, w)
>>> plt.plot(abs(wavelet))
>>> plt.show()
This example shows basic use of `morlet2` with `cwt` in time-frequency
analysis:
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t, dt = np.linspace(0, 1, 200, retstep=True)
>>> fs = 1/dt
>>> w = 6.
>>> sig = np.cos(2*np.pi*(50 + 10*t)*t) + np.sin(40*np.pi*t)
>>> freq = np.linspace(1, fs/2, 100)
>>> widths = w*fs / (2*freq*np.pi)
>>> cwtm = signal.cwt(sig, signal.morlet2, widths, w=w)
>>> plt.pcolormesh(t, freq, np.abs(cwtm),
    cmap='viridis', shading='gouraud')
>>> plt.show()

cusignal.wavelets.wavelets.qmf(hk)#

Return high-pass qmf filter from low-pass

Parameters:

hkarray_like: Coefficients of high-pass filter.

cusignal.wavelets.wavelets.ricker(points, a)#

Return a Ricker wavelet, also known as the “Mexican hat wavelet”.

It models the function:

A (1 - x^2/a^2) exp(-x^2/2 a^2),

where A = 2/sqrt(3a)pi^1/4.

Parameters:

pointsint: Number of points in vector. Will be centered around 0.
ascalar: Width parameter of the wavelet.

Returns:

vector(N,) ndarray: Array of length points in shape of ricker curve.

Examples

>>> import cusignal
>>> import cupy as cp
>>> import matplotlib.pyplot as plt

>>> points = 100
>>> a = 4.0
>>> vec2 = cusignal.ricker(points, a)
>>> print(len(vec2))
100
>>> plt.plot(cp.asnumpy(vec2))
>>> plt.show()

B-splines#

cusignal.bsplines.bsplines.cubic(x)#

A cubic B-spline.

This is a special case of bspline, and equivalent to bspline(x, 3).

cusignal.bsplines.bsplines.gauss_spline(x, n)#

Gaussian approximation to B-spline basis function of order n.

Parameters:

nint: The order of the spline. Must be nonnegative, i.e. n >= 0

References

[1]

Bouma H., Vilanova A., Bescos J.O., ter Haar Romeny B.M., Gerritsen F.A. (2007) Fast and Accurate Gaussian Derivatives Based on B-Splines. In: Sgallari F., Murli A., Paragios N. (eds) Scale Space and Variational Methods in Computer Vision. SSVM 2007. Lecture Notes in Computer Science, vol 4485. Springer, Berlin, Heidelberg

cusignal.bsplines.bsplines.quadratic(x)#

A quadratic B-spline.

This is a special case of bspline, and equivalent to bspline(x, 2).

Utilities#

Array Tools#

cusignal.utils.arraytools.from_pycuda(pycuda_arr, device=0)#

Read in gpuarray from PyCUDA and output CuPy array

Parameters:

pycuda_arrPyCUDA gpuarray
deviceint: GPU Device ID

Returns:

cupy_arrCuPy ndarray

cusignal.utils.arraytools.get_pinned_array(data)#

Return populated pinned memory.

Parameters:

datacupy.ndarray or numpy.ndarray: The array to be copied to shared buffer
strides: int or None
order: char

cusignal.utils.arraytools.get_pinned_mem(shape, dtype)#

Create a pinned memory allocation.

Parameters:

shapeint or tuple of ints: Output shape.
dtypedata-type: Output data type.

Returns:

retndarray: Pinned memory numpy array.

cusignal.utils.arraytools.get_shared_array(data, strides=None, order='C', stream=0, portable=False, wc=True)#

Return populated shared memory between GPU and CPU.

Parameters:

datacupy.ndarray or numpy.ndarray: The array to be copied to shared buffer
strides: int or None
order: char
streamint: Stream number (0 for default)
portablebool
wcbool

cusignal.utils.arraytools.get_shared_mem(shape, dtype=<class 'numpy.float32'>, strides=None, order='C', stream=0, portable=False, wc=True)#

Return shared memory between GPU and CPU. Similar to numpy.zeros

Parameters:

shapendarray.shape: Size of shared memory allocation
dtypecupy.dtype or numpy.dtype: Data type of allocation
strides: int or None
order: char
streamint: Stream number (0 for default)
portablebool
wcbool

FFTPack Helper#

cusignal.utils.fftpack_helper.next_fast_len(target)#

Find the next fast size of input data to fft, for zero-padding, etc.

SciPy’s FFTPACK has efficient functions for radix {2, 3, 4, 5}, so this returns the next composite of the prime factors 2, 3, and 5 which is greater than or equal to target. (These are also known as 5-smooth numbers, regular numbers, or Hamming numbers.)

Parameters:

targetint: Length to start searching from. Must be a positive integer.

Returns:

outint: The first 5-smooth number greater than or equal to target.

Notes

New in version 0.18.0.

Examples

On a particular machine, an FFT of prime length takes 133 ms:

>>> from scipy import fftpack
>>> min_len = 10007  # prime length is worst case for speed
>>> a = np.random.randn(min_len)
>>> b = fftpack.fft(a)

Zero-padding to the next 5-smooth length reduces computation time to 211 us, a speedup of 630 times:

>>> fftpack.helper.next_fast_len(min_len)
10125
>>> b = fftpack.fft(a, 10125)

Rounding up to the next power of 2 is not optimal, taking 367 us to compute, 1.7 times as long as the 5-smooth size:

>>> b = fftpack.fft(a, 16384)

IO#

Reader#

cusignal.io.reader.read_bin(file, buffer=None, dtype=<class 'numpy.uint8'>, num_samples=None, offset=0)#

Reads binary file into GPU memory. Can be used as a building blocks for custom unpack/pack data readers/writers.

Parameters:

filestr: A string of filename to be read to GPU.
bufferndarray, optional: Pinned memory buffer to use when copying data from GPU.
dtypedata-type, optional: Any object that can be interpreted as a numpy data type.
num_samplesint, optional: Number of samples to be loaded to GPU. If set to 0, read in all samples.
offsetint, optional: In the file, array data starts at this offset. Since offset is measured in bytes, it should normally be a multiple of the byte-size of dtype.
Returns
——-
outndarray: An 1-dimensional array containing binary data.

cusignal.io.reader.read_sigmf(data_file, meta_file=None, buffer=None, num_samples=None, offset=0)#

Read and unpack binary file, with SigMF spec, to GPU memory.

Parameters:

data_filestr: File contain sigmf data.
meta_filestr, optional: File contain sigmf meta.
bufferndarray, optional: Pinned memory buffer to use when copying data from GPU.
num_samplesint, optional: Number of samples to be loaded to GPU. If set to 0, read in all samples.
offsetint, optional: May be specified as a non-negative integer offset. It is the number of samples before loading ‘num_samples’. ‘offset’ must be a multiple of ALLOCATIONGRANULARITY which is equal to PAGESIZE on Unix systems.

Returns:

outndarray: An 1-dimensional array containing unpacked binary data.

cusignal.io.reader.unpack_bin(binary, dtype, endianness='L')#

Unpack binary file. If endianness is big-endian, it my be converted to little endian for NVIDIA GPU compatibility.

Parameters:

binaryndarray: The binary array to be unpack.
dtypedata-type, optional: Any object that can be interpreted as a numpy data type.
endianness{‘L’, ‘B’}, optional: Data set byte order

Returns:

outndarray: An 1-dimensional array containing unpacked binary data.

Writer#

cusignal.io.writer.pack_bin(in1)#

Pack binary arrary. Data will be packed with little endian for NVIDIA GPU compatibility.

Parameters:

in1ndarray: The ndarray to be pack at binary.

Returns:

outndarray: An 1-dimensional array containing packed binary data.

cusignal.io.writer.write_bin(file, binary, buffer=None, append=True)#

Writes binary array to file.

Parameters:

filestr: A string of filename to store output.
binaryndarray: Binary array to be written to file.
bufferndarray, optional: Pinned memory buffer to use when copying data from GPU.
appendbool, optional: Append to file if created.

Returns:

outndarray: An 1-dimensional array containing binary data.

cusignal.io.writer.write_sigmf(data_file, data, buffer=None, append=True)#

Pack and write binary array to file, with SigMF spec.

Parameters:

filestr: A string of filename to be read/unpacked to GPU.
binaryndarray: Binary array to be written to file.
bufferndarray, optional: Pinned memory buffer to use when copying data from GPU.
appendbool, optional: Append to file if created.

Returns:

Radar/Phased Array#

Radar Tools#

cusignal.radartools.radartools.ambgfun(x, fs, prf, y=None, cut='2d', cutValue=0)#

Calculates the normalized ambiguity function for the vector x

Parameters:

xndarray: Input pulse waveform
fs: int, float: Sampling rate in Hz
prf: int, float: Pulse repetition frequency in Hz
yndarray: Second input pulse waveform. If not given, y = x
cutstring: Direction of one-dimensional cut through ambiguity function
cutValueint, float: Time delay or doppler shift at which one-dimensional cut through ambiguity function is taken

Returns:

amfunndarray: Normalized magnitude of the ambiguity function

cusignal.radartools.radartools.ca_cfar(array, guard_cells, reference_cells, pfa=0.001)#

Computes the cell-averaged constant false alarm rate (CA CFAR) detector threshold and returns for a given array.

Parameters:

arrayndarray: Array containing data to be processed.
guard_cells_xint: One-sided guard cell count in the first dimension.
guard_cells_yint: One-sided guard cell count in the second dimension.
reference_cells_xint: one-sided reference cell count in the first dimension.
reference_cells_yint: one-sided referenc cell count in the second dimension.
pfafloat: Probability of false alarm.

Returns:

thresholdndarray: CFAR threshold
returnndarray: CFAR detections

cusignal.radartools.radartools.cfar_alpha(pfa, N)#

Computes the value of alpha corresponding to a given probability of false alarm and number of reference cells N.

Parameters:

pfafloat: Probability of false alarm.
Nint: Number of reference cells.

Returns:

alphafloat: Alpha value.

cusignal.radartools.radartools.pulse_compression(x, template, normalize=False, window=None, nfft=None)#

Pulse Compression is used to increase the range resolution and SNR by performing matched filtering of the transmitted pulse (template) with the received signal (x)

Parameters:

xndarray: Received signal, assume 2D array with [num_pulses, sample_per_pulse]
templatendarray: Transmitted signal, assume 1D array
normalizebool: Normalize transmitted signal
windowarray_like, callable, string, float, or tuple, optional: Specifies the window applied to the signal in the Fourier domain.
nfftint, size of FFT for pulse compression. Default is number of: samples per pulse

Returns:

compressedIQndarray: Pulse compressed output

cusignal.radartools.radartools.pulse_doppler(x, window=None, nfft=None)#

Pulse doppler processing yields a range/doppler data matrix that represents moving target data that’s separated from clutter. An estimation of the doppler shift can also be obtained from pulse doppler processing. FFT taken across slow-time (pulse) dimension.

Parameters:

xndarray: Received signal, assume 2D array with [num_pulses, sample_per_pulse]
windowarray_like, callable, string, float, or tuple, optional: Specifies the window applied to the signal in the Fourier domain.
nfftint, size of FFT for pulse compression. Default is number of: samples per pulse

Returns:

pd_dataMatrixndarray: Pulse-doppler output (range/doppler matrix)

cusignal.radartools.beamformers.mvdr(x, sv, calc_cov=True)#

Minimum variance distortionless response (MVDR) beamformer weights

Parameters:

xndarray: Received signal or input covariance matrix, assume 2D array with size [num_sensors, num_samples]
sv: ndarray: Steering vector, assume 1D array with size [num_sensors, 1]
calc_covbool: Determine whether to calculate covariance matrix. Simply put, calc_cov defines whether x input is made of sensor/observation data or is a precalculated covariance matrix
Note: Unlike MATLAB where input matrix x is of size MxN where N represents
the number of array elements, we assume row-major formatted data where each
row is assumed to be complex-valued data from a given sensor (i.e. NxM)

cuSignal API Reference#

Convolution#

Convolve#

Correlate#

See Also#

Notes#

Modulation/Demodulation#

Demodulation#

Estimation#

Kalman Filter#

Filtering#

Resample#

FIR Filters#

Parameters#

Returns#

Returns#

See also#

References#

Parameters#

Returns#

Returns#

References#

Channelizer#

Parameters#

Returns#

Notes#

Filter Design#

Resample#

Returns#

See also#

References#

Parameters#

Returns#

Returns#

References#

Peak Finding#

Peak Finding#

Window Functions#

Windows#

Returns#

See Also#

Examples#

Waveforms#

Waveform Generation#

Returns#

Examples#

Spectrum Analysis#

Spectral#

Returns#

See Also#

Examples#

Acoustics#

Returns#

Returns#

Wavelets#

Wavelets#

Parameters#

Returns#

Examples#

B-splines#

B-splines#

Utilities#

Array Tools#

FFTPack Helper#

IO#

Reader#

Writer#

Radar/Phased Array#

Radar Tools#