cuSignal API Reference

Convolution

Convolve

cusignal.convolution.convolve.choose_conv_method(in1, in2, mode='full', measure=False)

Find the fastest convolution/correlation method.

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

Parameters
in1array_like

The first argument passed into the convolution function.

in2array_like

The second argument passed into the convolution function.

modestr {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full

The output is the full discrete linear convolution of the inputs. (Default)

valid

The output consists only of those elements that do not rely on the zero-padding.

same

The output is the same size as in1, centered with respect to the ‘full’ output.

measurebool, optional

If True, run and time the convolution of in1 and in2 with both methods and return the fastest. If False (default), predict the fastest method using precomputed values.

Returns
methodstr

A string indicating which convolution method is fastest, either ‘direct’ or ‘fft’

timesdict, optional

A dictionary containing the times (in seconds) needed for each method. This value is only returned if measure=True.

See also

convolve
correlate

Examples

Estimate the fastest method for a given input:

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

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

>>> c = cp.random.randn(1000)
>>> d = cp.random.randn(1000000)
>>> # `method` works with correlate and convolve
>>> corr1 = cusignal.correlate(a, b, mode='same', method=method)
>>> corr2 = cusignal.correlate(c, d, mode='same', method=method)
>>> conv1 = cusignal.convolve(a, b, mode='same', method=method)
>>> conv2 = cusignal.convolve(c, d, mode='same', method=method)
cusignal.convolution.convolve.convolve(in1, in2, mode='full', method='auto', cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f1039050>, autosync=True)

Convolve two N-dimensional arrays.

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

Parameters
in1array_like

First input.

in2array_like

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

modestr {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full

The output is the full discrete linear convolution of the inputs. (Default)

valid

The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.

same

The output is the same size as in1, centered with respect to the ‘full’ output.

methodstr {‘auto’, ‘direct’, ‘fft’}, optional

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

direct

The convolution is determined directly from sums, the definition of convolution.

fft

The Fourier Transform is used to perform the convolution by calling fftconvolve.

auto

Automatically chooses direct or Fourier method based on an estimate of which is faster (default).

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize(). Default is True.

Returns
convolvearray

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

See also

choose_conv_method

chooses the fastest appropriate convolution method

fftconvolve

Notes

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

Examples

Smooth a square pulse using a Hann window:

>>> import cusignal
>>> import cupy as cp
>>> sig = cp.repeat(cp.asarray([0., 1., 0.]), 100)
>>> win = cusignal.hann(50)
>>> filtered = cusignal.convolve(sig, win, mode='same') / cp.sum(win)
>>> import matplotlib.pyplot as plt
>>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True)
>>> ax_orig.plot(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.convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0, cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f1039090>, autosync=True, use_numba=False)

Convolve two 2-dimensional arrays. Convolve in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue. Parameters ———- in1 : array_like

First input.

in2array_like

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

modestr {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output: full

The output is the full discrete linear convolution of the inputs. (Default)

valid

The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.

same

The output is the same size as in1, centered with respect to the ‘full’ output.

boundarystr {‘fill’, ‘wrap’, ‘symm’}, optional

A flag indicating how to handle boundaries: fill

pad input arrays with fillvalue. (default)

wrap

circular boundary conditions.

symm

symmetrical boundary conditions.

fillvaluescalar, optional

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

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize(). Default is True.

use_numbabool, optional

Option to use Numba CUDA kernel or raw CuPy kernel. Raw CuPy can yield performance gains over Numba. Default is False.

Returns
outndarray

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

Examples
Compute the gradient of an image by 2D convolution with a complex Scharr
operator. (Horizontal operator is real, vertical is imaginary.) Use
symmetric boundary condition to avoid creating edges at the image
boundaries.
>>> import cusignal
    ..
>>> import cupy as cp
    ..
>>> from scipy import misc
    ..
>>> ascent = cp.asarray(misc.ascent())
    ..
>>> scharr = cp.array([[ -3-3j, 0-10j,  +3 -3j],
    ..
… [-10+0j, 0+ 0j, +10 +0j],
… [ -3+3j, 0+10j, +3 +3j]]) # Gx + j*Gy
>>> grad = cusignal.convolve2d(ascent, scharr, boundary='symm',                 mode='same')
    ..
>>> import matplotlib.pyplot as plt
    ..
>>> fig, (ax_orig, ax_mag, ax_ang) = plt.subplots(3, 1, figsize=(6, 15))
    ..
>>> ax_orig.imshow(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 N-dimensional arrays using FFT.

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

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

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

Parameters
in1array_like

First input.

in2array_like

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

modestr {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full

The output is the full discrete linear convolution of the inputs. (Default)

valid

The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.

same

The output is the same size as in1, centered with respect to the ‘full’ output. axis : tuple, optional

axesint or array_like of ints or None, optional

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

Returns
outarray

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

Examples

Autocorrelation of white noise is an impulse.

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

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

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

Correlate

cusignal.convolution.correlate.correlate(in1, in2, mode='full', method='auto', cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f10390d0>, autosync=True)

Cross-correlate two N-dimensional arrays.

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

Parameters
in1array_like

First input.

in2array_like

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

modestr {‘full’, ‘valid’, ‘same’}, optional

A string indicating the size of the output:

full

The output is the full discrete linear cross-correlation of the inputs. (Default)

valid

The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.

same

The output is the same size as in1, centered with respect to the ‘full’ output.

methodstr {‘auto’, ‘direct’, ‘fft’}, optional

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

direct

The correlation is determined directly from sums, the definition of correlation.

fft

The Fast Fourier Transform is used to perform the correlation more quickly (only available for numerical arrays.)

auto

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

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize() Default is True.

Returns
correlatearray

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

See also

choose_conv_method

contains more documentation on method.

Notes

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

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

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

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

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

where \(||x||\) is the length of x, \(N = \max(||x||,||y||)\), and \(y_m\) is 0 when m is outside the range of y.

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

Examples

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

>>> import cusignal
>>> import cupy as cp
>>> sig = cp.repeat([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, (cp.asnumpy(ax_orig), cp.asnumpy(ax_noise), \
    cp.asnumpy(ax_corr)) = plt.subplots(3, 1, sharex=True)
>>> ax_orig.plot(sig)
>>> ax_orig.plot(clock, sig[clock], 'ro')
>>> ax_orig.set_title('Original signal')
>>> ax_noise.plot(cp.asnumpy(sig_noise))
>>> ax_noise.set_title('Signal with noise')
>>> ax_corr.plot(cp.asnumpy(corr))
>>> ax_corr.plot(cp.asnumpy(clock), cp.asnumpy(corr[clock]), 'ro')
>>> ax_corr.axhline(0.5, ls=':')
>>> ax_corr.set_title('Cross-correlated with rectangular pulse')
>>> ax_orig.margins(0, 0.1)
>>> fig.tight_layout()
>>> fig.show()
cusignal.convolution.correlate.correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0, cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f1039410>, autosync=True, use_numba=False)

Cross-correlate two 2-dimensional arrays. Cross correlate in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue. Parameters ———- 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 cross-correlation of the inputs. (Default)

valid

The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.

same

The output is the same size as in1, centered with respect to the ‘full’ output.

boundarystr {‘fill’, ‘wrap’, ‘symm’}, optional

A flag indicating how to handle boundaries: fill

pad input arrays with fillvalue. (default)

wrap

circular boundary conditions.

symm

symmetrical boundary conditions.

fillvaluescalar, optional

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

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize() Default is true.

use_numbabool, optional

Option to use Numba CUDA kernel or raw CuPy kernel. Raw CuPy can yield performance gains over Numba. Default is False.

Returns
correlate2dndarray

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

Examples
Use 2D cross-correlation to find the location of a template in a noisy
image:
>>> import cusignal
    ..
>>> import cupy as cp
    ..
>>> from scipy import misc
    ..
>>> face = cp.asarray(misc.face(gray=True) - misc.face(gray=True).mean())
    ..
>>> template = cp.copy(face[300:365, 670:750])  # right eye
    ..
>>> template -= template.mean()
    ..
>>> face = face + cp.random.randn(*face.shape) * 50  # add noise
    ..
>>> corr = cusignal.correlate2d(face, template, boundary='symm',         mode='same')
    ..
>>> y, x = cp.unravel_index(cp.argmax(corr), corr.shape)  # find the match
    ..
>>> import matplotlib.pyplot as plt
    ..
>>> fig, (cp.asnumpy(ax_orig), cp.asnumpy(ax_template),         cp.asnumpy(ax_corr)) = plt.subplots(3, 1, figsize=(6, 15))
    ..
>>> ax_orig.imshow(cp.asnumpy(face), cmap='gray')
    ..
>>> ax_orig.set_title('Original')
    ..
>>> ax_orig.set_axis_off()
    ..
>>> ax_template.imshow(cp.asnumpy(template), cmap='gray')
    ..
>>> ax_template.set_title('Template')
    ..
>>> ax_template.set_axis_off()
    ..
>>> ax_corr.imshow(cp.asnumpy(corr), cmap='gray')
    ..
>>> ax_corr.set_title('Cross-correlation')
    ..
>>> ax_corr.set_axis_off()
    ..
>>> ax_orig.plot(cp.asnumpy(x), cp.asnumpy(y), 'ro')
    ..
>>> fig.show()
    ..

Filtering

Resample

cusignal.filtering.resample.decimate(x, q, n=None, axis=-1, zero_phase=True, cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f7399950>, autosync=True, use_numba=False)

Downsample the signal after applying an anti-aliasing filter. Parameters ———- x : array_like

The signal to be downsampled, as an N-dimensional array.

qint

The downsampling factor.

nint or array_like, optional

The order of the filter (1 less than the length for FIR) to calculate, or the FIR filter coefficients to employ. Defaults to calculating the coefficients for 20 times the downsampling factor.

axisint, optional

The axis along which to decimate.

zero_phasebool, optional

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

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize(). Default is True.

use_numbabool, optional

Option to use Numba CUDA kernel or raw CuPy kernel. Raw CuPy can yield performance gains over Numba. Default is False.

yndarray

The down-sampled signal.

resample : Resample up or down using the FFT method. resample_poly : Resample using polyphase filtering and an FIR filter. Notes —– Only FIR filter types are currently supported in cuSignal.

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

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

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

Parameters
xarray_like

The data to be resampled.

numint

The number of samples in the resampled signal.

tarray_like, optional

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

axisint, optional

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

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

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

domainstring, optional

A string indicating the domain of the input x:

time

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

freq

Consider the input x as frequency-domain.

Returns
resampled_x or (resampled_x, resampled_t)

Either the resampled array, or, if t was given, a tuple containing the resampled array and the corresponding resampled positions.

See also

decimate

Downsample the signal after applying an FIR or IIR filter.

resample_poly

Resample using polyphase filtering and an FIR filter.

Notes

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

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

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

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

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

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

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

Examples

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

>>> import cusignal
>>> import cupy as cp
>>> x = cp.linspace(0, 10, 20, endpoint=False)
>>> y = cp.cos(-x**2/6.0)
>>> f = cusignal.resample(y, 100)
>>> xnew = cp.linspace(0, 10, 100, endpoint=False)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'go-', xnew, f, '.-', 10, 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), cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f7399990>, autosync=True, use_numba=False)

Resample x along the given axis using polyphase filtering.

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

Parameters
xarray_like

The data to be resampled.

upint

The upsampling factor.

downint

The downsampling factor.

axisint, optional

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

windowstring, tuple, or array_like, optional

Desired window to use to design the low-pass filter, or the FIR filter coefficients to employ. See below for details.

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize(). Default is True.

use_numbabool, optional

Option to use Numba CUDA kernel or raw CuPy kernel. Raw CuPy can yield performance gains over Numba. Default is False.

Returns
resampled_xarray

The resampled array.

See also

decimate

Downsample the signal after applying an FIR or IIR filter.

resample

Resample up or down using the FFT method.

Notes

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

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

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

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

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

Examples

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

>>> from scipy import signal
>>> 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(xnew, f_fft, 'b.-', xnew, f_poly, 'r.-')
>>> plt.plot(x, y, 'ko-')
>>> plt.plot(10, y[0], 'bo', 10, 0., 'ro')  # boundaries
>>> plt.legend(['resample', 'resamp_poly', 'data'], loc='best')
>>> plt.show()
cusignal.filtering.resample.upfirdn(h, x, up=1, down=1, axis=-1, cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f7399a10>, autosync=True, use_numba=False)

Upsample, FIR filter, and downsample Parameters ———- h : array_like

1-dimensional FIR (finite-impulse response) filter coefficients.

xarray_like

Input signal array.

upint, optional

Upsampling rate. Default is 1.

downint, optional

Downsampling rate. Default is 1.

axisint, optional

The axis of the input data array along which to apply the linear filter. The filter is applied to each subarray along this axis. Default is -1.

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize(). Default is True.

use_numbabool, optional

Option to use Numba CUDA kernel or raw CuPy kernel. Raw CuPy can yield performance gains over Numba. Default is False.

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.

The algorithm is an implementation of the block diagram shown on page 129 of the Vaidyanathan text [R9ac9eb812332-1] (Figure 4.3-8d). .. [R9ac9eb812332-1] P. P. Vaidyanathan, Multirate Systems and Filter Banks,

Prentice Hall, 1993.

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

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

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

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

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

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

FIR Filters

Filter Design

Resample

cusignal.filter_design.fir_filter_design.cmplx_sort(p)

Sort roots based on magnitude.

Parameters
parray_like

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

Returns
p_sortedndarray

Sorted roots.

indxndarray

Array of indices needed to sort the input p.

Examples

>>> from scipy import signal
>>> vals = [1, 4, 1+1.j, 3]
>>> p_sorted, indx = signal.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)

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.

Returns
h(numtaps,) ndarray

Coefficients of length numtaps FIR filter.

Raises
ValueError

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

See also

firwin2
firls
minimum_phase
remez

Examples

Low-pass from 0 to f:

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

Use a specific window function:

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

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

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

Band-pass:

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

Band-stop:

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

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

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

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

>>> signal.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, “Discrete-Time Signal Processing”, p.475-476.

Peak Finding

Peak Finding

cusignal.peak_finding.peak_finding.argrelextrema(data, comparator, axis=0, order=1)

Calculate the relative extrema of data.

Parameters
datandarray

Array in which to find the relative extrema.

comparatorcallable

Function to use to compare two data points. Should take two arrays as arguments.

axisint, optional

Axis over which to select from data. Default is 0.

orderint, optional

How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.

Returns
extrematuple of ndarrays

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

See also

argrelmin, argrelmax

Examples

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

Calculate the relative maxima of data.

Parameters
datandarray

Array in which to find the relative maxima.

axisint, optional

Axis over which to select from data. Default is 0.

orderint, optional

How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.

Returns
extrematuple of ndarrays

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

See also

argrelextrema, argrelmin, find_peaks

Notes

This function uses argrelextrema with 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 one-dimensional data find_peaks can be used to detect all local maxima, including flat ones.

Examples

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

Calculate the relative minima of data.

Parameters
datandarray

Array in which to find the relative minima.

axisint, optional

Axis over which to select from data. Default is 0.

orderint, optional

How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.

Returns
extrematuple of ndarrays

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

See also

argrelextrema, argrelmax, find_peaks

Notes

This function uses argrelextrema with 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 one-dimensional data find_peaks can be used to detect all local minima, including flat ones, by calling it with negated data.

Examples

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

Window Functions

Windows

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

Return a modified Bartlett-Hann window.

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.barthann(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Bartlett-Hann window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Bartlett-Hann window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.bartlett(M, sym=True)

Return a Bartlett window.

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

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

See also

triang

A triangular window that does not touch zero at the ends

Notes

The Bartlett window is defined as

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

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

References

1

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

2

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

3

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

4

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

5

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.bartlett(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Bartlett window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Bartlett window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.blackman(M, sym=True)

Return a Blackman window.

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

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The Blackman window is defined as

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

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

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

References

1

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

2

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

3

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.blackman(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Blackman window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(cp.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Blackman window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.blackmanharris(M, sym=True)

Return a minimum 4-term Blackman-Harris window.

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.blackmanharris(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Blackman-Harris window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Blackman-Harris window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.bohman(M, sym=True)

Return a Bohman window.

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.bohman(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Bohman window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Bohman window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.boxcar(M, sym=True)

Return a boxcar or rectangular window.

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

Parameters
Mint

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

symbool, optional

Whether the window is symmetric. (Has no effect for boxcar.)

Returns
wndarray

The window, with the maximum value normalized to 1.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.boxcar(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Boxcar window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the boxcar window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.chebwin(M, at, sym=True)

Return a Dolph-Chebyshev window.

Parameters
Mint

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

atfloat

Attenuation (in dB).

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

The window, with the maximum value always normalized to 1

Notes

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

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

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

where

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

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

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

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

References

1

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

2

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

3

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.chebwin(51, at=100)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Dolph-Chebyshev window (100 dB)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Dolph-Chebyshev window (100 dB)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.cosine(M, sym=True)

Return a window with a simple cosine shape.

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

New in version 0.13.0.

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.cosine(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Cosine window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the cosine window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
>>> plt.show()
cusignal.windows.windows.exponential(M, center=None, tau=1.0, sym=True)

Return an exponential (or Poisson) window.

Parameters
Mint

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

centerfloat, optional

Parameter defining the center location of the window function. The default value if not given is center = (M-1) / 2. This parameter must take its default value for symmetric windows.

taufloat, optional

Parameter defining the decay. For center = 0 use tau = -(M-1) / ln(x) if x is the fraction of the window remaining at the end.

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The Exponential window is defined as

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

References

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

Examples

Plot the symmetric window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> M = 51
>>> tau = 3.0
>>> window = cusignal.exponential(M, tau=tau)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Exponential Window (tau=3.0)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -35, 0])
>>> plt.title("Frequency response of the Exponential window (tau=3.0)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

This function can also generate non-symmetric windows:

>>> tau2 = -(M-1) / np.log(0.01)
>>> window2 = signal.exponential(M, 0, tau2, False)
>>> plt.figure()
>>> plt.plot(window2)
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
cusignal.windows.windows.flattop(M, sym=True)

Return a flat top window.

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

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

References

1

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.flattop(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Flat top window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the flat top window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.gaussian(M, std, sym=True)

Return a Gaussian window.

Parameters
Mint

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

stdfloat

The standard deviation, sigma.

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The Gaussian window is defined as

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.gaussian(51, std=7)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title(r"Gaussian window ($\sigma$=7)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title(r"Frequency response of the Gaussian window ($\sigma$=7)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.general_cosine(M, a, sym=True)

Generic weighted sum of cosine terms window

Parameters
Mint

Number of points in the output window

aarray_like

Sequence of weighting coefficients. This uses the convention of being centered on the origin, so these will typically all be positive numbers, not alternating sign.

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

References

1

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

2

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

Examples

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

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

where

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

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

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

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

>>> from cusignal.windows import general_cosine
>>> from cupy.fft import fft, fftshift
>>> import cupy as cp
>>> import matplotlib.pyplot as plt
>>> window = general_cosine(1000, HFT90D, sym=False)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("HFT90D window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 10000) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(cp.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-50/1000, 50/1000, -140, 0])
>>> plt.title("Frequency response of the HFT90D window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
>>> plt.axhline(-90.2, color='red')
>>> plt.show()
cusignal.windows.windows.general_gaussian(M, p, sig, sym=True)

Return a window with a generalized Gaussian shape.

Parameters
Mint

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

pfloat

Shape parameter. p = 1 is identical to gaussian, p = 0.5 is the same shape as the Laplace distribution.

sigfloat

The standard deviation, sigma.

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The generalized Gaussian window is defined as

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

the half-power point is at

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.general_gaussian(51, p=1.5, sig=7)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title(r"Generalized Gaussian window (p=1.5, $\sigma$=7)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title(r"Freq. resp. of the gen. Gaussian "
...           r"window (p=1.5, $\sigma$=7)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.general_hamming(M, alpha, sym=True)

Return a generalized Hamming window.

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

Parameters
Mint

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

alphafloat

The window coefficient, \(\alpha\)

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

See also

hamming, hann

Notes

The generalized Hamming window is defined as

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

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

References

1

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

2

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

3

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

4

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

Examples

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

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> fig1, spatial_plot = plt.subplots()
>>> spatial_plot.set_title("Generalized Hamming Windows")
>>> spatial_plot.set_ylabel("Amplitude")
>>> spatial_plot.set_xlabel("Sample")
>>> fig2, freq_plot = plt.subplots()
>>> freq_plot.set_title("Frequency Responses")
>>> freq_plot.set_ylabel("Normalized magnitude [dB]")
>>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]")
>>> for alpha in [0.75, 0.7, 0.52]:
...     window = general_hamming(41, alpha)
...     spatial_plot.plot(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(freq, response, label="{:.2f}".format(alpha))
>>> freq_plot.legend(loc="upper right")
>>> spatial_plot.legend(loc="upper right")
cusignal.windows.windows.get_window(window, Nx, fftbins=True)

Return a window of a given length and type.

Parameters
windowstring, float, or tuple

The type of window to create. See below for more details.

Nxint

The number of samples in the window.

fftbinsbool, optional

If True (default), create a “periodic” window, ready to use with ifftshift and be multiplied by the result of an FFT (see also fftpack.fftfreq). If False, create a “symmetric” window, for use in filter design.

Returns
get_windowndarray

Returns a window of length Nx and type window

Notes

Window types:

  • ~cusignal.windows.windows.boxcar

  • ~cusignal.windows.windows.triang

  • ~cusignal.windows.windows.blackman

  • ~cusignal.windows.windows.hamming

  • ~cusignal.windows.windows.hann

  • ~cusignal.windows.windows.bartlett

  • ~cusignal.windows.windows.flattop

  • ~cusignal.windows.windows.parzen

  • ~cusignal.windows.windows.bohman

  • ~cusignal.windows.windows.blackmanharris

  • ~cusignal.windows.windows.nuttall

  • ~cusignal.windows.windows.barthann

  • ~cusignal.windows.windows.kaiser (needs beta)

  • ~cusignal.windows.windows.gaussian (needs standard deviation)

  • ~cusignal.windows.windows.general_gaussian

    (needs power, width)

  • ~cusignal.windows.windows.slepian (needs width)

  • ~cusignal.windows.windows.dpss

    (needs normalized half-bandwidth)

  • ~cusignal.windows.windows.chebwin (needs attenuation)

  • ~cusignal.windows.windows.exponential (needs decay scale)

  • ~cusignal.windows.windows.tukey (needs taper fraction)

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

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

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

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

Examples

>>> import cusignal
>>> cusignal.get_window('triang', 7)
array([ 0.125,  0.375,  0.625,  0.875,  0.875,  0.625,  0.375])
>>> cusignal.get_window(('kaiser', 4.0), 9)
array([ 0.08848053,  0.29425961,  0.56437221,  0.82160913,  0.97885093,
        0.97885093,  0.82160913,  0.56437221,  0.29425961])
>>> cusignal.get_window(4.0, 9)
array([ 0.08848053,  0.29425961,  0.56437221,  0.82160913,  0.97885093,
        0.97885093,  0.82160913,  0.56437221,  0.29425961])
cusignal.windows.windows.hamming(M, sym=True)

Return a Hamming window.

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

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The Hamming window is defined as

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

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

References

1

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

2

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

3

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

4

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.hamming(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Hamming window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Hamming window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.hann(M, sym=True)

Return a Hann window.

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

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The Hann window is defined as

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

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

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

References

1

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

2

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

3

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

4

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.hann(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Hann window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(np.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Hann window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.kaiser(M, beta, sym=True)

Return a Kaiser window.

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

Parameters
Mint

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

betafloat

Shape parameter, determines trade-off between main-lobe width and side lobe level. As beta gets large, the window narrows.

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

Notes

The Kaiser window is defined as

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

with

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

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

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

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

beta

Window shape

0

Rectangular

5

Similar to a Hamming

6

Similar to a Hann

8.6

Similar to a Blackman

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

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

References

1

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

2

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

3

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

4

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.kaiser(51, beta=14)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title(r"Kaiser window ($\beta$=14)")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title(r"Frequency response of the Kaiser window ($\beta$=14)")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.nuttall(M, sym=True)

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

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

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

References

1

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

2

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.nuttall(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Nuttall window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Nuttall window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.parzen(M, sym=True)

Return a Parzen 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).

References

1

E. Parzen, “Mathematical Considerations in the Estimation of Spectra”, Technometrics, Vol. 3, No. 2 (May, 1961), pp. 167-190

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.parzen(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Parzen 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 Parzen window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.triang(M, sym=True)

Return a triangular window.

Parameters
Mint

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

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

See also

bartlett

A triangular window that touches zero

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.triang(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Triangular window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = cp.abs(fftshift(A / cp.abs(A).max()))
>>> response = 20 * cp.log10(cp.maximum(response, 1e-10))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the triangular window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")
cusignal.windows.windows.tukey(M, alpha=0.5, sym=True)

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

Parameters
Mint

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

alphafloat, optional

Shape parameter of the Tukey window, representing the fraction of the window inside the cosine tapered region. If zero, the Tukey window is equivalent to a rectangular window. If one, the Tukey window is equivalent to a Hann window.

symbool, optional

When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns
wndarray

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

References

1

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

2

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

Examples

Plot the window and its frequency response:

>>> import cusignal
>>> import cupy as cp
>>> from cupy.fft import fft, fftshift
>>> import matplotlib.pyplot as plt
>>> window = cusignal.tukey(51)
>>> plt.plot(cp.asnumpy(window))
>>> plt.title("Tukey window")
>>> plt.ylabel("Amplitude")
>>> plt.xlabel("Sample")
>>> plt.ylim([0, 1.1])
>>> plt.figure()
>>> A = fft(window, 2048) / (len(window)/2.0)
>>> freq = cp.linspace(-0.5, 0.5, len(A))
>>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max())))
>>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response))
>>> plt.axis([-0.5, 0.5, -120, 0])
>>> plt.title("Frequency response of the Tukey window")
>>> plt.ylabel("Normalized magnitude [dB]")
>>> plt.xlabel("Normalized frequency [cycles per sample]")

Waveforms

Waveform Generation

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

Frequency-swept cosine generator.

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

Parameters
tarray_like

Times at which to evaluate the waveform.

f0float

Frequency (e.g. Hz) at time t=0.

t1float

Time at which f1 is specified.

f1float

Frequency (e.g. Hz) of the waveform at time t1.

method{‘linear’, ‘quadratic’, ‘logarithmic’, ‘hyperbolic’}, optional

Kind of frequency sweep. If not given, linear is assumed. See Notes below for more details.

phifloat, optional

Phase offset, in degrees. Default is 0.

vertex_zerobool, optional

This parameter is only used when method is ‘quadratic’. It determines whether the vertex of the parabola that is the graph of the frequency is at t=0 or t=t1.

Returns
yndarray

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

Examples

The following will be used in the examples:

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

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

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

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

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

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

>>> w = chirp(t, f0=1500, f1=250, t1=10, method='quadratic')
>>> ff, tt, Sxx = spectrogram(w, fs=fs, noverlap=256, nperseg=512,
...                           nfft=2048)
>>> plt.pcolormesh(cp.asnumpy(tt), cp.asnumpy(ff[:513]),
                   cp.asnumpy(Sxx[:513]), cmap='gray_r')
>>> plt.title('Quadratic Chirp, f(0)=1500, f(10)=250')
>>> plt.xlabel('t (sec)')
>>> plt.ylabel('Frequency (Hz)')
>>> plt.grid()
>>> plt.show()
cusignal.waveforms.waveforms.gausspulse(t, fc=1000, bw=0.5, bwr=- 6, tpr=- 60, retquad=False, retenv=False)

Return a Gaussian modulated sinusoid:

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

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

Parameters
tndarray or the string ‘cutoff’

Input array.

fcint, optional

Center frequency (e.g. Hz). Default is 1000.

bwfloat, optional

Fractional bandwidth in frequency domain of pulse (e.g. Hz). Default is 0.5.

bwrfloat, optional

Reference level at which fractional bandwidth is calculated (dB). Default is -6.

tprfloat, optional

If t is ‘cutoff’, then the function returns the cutoff time for when the pulse amplitude falls below tpr (in dB). Default is -60.

retquadbool, optional

If True, return the quadrature (imaginary) as well as the real part of the signal. Default is False.

retenvbool, optional

If True, return the envelope of the signal. Default is False.

Returns
yIndarray

Real part of signal. Always returned.

yQndarray

Imaginary part of signal. Only returned if retquad is True.

yenvndarray

Envelope of signal. Only returned if retenv is True.

See also

cusignal.morlet

Examples

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

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

Return a periodic square-wave waveform.

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

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

Parameters
tarray_like

The input time array.

dutyarray_like, optional

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

Returns
yndarray

Output array containing the square waveform.

Examples

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

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

A pulse-width modulated sine wave:

>>> plt.figure()
>>> sig = cp.sin(2 * cp.pi * t)
>>> pwm = cusignal.square(2 * cp.pi * 30 * t, duty=(sig + 1)/2)
>>> plt.subplot(2, 1, 1)
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(sig))
>>> plt.subplot(2, 1, 2)
>>> plt.plot(cp.asnumpy(t), cp.asnumpy(pwm))
>>> plt.ylim(-1.5, 1.5)
cusignal.waveforms.waveforms.unit_impulse(shape, idx=None, dtype=<class 'float'>)

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

Parameters
shapeint or tuple of int

Number of samples in the output (1-D), or a tuple that represents the shape of the output (N-D).

idxNone or int or tuple of int or ‘mid’, optional

Index at which the value is 1. If None, defaults to the 0th element. If idx='mid', the impulse will be centered at shape // 2 in all dimensions. If an int, the impulse will be at idx in all dimensions.

dtypedata-type, optional

The desired data-type for the array, e.g., numpy.int8. Default is numpy.float64.

Returns
yndarray

Output array containing an impulse signal.

Notes

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

Examples

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

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

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

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

2-dimensional impulse, centered:

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

Impulse at (2, 2), using broadcasting:

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

Spectrum Analysis

Spectral

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

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

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

Parameters
xarray_like

Time series of measurement values

yarray_like

Time series of measurement values

fsfloat, optional

Sampling frequency of the x and y time series. Defaults to 1.0.

windowstr or tuple or array_like, optional

Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.

npersegint, optional

Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.

noverlap: int, optional

Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None.

nfftint, optional

Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.

detrendstr or function or False, optional

Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.

axisint, optional

Axis along which the coherence is computed for both inputs; the default is over the last axis (i.e. axis=-1).

Returns
fndarray

Array of sample frequencies.

Cxyndarray

Magnitude squared coherence of x and y.

See also

periodogram

Simple, optionally modified periodogram

lombscargle

Lomb-Scargle periodogram for unevenly sampled data

welch

Power spectral density by Welch’s method.

csd

Cross spectral density by Welch’s method.

Notes

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

References

1

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

2

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

Examples

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

Generate two test signals with some common features.

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

Compute and plot the coherence.

>>> f, Cxy = cusignal.coherence(x, y, fs, nperseg=1024)
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Cxy))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Coherence')
>>> plt.show()
cusignal.spectral_analysis.spectral.csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=- 1, average='mean')

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

Parameters
xarray_like

Time series of measurement values

yarray_like

Time series of measurement values

fsfloat, optional

Sampling frequency of the x and y time series. Defaults to 1.0.

windowstr or tuple or array_like, optional

Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.

npersegint, optional

Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.

noverlap: int, optional

Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None.

nfftint, optional

Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.

detrendstr or function or False, optional

Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.

return_onesidedbool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

scaling{ ‘density’, ‘spectrum’ }, optional

Selects between computing the cross spectral density (‘density’) where Pxy has units of V**2/Hz and computing the cross spectrum (‘spectrum’) where Pxy has units of V**2, if x and y are measured in V and fs is measured in Hz. Defaults to ‘density’

axisint, optional

Axis along which the CSD is computed for both inputs; the default is over the last axis (i.e. axis=-1).

average{ ‘mean’, ‘median’ }, optional

Method to use when averaging periodograms. Defaults to ‘mean’.

Returns
fndarray

Array of sample frequencies.

Pxyndarray

Cross spectral density or cross power spectrum of x,y.

See also

periodogram

Simple, optionally modified periodogram

lombscargle

Lomb-Scargle periodogram for unevenly sampled data

welch

Power spectral density by Welch’s method. [Equivalent to csd(x,x)]

coherence

Magnitude squared coherence by Welch’s method.

Notes

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

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

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

References

1

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

2

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

Examples

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

Generate two test signals with some common features.

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

Compute and plot the magnitude of the cross spectral density.

>>> f, Pxy = cusignal.csd(x, y, fs, nperseg=1024)
>>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.abs(Pxy)))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('CSD [V**2/Hz]')
>>> plt.show()
cusignal.spectral_analysis.spectral.lombscargle(x, y, freqs, precenter=False, normalize=False, cp_stream=<cupy.cuda.stream.Stream object at 0x7f88f0fcbb90>, autosync=True, use_numba=False)

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

Sample times.

yarray_like

Measurement values.

freqsarray_like

Angular frequencies for output periodogram.

precenterbool, optional

Pre-center amplitudes by subtracting the mean.

normalizebool, optional

Compute normalized periodogram.

cp_streamCuPy stream, optional

Option allows upfirdn to run in a non-default stream. The use of multiple non-default streams allow multiple kernels to run concurrently. Default is cp.cuda.stream.Stream(null=True) or default stream.

autosyncbool, optional

Option to automatically synchronize cp_stream. This will block the host code until kernel is finished on the GPU. Setting to false will allow asynchronous operation but might required manual synchronize later cp_stream.synchronize(). Default is True.

use_numbabool, optional

Option to use Numba CUDA kernel or raw CuPy kernel. Raw CuPy can yield performance gains over Numba. Default is False.

Returns
pgramarray_like

Lomb-Scargle periodogram.

Raises
ValueError

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

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

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

2

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

3

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

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

Estimate power spectral density using a periodogram.

Parameters
xarray_like

Time series of measurement values

fsfloat, optional

Sampling frequency of the x time series. Defaults to 1.0.

windowstr or tuple or array_like, optional

Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to ‘boxcar’.

nfftint, optional

Length of the FFT used. If None the length of x will be used.

detrendstr or function or False, optional

Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.

return_onesidedbool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

scaling{ ‘density’, ‘spectrum’ }, optional

Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’

axisint, optional

Axis along which the periodogram is computed; the default is over the last axis (i.e. axis=-1).

Returns
fndarray

Array of sample frequencies.

Pxxndarray

Power spectral density or power spectrum of x.

See also

welch

Estimate power spectral density using Welch’s method

lombscargle

Lomb-Scargle periodogram for unevenly sampled data

Examples

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

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

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

Compute and plot the power spectral density.

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

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

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

Now compute and plot the power spectrum.

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

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

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

Compute a spectrogram with consecutive Fourier transforms.

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

Parameters
xarray_like

Time series of measurement values

fsfloat, optional

Sampling frequency of the x time series. Defaults to 1.0.

windowstr or tuple or array_like, optional

Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Tukey window with shape parameter of 0.25.

npersegint, optional

Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.

noverlapint, optional

Number of points to overlap between segments. If None, noverlap = nperseg // 8. Defaults to None.

nfftint, optional

Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.

detrendstr or function or False, optional

Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.

return_onesidedbool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

scaling{ ‘density’, ‘spectrum’ }, optional

Selects between computing the power spectral density (‘density’) where Sxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Sxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’.

axisint, optional

Axis along which the spectrogram is computed; the default is over the last axis (i.e. axis=-1).

modestr, optional

Defines what kind of return values are expected. Options are [‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’]. ‘complex’ is equivalent to the output of stft with no padding or boundary extension. ‘magnitude’ returns the absolute magnitude of the STFT. ‘angle’ and ‘phase’ return the complex angle of the STFT, with and without unwrapping, respectively.

Returns
fndarray

Array of sample frequencies.

tndarray

Array of segment times.

Sxxndarray

Spectrogram of x. By default, the last axis of Sxx corresponds to the segment times.

See also

periodogram

Simple, optionally modified periodogram

lombscargle

Lomb-Scargle periodogram for unevenly sampled data

welch

Power spectral density by Welch’s method.

csd

Cross spectral density by Welch’s method.

Notes

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

New in version 0.16.0.

References

1

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

Examples

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

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

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

Compute and plot the spectrogram.

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

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

>>> f, t, Sxx = cusignal.spectrogram(x, fs, return_onesided=False)
>>> plt.pcolormesh(cp.asnumpy(t), cp.fft.fftshift(f),         cp.fft.fftshift(Sxx, axes=0))
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()
cusignal.spectral_analysis.spectral.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend=False, return_onesided=True, boundary='zeros', padded=True, axis=- 1)

Compute the Short Time Fourier Transform (STFT).

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

Parameters
xarray_like

Time series of measurement values

fsfloat, optional

Sampling frequency of the x time series. Defaults to 1.0.

windowstr or tuple or array_like, optional

Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.

npersegint, optional

Length of each segment. Defaults to 256.

noverlapint, optional

Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None. When specified, the COLA constraint must be met (see Notes below).

nfftint, optional

Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.

detrendstr or function or False, optional

Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to False.

return_onesidedbool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

boundarystr or None, optional

Specifies whether the input signal is extended at both ends, and how to generate the new values, in order to center the first windowed segment on the first input point. This has the benefit of enabling reconstruction of the first input point when the employed window function starts at zero. Valid options are ['even', 'odd', 'constant', 'zeros', None]. Defaults to ‘zeros’, for zero padding extension. I.e. [1, 2, 3, 4] is extended to [0, 1, 2, 3, 4, 0] for nperseg=3.

paddedbool, optional

Specifies whether the input signal is zero-padded at the end to make the signal fit exactly into an integer number of window segments, so that all of the signal is included in the output. Defaults to True. Padding occurs after boundary extension, if boundary is not None, and padded is True, as is the default.

axisint, optional

Axis along which the STFT is computed; the default is over the last axis (i.e. axis=-1).

Returns
fndarray

Array of sample frequencies.

tndarray

Array of segment times.

Zxxndarray

STFT of x. By default, the last axis of Zxx corresponds to the segment times.

See also

welch

Power spectral density by Welch’s method.

spectrogram

Spectrogram by Welch’s method.

csd

Cross spectral density by Welch’s method.

lombscargle

Lomb-Scargle periodogram for unevenly sampled data

Notes

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

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

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

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

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

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

References

1

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

2

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

Examples

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

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

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

Compute and plot the STFT’s magnitude.

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

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

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

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

Parameters
events1D array_like

An array of time points containing the timing of the events.

periodfloat or array_like

The period of the signal that the events should synchronize to. The period is in the same units as events. It can also be an array of periods, in which case the outputs are arrays of the same length.

Returns
strengthfloat or 1D array

The strength of the synchronization. 1.0 is perfect synchronization and 0.0 is no synchronization. If period is an array, this is also an array with each element containing the vector strength at the corresponding period.

phasefloat or array

The phase that the events are most strongly synchronized to in radians. If period is an array, this is also an array with each element containing the phase for the corresponding period.

References

van Hemmen, JL, Longtin, A, and Vollmayr, AN. Testing resonating vector

strength: Auditory system, electric fish, and noise. Chaos 21, 047508 (2011); :doi:`10.1063/1.3670512`.

van Hemmen, JL. Vector strength after Goldberg, Brown, and von Mises:

biological and mathematical perspectives. Biol Cybern. 2013 Aug;107(4):385-96. :doi:`10.1007/s00422-013-0561-7`.

van Hemmen, JL and Vollmayr, AN. Resonating vector strength: what happens

when we vary the “probing” frequency while keeping the spike times fixed. Biol Cybern. 2013 Aug;107(4):491-94. :doi:`10.1007/s00422-013-0560-8`.

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

Estimate power spectral density using Welch’s method.

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

Parameters
xarray_like

Time series of measurement values

fsfloat, optional

Sampling frequency of the x time series. Defaults to 1.0.

windowstr or tuple or array_like, optional

Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.

npersegint, optional

Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.

noverlapint, optional

Number of points to overlap between segments. If None, noverlap = nperseg // 2. Defaults to None.

nfftint, optional

Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.

detrendstr or function or False, optional

Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.

return_onesidedbool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

scaling{ ‘density’, ‘spectrum’ }, optional

Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’

axisint, optional

Axis along which the periodogram is computed; the default is over the last axis (i.e. axis=-1).

average{ ‘mean’, ‘median’ }, optional

Method to use when averaging periodograms. Defaults to ‘mean’.

Returns
fndarray

Array of sample frequencies.

Pxxndarray

Power spectral density or power spectrum of x.

See also

periodogram

Simple, optionally modified periodogram

lombscargle

Lomb-Scargle periodogram for unevenly sampled data

Notes

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

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

New in version 0.12.0.

References

1

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

2

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

Examples

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

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

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

Compute and plot the power spectral density.

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

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

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

Now compute and plot the power spectrum.

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

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

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

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

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

Acoustics

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

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

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

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

nint

Size of Fourier Transform; If none, will use length of input array

axis: int

Direction for cepstrum calculation

cepsndarray

Complex cepstrum result

cusignal.acoustics.cepstrum.cceps_unwrap(x)

Unwrap phase for complex cepstrum calculation; helper function

cusignal.acoustics.cepstrum.rceps(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 ———- 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

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 high-pass qmf filter from low-pass

Parameters
hkarray_like

Coefficients of high-pass filter.

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

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

It models the function:

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

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

Parameters
pointsint

Number of points in vector. Will be centered around 0.

ascalar

Width parameter of the wavelet.

Returns
vector(N,) ndarray

Array of length points in shape of ricker curve.

Examples

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

B-splines

B-splines

cusignal.bsplines.bsplines.cspline1d(signal, lamb=0.0)

Compute cubic spline coefficients for rank-1 array.

Find the cubic spline coefficients for a 1-D signal assuming mirror-symmetric boundary conditions. To obtain the signal back from the spline representation mirror-symmetric-convolve these coefficients with a length 3 FIR window [1.0, 4.0, 1.0]/ 6.0 .

Parameters
signalndarray

A rank-1 array representing samples of a signal.

lambfloat, optional

Smoothing coefficient, default is 0.0.

Returns
cndarray

Cubic spline coefficients.

cusignal.bsplines.bsplines.cubic(x)

A cubic B-spline.

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

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

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

Parameters
nint

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

References

1

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

cusignal.bsplines.bsplines.quadratic(x)

A quadratic B-spline.

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

Utilities

Array Tools

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

Return populated shared memory between GPU and CPU.

Parameters
datacupy.ndarray or numpy.ndarray

The array to be copied to shared buffer

strides: int or None
order: char
streamint

Stream number (0 for default)

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

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

Parameters
shapendarray.shape

Size of shared memory allocation

dtypecupy.dtype or numpy.dtype

Data type of allocation

strides: int or None
order: char
streamint

Stream number (0 for default)

portablebool
wcbool

FFTPack Helper

cusignal.utils.fftpack_helper.next_fast_len(target)

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

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

Parameters
targetint

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

Returns
outint

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

Notes

New in version 0.18.0.

Examples

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

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

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

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

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

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