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.

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.

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.

See Also#

correlate : Compute the N-dimensional cross-correlation. Returns ——- lags : array

Returns an array containing cross-correlation lag/displacement indices. Indices can be indexed with the np.argmax of the correlation to return the lag/displacement.

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.

See also#

firls firwin minimum_phase remez Notes —– From the given set of frequencies and gains, the desired response is constructed in the frequency domain. The inverse FFT is applied to the desired response to create the associated convolution kernel, and the first numtaps coefficients of this kernel, scaled by window, are returned. The FIR filter will have linear phase. The type of filter is determined by the value of ‘numtaps` and antisymmetric flag. There are four possible combinations:

  • odd numtaps, antisymmetric is False, type I filter is produced

  • even numtaps, antisymmetric is False, type II filter is produced

  • odd numtaps, antisymmetric is True, type III filter is produced

  • even numtaps, antisymmetric is True, type IV filter is produced

Magnitude response of all but type I filters are subjects to following constraints:

  • type II – zero at the Nyquist frequency

  • type III – zero at zero and Nyquist frequencies

  • type IV – zero at zero frequency

New in version 0.9.0.

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.

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.

See also#

firls firwin minimum_phase remez Notes —– From the given set of frequencies and gains, the desired response is constructed in the frequency domain. The inverse FFT is applied to the desired response to create the associated convolution kernel, and the first numtaps coefficients of this kernel, scaled by window, are returned. The FIR filter will have linear phase. The type of filter is determined by the value of ‘numtaps` and antisymmetric flag. There are four possible combinations:

  • odd numtaps, antisymmetric is False, type I filter is produced

  • even numtaps, antisymmetric is False, type II filter is produced

  • odd numtaps, antisymmetric is True, type III filter is produced

  • even numtaps, antisymmetric is True, type IV filter is produced

Magnitude response of all but type I filters are subjects to following constraints:

  • type II – zero at the Nyquist frequency

  • type III – zero at zero and Nyquist frequencies

  • type IV – zero at zero frequency

New in version 0.9.0.

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#

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]

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

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

See Also#

chebwin, kaiser, bartlett, blackman, hamming, hanning References ———- .. [Rc5e098a1a171-1] W. Carrara, R. Goodman, and R. Majewski, “Spotlight Synthetic

Aperture Radar: Signal Processing Algorithms” Pages 512-513, July 1995.

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.

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]

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.

See Also#

stft: Short Time Fourier Transform check_COLA: Check whether the Constant OverLap Add (COLA) constraint

is met

check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met Notes —– In order to enable inversion of an STFT via the inverse STFT with istft, the signal windowing must obey the constraint of “nonzero overlap add” (NOLA): .. math:: sum_{t}w^{2}[n-tH] ne 0 This ensures that the normalization factors that appear in the denominator of the overlap-add reconstruction equation .. math:: x[n]=frac{sum_{t}x_{t}[n]w[n-tH]}{sum_{t}w^{2}[n-tH]} are not zero. The NOLA constraint can be checked with the check_NOLA function. An STFT which has been modified (via masking or otherwise) is not guaranteed to correspond to a exactly realizible signal. This function implements the iSTFT via the least-squares estimation algorithm detailed in [2], which produces a signal that minimizes the mean squared error between the STFT of the returned signal and the modified STFT. .. versionadded:: 0.19.0 References ———- .. [Rb890c6d0cb06-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#

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

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.

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

New in version 0.12.0.

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]

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#

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#

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)