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 realworld 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 zeropadding.
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 Ndimensional 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 zeropadding. 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 Ndimensional 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 precomputed values (choose_conv_method can also measure realworld 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 1dimensional 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 zeropadding. 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 1dimensional array containing a subset of the discrete linear convolution of in1 with in2.
See also
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.convolve1d3o(a,b)

cusignal.convolution.convolve.
convolve1d3o
(in1, in2, mode='valid', method='direct')¶ Convolve a 1dimensional 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 zeropadding. 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 1dimensional array containing a subset of the discrete linear convolution of in1 with in2.
See also
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 2dimensional 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 zeropadding. 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 2dimensional 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([[ 33j, 010j, +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.asarray(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 Ndimensional 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 zeropadding. 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 Ndimensional 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(np.arange(len(sig)+1,len(sig)), autocorr) >>> ax_mag.set_title('Autocorrelation') >>> fig.tight_layout() >>> fig.show()
Gaussian blur implemented using FFT convolution. Notice the dark borders around the image, due to the zeropadding 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.fftconvolve(face, kernel, mode='same')
>>> fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1, ... figsize=(6, 15)) >>> ax_orig.imshow(face, cmap='gray') >>> ax_orig.set_title('Original') >>> ax_orig.set_axis_off() >>> ax_kernel.imshow(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')¶ Crosscorrelate two Ndimensional arrays.
Crosscorrelate 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 crosscorrelation of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zeropadding. 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 Ndimensional array containing a subset of the discrete linear crosscorrelation of in1 with in2.
See also
choose_conv_method
contains more documentation on method.
Notes
The correlation z of two ddimensional 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 1D arrays and
z = correlate(x, y, 'full')
then\[z[k] = (x * y)(k  N + 1) = \sum_{l=0}^{x1}x_l y_{lk+N1}^{*}\]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 crosscorrelation, 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('Crosscorrelated 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)¶ Crosscorrelate two 2dimensional 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 crosscorrelation of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zeropadding. 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 2dimensional array containing a subset of the discrete linear crosscorrelation of in1 with in2.
Examples
Use 2D crosscorrelation 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(‘Crosscorrelation’) >>> ax_corr.set_axis_off() >>> ax_orig.plot(cp.asnumpy(x), cp.asnumpy(y), ‘ro’) >>> fig.show()
Estimation¶
Kalman Filter¶

class
cusignal.estimation.filters.
KalmanFilter
(dim_x, dim_z, dim_u=0, points=1, dtype=<class 'numpy.float32'>)¶ This is a multipoint Kalman Filter implementation of https://github.com/rlabbe/filterpy/blob/master/filterpy/kalman/kalman_filter.py, 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. 208212. (2006)
 2
Roger Labbe. “Kalman and Bayesian Filters in Python” https://github.com/rlabbe/KalmanandBayesianFiltersinPython
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
 xarray(points, dim_x, 1)
Current state estimate. Any call to update() or predict() updates this variable.
 Parray(points, dim_x, dim_x)
Current state covariance matrix. Any call to update() or predict() updates this variable.
 zarray(points, dim_z, 1)
Last measurement used in update(). Read only.
 Rarray(points, dim_z, dim_z)
Measurement noise matrix
 Qarray(points, dim_x, dim_x)
Process noise matrix
 Farray(points, dim_x, dim_x)
State Transition matrix
 Harray(points, dim_z, dim_x)
Measurement function
 _alpha_sqfloat (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 antialiasing filter. Parameters ——— x : array_like
The signal to be downsampled, as an Ndimensional 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 downsampled signal.
 See Also
 resampleResample up or down using the FFT method.
 resample_polyResample 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 timedomain. (Default)
freq
Consider the input x as frequencydomain.
 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 Fourierdomain window that tapers the Fourier spectrum before zeropadding to alleviate ringing in the resampled values for sampled signals you didn’t intend to be interpreted as bandlimited.
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 lowfrequency 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
todx * 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 zerophase lowpass 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 lowpass 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
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 lowpass 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 zerophase 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
todx * 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
1dimensional FIR (finiteimpulse 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.38d).
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 sampleandhold 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 onedimension 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 1D sequence.
 xarray_like
An Ndimensional 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 Ndimensional 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.
Examples¶
>>> from scipy import signal
>>> import cupy as cp
>>> import cusignal
>>> [b, a] = signal.butter(3, 0.5)
>>> b = cp.asarray(b)
>>> x = cp.random.randn(2**8)
>>> y = cusignal.firfilter(b, x)
Channelizer¶
Polyphase channelize signal into n channels
Parameters¶
 xarray_like
The input data to be channelized
 harray_like
The 1D 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 1D 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=1.0, 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
Lowpass 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.56607041e04, 9.99286786e01, 3.56607041e04])
Highpass (‘stop’ from 0 to f):
>>> cusignal.firwin(numtaps, f, pass_zero=False) array([0.00859313, 0.98281375, 0.00859313])
Bandpass:
>>> f1, f2 = 0.1, 0.2 >>> cusignal.firwin(numtaps, [f1, f2], pass_zero=False) array([ 0.06301614, 0.88770441, 0.06301614])
Bandstop:
>>> cusignal.firwin(numtaps, [f1, f2]) array([0.00801395, 1.0160279 , 0.00801395])
Multiband (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])
Multiband (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.
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
 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.
 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.
 betafloat
The beta parameter to be used in the formula for a Kaiser window.
Oppenheim, Schafer, “DiscreteTime Signal Processing”, p.475476.
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 onedimensional.
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([0, 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, 0, 2]), array([0, 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 onedimensional.
See also
argrelextrema
,argrelmin
,find_peaks
Notes
This function uses argrelextrema with np.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 onedimensional 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([0, 3, 6]),) >>> y = cp.array([[1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]]) ... >>> argrelmax(y, axis=1) (array([0, 0, 2]), array([1 ,3, 0]))

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 onedimensional.
See also
argrelextrema
,argrelmax
,find_peaks
Notes
This function uses argrelextrema with np.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 onedimensional 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, 7]),) >>> y = cp.array([[1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]]) ... >>> argrelmin(y, axis=1) (array([0, 0, 2]), array([0, 2, 1]))
Window Functions¶
Windows¶

cusignal.windows.windows.
barthann
(M, sym=True)¶ Return a modified BartlettHann 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("BartlettHann 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 BartlettHann 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}{M1} \left( \frac{M1}{2}  \leftn  \frac{M1}{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, 116, 1950.
 2
E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 109110.
 3
A.V. Oppenheim and R.W. Schafer, “DiscreteTime Signal Processing”, PrenticeHall, 1999, pp. 468471.
 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 falloff. 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 falloff 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. DiscreteTime Signal Processing. Upper Saddle River, NJ: PrenticeHall, 1999, pp. 468471.
 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): 5183. :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, 1e10)) >>> 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 4term BlackmanHarris 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("BlackmanHarris 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 BlackmanHarris 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 DolphChebyshev 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 DolphChebyshev 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) <= M1. A is the attenuation in decibels (at).
The time domain window is then generated using the IFFT, so poweroftwo 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 sidelobe level”, Proceedings of the IEEE, Vol. 34, Issue 6
 2
Peter Lynch, “The DolphChebyshev 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("DolphChebyshev 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 DolphChebyshev 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 = (M1) / 2
. This parameter must take its default value for symmetric windows. taufloat, optional
Parameter defining the decay. For
center = 0
usetau = (M1) / ln(x)
ifx
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^{ncenter / \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 nonsymmetric windows:
>>> tau2 = (M1) / 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 5thorder 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/0387286667`.
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. 8491, 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 flattop windows”, February 15, 2002 https://holometer.fnal.gov/GH_FFT.pdf
Examples
Heinzel describes a flattop 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, 1e10)) >>> 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 halfpower 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).
Notes
The generalized Hamming window is defined as
\[w(n) = \alpha  \left(1  \alpha\right) \cos\left(\frac{2\pi{n}}{M1}\right) \qquad 0 \leq n \leq M1\]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, “Sentinel1 Level 1 Detailed Algorithm Definition”, https://sentinel.esa.int/documents/247904/1877131/Sentinel1Level1DetailedAlgorithmDefinition
 4
Matthieu Bourbigot ESA, “Sentinel1 Product Definition”, https://sentinel.esa.int/documents/247904/1877131/Sentinel1ProductDefinition
Examples
The Sentinel1A/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 halfbandwidth)
~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 nonzero 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}}{M1}\right) \qquad 0 \leq n \leq M1\]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. 109110.
 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 sinesquared 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}}{M1}\right) \qquad 0 \leq n \leq M1\]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. 106108.
 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, 1e10)) >>> 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 tradeoff between mainlobe 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}{(M1)^2}} \right)/I_0(\beta)\]with
\[\quad \frac{M1}{2} \leq n \leq \frac{M1}{2},\]where \(I_0\) is the modified zerothorder 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 218285. John Wiley and Sons, New York, (1966).
 2
E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 177178.
 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. 5183, 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 4term BlackmanHarris 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. 8491, 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 flattop 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.
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, 1e10)) >>> 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): 5183. :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)¶ Frequencyswept 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.
 Returns
 yndarray
A numpy array containing the signal evaluated at t with the requested timevarying frequency. More precisely, the function returns
cos(phase + (pi/180)*phi)
where phase is the integral (from 0 to t) of2*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 (inphase 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.
square
(t, duty=0.5)¶ Return a periodic squarewave waveform.
The square wave has a period
2*pi
, has value +1 from 0 to2*pi*duty
and 1 from2*pi*duty
to2*pi
. duty must be in the interval [0,1].Note that this is not bandlimited. 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(t, cp.asnumpy(cusignal.square(2 * cp.pi * 5 * t))) >>> plt.ylim(2, 2)
A pulsewidth 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 (1D), or a tuple that represents the shape of the output (ND).
 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 atshape // 2
in all dimensions. If an int, the impulse will be at idx in all dimensions. dtypedatatype, optional
The desired datatype for the array, e.g.,
numpy.int8
. Default isnumpy.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[n2]\)):
>>> cusignal.unit_impulse(7, 2) array([ 0., 0., 1., 0., 0., 0., 0.])
2dimensional 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 discretetime 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 DFTeven 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
LombScargle 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. 7073, 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 DFTeven 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 onesided spectrum for real data. If False return a twosided spectrum. Defaults to True, but for complex data, a twosided 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
LombScargle 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 zeropadded 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. 7073, 1967.
 2
Rabiner, Lawrence R., and B. Gold. “Theory and Application of Digital Signal Processing” PrenticeHall, pp. 414419, 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.
lombscargle
(x, y, freqs)¶ Computes the LombScargle periodogram. The LombScargle 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 onedimensional and will be cast to float64. Parameters
 xarray_like
Sample times.
 yarray_like
Measurement values.
 freqsarray_like
Angular frequencies for output periodogram.
 precenterbool, optional
Precenter amplitudes by subtracting the mean.
 normalizebool, optional
Compute normalized periodogram.
 Returns
 pgramarray_like
LombScargle 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 “Leastsquares frequency analysis of unequally spaced data”, Astrophysics and Space Science, vol 39, pp. 447462, 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. 835853, 1982
 3
R.H.D. Townsend, “Fast calculation of the LombScargle periodogram using graphics processing units.”, The Astrophysical Journal Supplement Series, vol 191, pp. 247253, 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 LombScargle 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 DFTeven 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 onesided spectrum for real data. If False return a twosided spectrum. Defaults to True, but for complex data, a twosided 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
LombScargle 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([1e7, 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([1e4, 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 DFTeven 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 onesided spectrum for real data. If False return a twosided spectrum. Defaults to True, but for complex data, a twosided 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
LombScargle 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 “DiscreteTime 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 DFTeven 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 onesided spectrum for real data. If False return a twosided spectrum. Defaults to True, but for complex data, a twosided 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]
fornperseg=3
. paddedbool, optional
Specifies whether the input signal is zeropadded 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
LombScargle 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) % (npersegnoverlap) == 0
). The padded argument may be used to accomplish this.Given a timedomain 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[ntH]\]The overlapadd (OLA) reconstruction equation is given by
\[x[n]=\frac{\sum_{t}x_{t}[n]w[ntH]}{\sum_{t}w^{2}[ntH]}\]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 “DiscreteTime Signal Processing”, Prentice Hall, 1999.
 2
Daniel W. Griffin, Jae S. Lim “Signal Estimation from Modified ShortTime 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):38596. :doi:`10.1007/s0042201305617`.
 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):49194. :doi:`10.1007/s0042201305608`.

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 DFTeven 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 onesided spectrum for real data. If False return a twosided spectrum. Defaults to True, but for complex data, a twosided 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
LombScargle 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. 7073, 1967.
 2
M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika, vol. 37, pp. 116, 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.5e3, 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.5e3, 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 zerophase 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
 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.
 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(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.
qmf
(hk)¶ Return highpass qmf filter from lowpass
 Parameters
 hkarray_like
Coefficients of highpass 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()
Bsplines¶
Bsplines¶

cusignal.bsplines.bsplines.
cubic
(x)¶ A cubic Bspline.
This is a special case of bspline, and equivalent to
bspline(x, 3)
.

cusignal.bsplines.bsplines.
gauss_spline
(x, n)¶ Gaussian approximation to Bspline 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 BSplines. 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 Bspline.
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
 sizeint or tuple of ints
Output shape.
 dtypedatatype
Output data type.
 Returns
 outndarray
Pinned memory numpy array.
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
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 zeropadding, 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 5smooth numbers, regular numbers, or Hamming numbers.)
 Parameters
 targetint
Length to start searching from. Must be a positive integer.
 Returns
 outint
The first 5smooth 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)
Zeropadding to the next 5smooth 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 5smooth 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.
 dtypedatatype, 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 bytesize of dtype.
 Returns
 ——
 outndarray
An 1dimensional 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 nonnegative 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 1dimensional array containing unpacked binary data.

cusignal.io.reader.
unpack_bin
(binary, dtype, endianness='L')¶ Unpack binary file. If endianness is bigendian, it my be converted to little endian for NVIDIA GPU compatibility.
 Parameters
 binaryndarray
The binary array to be unpack.
 dtypedatatype, optional
Any object that can be interpreted as a numpy data type.
 endianness{‘L’, ‘B’}, optional
Data set byte order
 Returns
 outndarray
An 1dimensional 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 1dimensional 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 1dimensional 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.