cuSignal API Reference#
Convolution#
Convolve#
- cusignal.convolution.convolve.choose_conv_method(in1, in2, mode='full', measure=False)#
Find the fastest convolution/correlation method.
This primarily exists to be called during the
method='auto'
option in convolve and correlate, but can also be used when performing many convolutions of the same input shapes and dtypes, determining which method to use for all of them, either to avoid the overhead of the ‘auto’ option or to use accurate real-world measurements.- Parameters:
- in1array_like
The first argument passed into the convolution function.
- in2array_like
The second argument passed into the convolution function.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- measurebool, optional
If True, run and time the convolution of in1 and in2 with both methods and return the fastest. If False (default), predict the fastest method using precomputed values.
- Returns:
- methodstr
A string indicating which convolution method is fastest, either ‘direct’ or ‘fft’
- timesdict, optional
A dictionary containing the times (in seconds) needed for each method. This value is only returned if
measure=True
.
See also
convolve
correlate
Examples
Estimate the fastest method for a given input:
>>> import cusignal >>> import cupy as cp >>> a = cp.random.randn(1000) >>> b = cp.random.randn(1000000) >>> method = cusignal.choose_conv_method(a, b, mode='same') >>> method 'fft'
This can then be applied to other arrays of the same dtype and shape:
>>> c = cp.random.randn(1000) >>> d = cp.random.randn(1000000) >>> # `method` works with correlate and convolve >>> corr1 = cusignal.correlate(a, b, mode='same', method=method) >>> corr2 = cusignal.correlate(c, d, mode='same', method=method) >>> conv1 = cusignal.convolve(a, b, mode='same', method=method) >>> conv2 = cusignal.convolve(c, d, mode='same', method=method)
- cusignal.convolution.convolve.convolve(in1, in2, mode='full', method='auto')#
Convolve two N-dimensional arrays.
Convolve in1 and in2, with the output size determined by the mode argument.
- Parameters:
- in1array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- methodstr {‘auto’, ‘direct’, ‘fft’}, optional
A string indicating which method to use to calculate the convolution.
direct
The convolution is determined directly from sums, the definition of convolution.
fft
The Fourier Transform is used to perform the convolution by calling fftconvolve.
auto
Automatically chooses direct or Fourier method based on an estimate of which is faster (default).
- Returns:
- convolvearray
An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.
See also
choose_conv_method
chooses the fastest appropriate convolution method
fftconvolve
Notes
By default, convolve and correlate use
method='auto'
, which calls choose_conv_method to choose the fastest method using pre-computed values (choose_conv_method can also measure real-world timing with a keyword argument). Because fftconvolve relies on floating point numbers, there are certain constraints that may force method=direct (more detail in choose_conv_method docstring).Examples
Smooth a square pulse using a Hann window:
>>> import cusignal >>> import cupy as cp >>> sig = cp.repeat(cp.asarray([0., 1., 0.]), 100) >>> win = cusignal.hann(50) >>> filtered = cusignal.convolve(sig, win, mode='same') / cp.sum(win)
>>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True) >>> ax_orig.plot(cp.asnumpy(sig)) >>> ax_orig.set_title('Original pulse') >>> ax_orig.margins(0, 0.1) >>> ax_win.plot(cp.asnumpy(win)) >>> ax_win.set_title('Filter impulse response') >>> ax_win.margins(0, 0.1) >>> ax_filt.plot(cp.asnumpy(filtered)) >>> ax_filt.set_title('Filtered signal') >>> ax_filt.margins(0, 0.1) >>> fig.tight_layout() >>> fig.show()
- cusignal.convolution.convolve.convolve1d2o(in1, in2, mode='valid', method='direct')#
Convolve a 1-dimensional arrays with a 2nd order filter. This results in a second order convolution.
Convolve in1 and in2, with the output size determined by the mode argument.
- Parameters:
- in1array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- methodstr {‘auto’, ‘direct’, ‘fft’}, optional
A string indicating which method to use to calculate the convolution.
direct
The convolution is determined directly from sums, the definition of convolution.
fft
The Fourier Transform is used to perform the convolution by calling fftconvolve.
auto
Automatically chooses direct or Fourier method based on an estimate of which is faster (default).
- Returns:
- outndarray
A 1-dimensional array containing a subset of the discrete linear convolution of in1 with in2.
See also
Examples
Convolution of a 2nd order filter on a 1d signal
>>> import cusignal as cs >>> import numpy as np >>> d = 50 >>> a = np.random.uniform(-1,1,(200)) >>> b = np.random.uniform(-1,1,(d,d)) >>> c = cs.convolve1d2o(a,b)
- cusignal.convolution.convolve.convolve1d3o(in1, in2, mode='valid', method='direct')#
Convolve a 1-dimensional array with a 3rd order filter. This results in a second order convolution.
Convolve in1 and in2, with the output size determined by the mode argument.
- Parameters:
- in1array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- methodstr {‘auto’, ‘direct’, ‘fft’}, optional
A string indicating which method to use to calculate the convolution.
direct
The convolution is determined directly from sums, the definition of convolution.
fft
The Fourier Transform is used to perform the convolution by calling fftconvolve.
auto
Automatically chooses direct or Fourier method based on an estimate of which is faster (default).
- Returns:
- outndarray
A 1-dimensional array containing a subset of the discrete linear convolution of in1 with in2.
See also
Examples
Convolution of a 3rd order filter on a 1d signal
>>> import cusignal as cs >>> import numpy as np >>> d = 50 >>> a = np.random.uniform(-1,1,(200)) >>> b = np.random.uniform(-1,1,(d,d,d)) >>> c = cs.convolve1d3o(a,b)
- cusignal.convolution.convolve.convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0)#
Convolve two 2-dimensional arrays. Convolve in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue. Parameters ———- in1 : array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- boundarystr {‘fill’, ‘wrap’, ‘symm’}, optional
A flag indicating how to handle boundaries:
fill
pad input arrays with fillvalue. (default)
wrap
circular boundary conditions.
symm
symmetrical boundary conditions.
- fillvaluescalar, optional
Value to fill pad input arrays with. Default is 0.
- Returns:
- outndarray
A 2-dimensional array containing a subset of the discrete linear convolution of in1 with in2.
Examples
Compute the gradient of an image by 2D convolution with a complex Scharr operator. (Horizontal operator is real, vertical is imaginary.) Use symmetric boundary condition to avoid creating edges at the image boundaries.
>>> import cusignal >>> import cupy as cp >>> from scipy import misc >>> ascent = cp.asarray(misc.ascent()) >>> scharr = cp.array([[ -3-3j, 0-10j, +3 -3j], ... [-10+0j, 0+ 0j, +10 +0j], ... [ -3+3j, 0+10j, +3 +3j]]) # Gx + j*Gy >>> grad = cusignal.convolve2d(ascent, scharr, boundary='symm', mode='same') >>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_mag, ax_ang) = plt.subplots(3, 1, figsize=(6, 15)) >>> ax_orig.imshow(cp.asnumpy(ascent), cmap='gray') >>> ax_orig.set_title('Original') >>> ax_orig.set_axis_off() >>> ax_mag.imshow(cp.asnumpy(cp.absolute(grad)), cmap='gray') >>> ax_mag.set_title('Gradient magnitude') >>> ax_mag.set_axis_off() >>> ax_ang.imshow(cp.asnumpy(cp.angle(grad)), cmap='hsv') >>> ax_ang.set_title('Gradient orientation') >>> ax_ang.set_axis_off() >>> fig.show()
- cusignal.convolution.convolve.fftconvolve(in1, in2, mode='full', axes=None)#
Convolve two N-dimensional arrays using FFT.
Convolve in1 and in2 using the fast Fourier transform method, with the output size determined by the mode argument.
This is generally much faster than convolve for large arrays (n > ~500), but can be slower when only a few output values are needed, and can only output float arrays (int or object array inputs will be cast to float).
As of v0.19, convolve automatically chooses this method or the direct method based on an estimation of which is faster.
- Parameters:
- in1array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear convolution of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output. axis : tuple, optional
- axesint or array_like of ints or None, optional
Axes over which to compute the convolution. The default is over all axes.
- Returns:
- outarray
An N-dimensional array containing a subset of the discrete linear convolution of in1 with in2.
Examples
Autocorrelation of white noise is an impulse.
>>> import cusignal >>> import cupy as cp >>> import numpy as np >>> sig = cp.random.randn(1000) >>> autocorr = cusignal.fftconvolve(sig, sig[::-1], mode='full')
>>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1) >>> ax_orig.plot(cp.asnumpy(sig)) >>> ax_orig.set_title('White noise') >>> ax_mag.plot(cp.asnumpy(cp.arange(-len(sig)+1,len(sig))), cp.asnumpy(autocorr)) >>> ax_mag.set_title('Autocorrelation') >>> fig.tight_layout() >>> fig.show()
Gaussian blur implemented using FFT convolution. Notice the dark borders around the image, due to the zero-padding beyond its boundaries. The convolve2d function allows for other types of image boundaries, but is far slower.
>>> from scipy import misc >>> face = misc.face(gray=True) >>> kernel = cp.outer(cusignal.gaussian(70, 8), cusignal.gaussian(70, 8)) >>> blurred = cusignal.convolve(face, kernel, mode='same')
>>> fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1, ... figsize=(6, 15)) >>> ax_orig.imshow(face, cmap='gray') >>> ax_orig.set_title('Original') >>> ax_orig.set_axis_off() >>> ax_kernel.imshow(cp.asnumpy(kernel), cmap='gray') >>> ax_kernel.set_title('Gaussian kernel') >>> ax_kernel.set_axis_off() >>> ax_blurred.imshow(cp.asnumpy(blurred), cmap='gray') >>> ax_blurred.set_title('Blurred') >>> ax_blurred.set_axis_off() >>> fig.show()
Correlate#
- cusignal.convolution.correlate.correlate(in1, in2, mode='full', method='auto')#
Cross-correlate two N-dimensional arrays.
Cross-correlate in1 and in2, with the output size determined by the mode argument.
- Parameters:
- in1array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear cross-correlation of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- methodstr {‘auto’, ‘direct’, ‘fft’}, optional
A string indicating which method to use to calculate the correlation.
direct
The correlation is determined directly from sums, the definition of correlation.
fft
The Fast Fourier Transform is used to perform the correlation more quickly (only available for numerical arrays.)
auto
Automatically chooses direct or Fourier method based on an estimate of which is faster (default). See convolve Notes for more detail.
- Returns:
- correlatearray
An N-dimensional array containing a subset of the discrete linear cross-correlation of in1 with in2.
See also
choose_conv_method
contains more documentation on method.
Notes
The correlation z of two d-dimensional arrays x and y is defined as:
z[...,k,...] = sum[..., i_l, ...] x[..., i_l,...] * conj(y[..., i_l - k,...])
This way, if x and y are 1-D arrays and
z = correlate(x, y, 'full')
then\[z[k] = (x * y)(k - N + 1) = \sum_{l=0}^{||x||-1}x_l y_{l-k+N-1}^{*}\]for \(k = 0, 1, ..., ||x|| + ||y|| - 2\)
where \(||x||\) is the length of
x
, \(N = \max(||x||,||y||)\), and \(y_m\) is 0 when m is outside the range of y.method='fft'
only works for numerical arrays as it relies on fftconvolve. In certain cases (i.e., arrays of objects or when rounding integers can lose precision),method='direct'
is always used.Examples
Implement a matched filter using cross-correlation, to recover a signal that has passed through a noisy channel.
>>> import cusignal >>> import cupy as cp >>> sig = cp.repeat(cp.array([0., 1., 1., 0., 1., 0., 0., 1.]), 128) >>> sig_noise = sig + cp.random.randn(len(sig)) >>> corr = cusignal.correlate(sig_noise, cp.ones(128), mode='same') / 128
>>> import matplotlib.pyplot as plt >>> clock = cp.arange(64, len(sig), 128) >>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True) >>> ax_orig.plot(cp.asnumpy(sig)) >>> ax_orig.plot(cp.asnumpy(clock), cp.asnumpy(sig[clock]), 'ro') >>> ax_orig.set_title('Original signal') >>> ax_noise.plot(cp.asnumpy(sig_noise)) >>> ax_noise.set_title('Signal with noise') >>> ax_corr.plot(cp.asnumpy(corr)) >>> ax_corr.plot(cp.asnumpy(clock), cp.asnumpy(corr[clock]), 'ro') >>> ax_corr.axhline(0.5, ls=':') >>> ax_corr.set_title('Cross-correlated with rectangular pulse') >>> ax_orig.margins(0, 0.1) >>> fig.tight_layout() >>> fig.show()
- cusignal.convolution.correlate.correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0)#
Cross-correlate two 2-dimensional arrays. Cross correlate in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue.
- Parameters:
- in1array_like
First input.
- in2array_like
Second input. Should have the same number of dimensions as in1.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output:
full
The output is the full discrete linear cross-correlation of the inputs. (Default)
valid
The output consists only of those elements that do not rely on the zero-padding. In ‘valid’ mode, either in1 or in2 must be at least as large as the other in every dimension.
same
The output is the same size as in1, centered with respect to the ‘full’ output.
- boundarystr {‘fill’, ‘wrap’, ‘symm’}, optional
A flag indicating how to handle boundaries:
fill
pad input arrays with fillvalue. (default)
wrap
circular boundary conditions.
symm
symmetrical boundary conditions.
- fillvaluescalar, optional
Value to fill pad input arrays with. Default is 0.
- Returns:
- correlate2dndarray
A 2-dimensional array containing a subset of the discrete linear cross-correlation of in1 with in2.
Examples
Use 2D cross-correlation to find the location of a template in a noisy image: >>> import cusignal >>> import cupy as cp >>> from scipy import misc >>> face = cp.asarray(misc.face(gray=True) - misc.face(gray=True).mean()) >>> template = cp.copy(face[300:365, 670:750]) # right eye >>> template -= template.mean() >>> face = face + cp.random.randn(*face.shape) * 50 # add noise >>> corr = cusignal.correlate2d(face, template, boundary=’symm’, mode=’same’) >>> y, x = cp.unravel_index(cp.argmax(corr), corr.shape) # find the match >>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_template, ax_corr) = … plt.subplots(3, 1, figsize=(6, 15)) >>> ax_orig.imshow(cp.asnumpy(face), cmap=’gray’) >>> ax_orig.set_title(‘Original’) >>> ax_orig.set_axis_off() >>> ax_template.imshow(cp.asnumpy(template), cmap=’gray’) >>> ax_template.set_title(‘Template’) >>> ax_template.set_axis_off() >>> ax_corr.imshow(cp.asnumpy(corr), cmap=’gray’) >>> ax_corr.set_title(‘Cross-correlation’) >>> ax_corr.set_axis_off() >>> ax_orig.plot(cp.asnumpy(x), cp.asnumpy(y), ‘ro’) >>> fig.show()
- cusignal.convolution.correlate.correlation_lags(in1_len, in2_len, mode='full')#
Calculates the lag / displacement indices array for 1D cross-correlation. Parameters ———- in1_size : int
First input size.
- in2_sizeint
Second input size.
- modestr {‘full’, ‘valid’, ‘same’}, optional
A string indicating the size of the output. See the documentation correlate for more information.
See Also#
correlate : Compute the N-dimensional cross-correlation. Returns ——- lags : array
Returns an array containing cross-correlation lag/displacement indices. Indices can be indexed with the np.argmax of the correlation to return the lag/displacement.
Notes#
Cross-correlation for continuous functions \(f\) and \(g\) is defined as: .. math:
\left ( f\star g \right )\left ( \tau \right ) \triangleq \int_{t_0}^{t_0 +T} \overline{f\left ( t \right )}g\left ( t+\tau \right )dt
Where \(\tau\) is defined as the displacement, also known as the lag. Cross correlation for discrete functions \(f\) and \(g\) is defined as: .. math:
\left ( f\star g \right )\left [ n \right ] \triangleq \sum_{-\infty}^{\infty} \overline{f\left [ m \right ]}g\left [ m+n \right ]
Where \(n\) is the lag. Examples ——– Cross-correlation of a signal with its time-delayed self. >>> import cusignal >>> import cupy as cp >>> from cupy.random import default_rng >>> rng = default_rng() >>> x = rng.standard_normal(1000) >>> y = cp.concatenate([rng.standard_normal(100), x]) >>> correlation = cusignal.correlate(x, y, mode=”full”) >>> lags = cusignal.correlation_lags(x.size, y.size, mode=”full”) >>> lag = lags[cp.argmax(correlation)]
Modulation/Demodulation#
Demodulation#
- cusignal.demod.demod.fm_demod(x, axis=-1)#
Demodulate Frequency Modulated Signal
- Parameters:
- xndarray
Received complex valued signal or batch of signals
- Returns:
- yndarray
The demodulated output with the same shape as x.
Estimation#
Kalman Filter#
- class cusignal.estimation.filters.KalmanFilter(dim_x, dim_z, dim_u=0, points=1, dtype=<class 'numpy.float32'>)#
This is a multi-point Kalman Filter implementation of rlabbe/filterpy, with a subset of functionality.
All Kalman Filter matrices are stack on the X axis. This is to allow for optimal global accesses on the GPU.
- Parameters:
- dim_xint
Number of state variables for the Kalman filter. For example, if you are tracking the position and velocity of an object in two dimensions, dim_x would be 4. This is used to set the default size of P, Q, and u
- dim_zint
Number of of measurement inputs. For example, if the sensor provides you with position in (x,y), dim_z would be 2.
- dim_uint (optional)
Size of the control input, if it is being used. Default value of 0 indicates it is not used.
- pointsint (optional)
Number of Kalman Filter points to track.
- dtypedtype (optional)
Data type of compute.
References
[1]Dan Simon. “Optimal State Estimation.” John Wiley & Sons. p. 208-212. (2006)
[2]Roger Labbe. “Kalman and Bayesian Filters in Python” rlabbe/Kalman-and-Bayesian-Filters-in-Python
Examples
Here is a filter that tracks position and velocity using a sensor that only reads position.
First construct the object with the required dimensionality, number of points, and data type.
import cupy as cp import numpy as np from cusignal import KalmanFilter points = 1024 kf = KalmanFilter(dim_x=4, dim_z=2, points=points, dtype=cp.float64)
Assign the initial value for the state (position and velocity) for all Kalman Filter points.
initial_location = np.array( [[10.0, 10.0, 0.0, 0.0]], dtype=dt ).T # x, y, v_x, v_y kf.x = cp.repeat( cp.asarray(initial_location[cp.newaxis, :, :]), points, axis=0 )
Define the state transition matrix for all Kalman Filter points:
F = np.array( [ [1.0, 0.0, 1.0, 0.0], # x = x0 + v_x*dt [0.0, 1.0, 0.0, 1.0], # y = y0 + v_y*dt [0.0, 0.0, 1.0, 0.0], # dx = v_x [1.0, 0.0, 0.0, 1.0], ], # dy = v_y dtype=dt, ) kf.F = cp.repeat(cp.asarray(F[cp.newaxis, :, :]), points, axis=0)
Define the measurement function for all Kalman Filter points:
H = np.array( [[1.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 1.0]], dtype=dt, # x_0 # y_0 ) kf.H = cp.repeat(cp.asarray(H[cp.newaxis, :, :]), points, axis=0)
Define the covariance matrix for all Kalman Filter points:
initial_estimate_error = np.eye(dim_x, dtype=dt) * np.array( [1.0, 1.0, 2.0, 2.0], dtype=dt ) kf.P = cp.repeat( cp.asarray(initial_estimate_error[cp.newaxis, :, :]), points, axis=0, )
Define the measurement noise for all Kalman Filter points:
measurement_noise = np.eye(dim_z, dtype=dt) * 0.01 kf.R = cp.repeat( cp.asarray(measurement_noise[cp.newaxis, :, :]), points, axis=0 )
Define the process noise for all Kalman Filter points:
Now just perform the standard predict/update loop: Note: This example just uses the same sensor reading for all points
kf.predict() z = get_sensor_reading() (dim_z, 1) kf.z = cp.repeat(z[cp.newaxis, :, :], points, axis=0) kf.update()
Results are in:
Attributes
x
(array(points, dim_x, 1)) Current state estimate. Any call to update() or predict() updates this variable.
P
(array(points, dim_x, dim_x)) Current state covariance matrix. Any call to update() or predict() updates this variable.
z
(array(points, dim_z, 1)) Last measurement used in update(). Read only.
R
(array(points, dim_z, dim_z)) Measurement noise matrix
Q
(array(points, dim_x, dim_x)) Process noise matrix
F
(array(points, dim_x, dim_x)) State Transition matrix
H
(array(points, dim_z, dim_x)) Measurement function
_alpha_sq
(float (points, 1, 1)) Fading memory setting. 1.0 gives the normal Kalman filter, and values slightly larger than 1.0 (such as 1.02) give a fading memory effect - previous measurements have less influence on the filter’s estimates. This formulation of the Fading memory filter (there are many) is due to Dan Simon [1].
Methods
predict
([u, B, F, Q])Predict next state (prior) using the Kalman filter state propagation equations.
update
(z[, R, H])Add a new measurement (z) to the Kalman filter.
- predict(u=None, B=None, F=None, Q=None)#
Predict next state (prior) using the Kalman filter state propagation equations.
- Parameters:
- unarray, default 0
Optional control vector.
- Barray(points, dim_x, dim_u), or None
Optional control transition matrix; a value of None will cause the filter to use self.B.
- Farray(points, dim_x, dim_x), or None
Optional state transition matrix; a value of None will cause the filter to use self.F.
- Qarray(points, dim_x, dim_x), scalar, or None
Optional process noise matrix; a value of None will cause the filter to use self.Q.
- update(z, R=None, H=None)#
Add a new measurement (z) to the Kalman filter.
- Parameters:
- zarray(points, dim_z, 1)
measurement for this update. z can be a scalar if dim_z is 1, otherwise it must be convertible to a column vector. If you pass in a value of H, z must be a column vector the of the correct size.
- Rarray(points, dim_z, dim_z), scalar, or None
Optionally provide R to override the measurement noise for this one call, otherwise self.R will be used.
- Harray(points, dim_z, dim_x), or None
Optionally provide H to override the measurement function for this one call, otherwise self.H will be used.
Filtering#
Resample#
- cusignal.filtering.resample.decimate(x, q, n=None, axis=-1, zero_phase=True, gpupath=True)#
Downsample the signal after applying an anti-aliasing filter.
- Parameters:
- xarray_like
The signal to be downsampled, as an N-dimensional array.
- qint
The downsampling factor.
- nint or array_like, optional
The order of the filter (1 less than the length for FIR) to calculate, or the FIR filter coefficients to employ. Defaults to calculating the coefficients for 20 times the downsampling factor.
- axisint, optional
The axis along which to decimate.
- zero_phasebool, optional
Prevent shifting the outputs back by the filter’s group delay when using an FIR filter. The default value of
True
is recommended, since a phase shift is generally not desired.- gpupathbool, Optional
Optional path for filter design. gpupath == False may be desirable if filter sizes are small.
- Returns:
- yndarray
The down-sampled signal.
See also
resample
Resample up or down using the FFT method.
resample_poly
Resample using polyphase filtering and an FIR filter.
Notes
Only FIR filter types are currently supported in cuSignal.
- cusignal.filtering.resample.resample(x, num, t=None, axis=0, window=None, domain='time')#
Resample x to num samples using Fourier method along the given axis.
The resampled signal starts at the same value as x but is sampled with a spacing of
len(x) / num * (spacing of x)
. Because a Fourier method is used, the signal is assumed to be periodic.- Parameters:
- xarray_like
The data to be resampled.
- numint
The number of samples in the resampled signal.
- tarray_like, optional
If t is given, it is assumed to be the sample positions associated with the signal data in x.
- axisint, optional
The axis of x that is resampled. Default is 0.
- windowarray_like, callable, string, float, or tuple, optional
Specifies the window applied to the signal in the Fourier domain. See below for details.
- domainstring, optional
A string indicating the domain of the input x:
time
Consider the input x as time-domain. (Default)
freq
Consider the input x as frequency-domain.
- Returns:
- resampled_x or (resampled_x, resampled_t)
Either the resampled array, or, if t was given, a tuple containing the resampled array and the corresponding resampled positions.
See also
decimate
Downsample the signal after applying an FIR or IIR filter.
resample_poly
Resample using polyphase filtering and an FIR filter.
Notes
The argument window controls a Fourier-domain window that tapers the Fourier spectrum before zero-padding to alleviate ringing in the resampled values for sampled signals you didn’t intend to be interpreted as band-limited.
If window is a function, then it is called with a vector of inputs indicating the frequency bins (i.e. fftfreq(x.shape[axis]) ).
If window is an array of the same length as x.shape[axis] it is assumed to be the window to be applied directly in the Fourier domain (with dc and low-frequency first).
For any other type of window, the function cusignal.get_window is called to generate the window.
The first sample of the returned vector is the same as the first sample of the input vector. The spacing between samples is changed from
dx
todx * len(x) / num
.If t is not None, then it represents the old sample positions, and the new sample positions will be returned as well as the new samples.
As noted, resample uses FFT transformations, which can be very slow if the number of input or output samples is large and prime; see scipy.fftpack.fft.
Examples
Note that the end of the resampled data rises to meet the first sample of the next cycle:
>>> import cusignal >>> import cupy as cp
>>> x = cp.linspace(0, 10, 20, endpoint=False) >>> y = cp.cos(-x**2/6.0) >>> f = cusignal.resample(y, 100) >>> xnew = cp.linspace(0, 10, 100, endpoint=False)
>>> import matplotlib.pyplot as plt >>> plt.plot(cp.asnumpy(x), cp.asnumpy(y), 'go-', cp.asnumpy(xnew), cp.asnumpy(f), '.-', 10, cp.asnumpy(y[0]), 'ro') >>> plt.legend(['data', 'resampled'], loc='best') >>> plt.show()
- cusignal.filtering.resample.resample_poly(x, up, down, axis=0, window=('kaiser', 5.0), gpupath=True)#
Resample x along the given axis using polyphase filtering.
The signal x is upsampled by the factor up, a zero-phase low-pass FIR filter is applied, and then it is downsampled by the factor down. The resulting sample rate is
up / down
times the original sample rate. Values beyond the boundary of the signal are assumed to be zero during the filtering step.- Parameters:
- xarray_like
The data to be resampled.
- upint
The upsampling factor.
- downint
The downsampling factor.
- axisint, optional
The axis of x that is resampled. Default is 0.
- windowstring, tuple, or array_like, optional
Desired window to use to design the low-pass filter, or the FIR filter coefficients to employ. See below for details.
- gpupathbool, Optional
Optional path for filter design. gpupath == False may be desirable if filter sizes are small.
- Returns:
- resampled_xarray
The resampled array.
See also
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
todx * down / float(up)
.Examples
Note that the end of the resampled data rises to meet the first sample of the next cycle for the FFT method, and gets closer to zero for the polyphase method:
>>> import cusignal >>> import cupy as cp
>>> x = cp.linspace(0, 10, 20, endpoint=False) >>> y = cp.cos(-x**2/6.0) >>> f_fft = cusignal.resample(y, 100) >>> f_poly = cusignal.resample_poly(y, 100, 20) >>> xnew = cp.linspace(0, 10, 100, endpoint=False)
>>> import matplotlib.pyplot as plt >>> plt.plot(cp.asnumpy(xnew), cp.asnumpy(f_fft), 'b.-', cp.asnumpy(xnew), cp.asnumpy(f_poly), 'r.-') >>> plt.plot(cp.asnumpy(x), cp.asnumpy(y), 'ko-') >>> plt.plot(10, cp.asnumpy(y[0]), 'bo', 10, 0., 'ro') # boundaries >>> plt.legend(['resample', 'resamp_poly', 'data'], loc='best') >>> plt.show()
- cusignal.filtering.resample.upfirdn(h, x, up=1, down=1, axis=-1)#
Upsample, FIR filter, and downsample.
- Parameters:
- harray_like
1-dimensional FIR (finite-impulse response) filter coefficients.
- xarray_like
Input signal array.
- upint, optional
Upsampling rate. Default is 1.
- downint, optional
Downsampling rate. Default is 1.
- axisint, optional
The axis of the input data array along which to apply the linear filter. The filter is applied to each subarray along this axis. Default is -1.
- Returns:
- yndarray
The output signal array. Dimensions will be the same as x except for along axis, which will change size according to the h, up, and down parameters.
Notes
The algorithm is an implementation of the block diagram shown on page 129 of the Vaidyanathan text [1] (Figure 4.3-8d).
The direct approach of upsampling by factor of P with zero insertion, FIR filtering of length
N
, and downsampling by factor of Q is O(N*Q) per output sample. The polyphase implementation used here is O(N/P).References
[1]P. P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice Hall, 1993.
Examples
Simple operations: >>> from cusignal import upfirdn >>> upfirdn([1, 1, 1], [1, 1, 1]) # FIR filter array([ 1., 2., 3., 2., 1.]) >>> upfirdn([1], [1, 2, 3], 3) # upsampling with zeros insertion array([ 1., 0., 0., 2., 0., 0., 3., 0., 0.]) >>> upfirdn([1, 1, 1], [1, 2, 3], 3) # upsampling with sample-and-hold array([ 1., 1., 1., 2., 2., 2., 3., 3., 3.]) >>> upfirdn([.5, 1, .5], [1, 1, 1], 2) # linear interpolation array([ 0.5, 1. , 1. , 1. , 1. , 1. , 0.5, 0. ]) >>> upfirdn([1], cp.arange(10), 1, 3) # decimation by 3 array([ 0., 3., 6., 9.]) >>> upfirdn([.5, 1, .5], cp.arange(10), 2, 3) # linear interp, rate 2/3 array([ 0. , 1. , 2.5, 4. , 5.5, 7. , 8.5, 0. ]) Apply a single filter to multiple signals: >>> x = cp.reshape(cp.arange(8), (4, 2)) >>> x array([[0, 1],
[2, 3], [4, 5], [6, 7]])
Apply along the last dimension of
x
: >>> h = [1, 1] >>> upfirdn(h, x, 2) array([[ 0., 0., 1., 1.],[ 2., 2., 3., 3.], [ 4., 4., 5., 5.], [ 6., 6., 7., 7.]])
Apply along the 0th dimension of
x
: >>> upfirdn(h, x, 2, axis=0) array([[ 0., 1.],[ 0., 1.], [ 2., 3.], [ 2., 3.], [ 4., 5.], [ 4., 5.], [ 6., 7.], [ 6., 7.]])
FIR Filters#
Filter data along one-dimension with an FIR filter.
Filter a data sequence, x, using a digital filter. This works for many fundamental data types (including Object type). Please note, cuSignal doesn’t support IIR filters presently, and this implementation is optimized for large filtering operations (and inherently depends on fftconvolve)
Parameters#
- barray_like
The numerator coefficient vector in a 1-D sequence.
- xarray_like
An N-dimensional input array.
- axisint, optional
The axis of the input data array along which to apply the linear filter. The filter is applied to each subarray along this axis. Default is -1.
- ziarray_like, optional
Initial conditions for the filter delays. It is a vector (or array of vectors for an N-dimensional input) of length
max(len(a), len(b)) - 1
. If zi is None or is not given then initial rest is assumed. See lfiltic for more information.
Returns#
- yarray
The output of the digital filter.
- zfarray, optional
If zi is None, this is not returned, otherwise, zf holds the final filter delay values.
- cusignal.filter_design.fir_filter_design.cmplx_sort(p)#
Sort roots based on magnitude.
- Parameters:
- parray_like
The roots to sort, as a 1-D array.
- Returns:
- p_sortedndarray
Sorted roots.
- indxndarray
Array of indices needed to sort the input p.
Examples
>>> import cusignal >>> vals = [1, 4, 1+1.j, 3] >>> p_sorted, indx = cusignal.cmplx_sort(vals) >>> p_sorted array([1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j]) >>> indx array([0, 2, 3, 1])
- cusignal.filter_design.fir_filter_design.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None, gpupath=True)#
FIR filter design using the window method.
This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if numtaps is odd and Type II if numtaps is even.
Type II filters always have zero response at the Nyquist frequency, so a ValueError exception is raised if firwin is called with numtaps even and having a passband whose right end is at the Nyquist frequency.
- Parameters:
- numtapsint
Length of the filter (number of coefficients, i.e. the filter order + 1). numtaps must be odd if a passband includes the Nyquist frequency.
- cutofffloat or 1D array_like
Cutoff frequency of filter (expressed in the same units as fs) OR an array of cutoff frequencies (that is, band edges). In the latter case, the frequencies in cutoff should be positive and monotonically increasing between 0 and fs/2. The values 0 and fs/2 must not be included in cutoff.
- widthfloat or None, optional
If width is not None, then assume it is the approximate width of the transition region (expressed in the same units as fs) for use in Kaiser FIR filter design. In this case, the window argument is ignored.
- windowstring or tuple of string and parameter values, optional
Desired window to use. See cusignal.get_window for a list of windows and required parameters.
- pass_zero{True, False, ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’},
optional If True, the gain at the frequency 0 (i.e. the “DC gain”) is 1. If False, the DC gain is 0. Can also be a string argument for the desired filter type (equivalent to
btype
in IIR design functions).New in version 1.3.0: Support for string arguments.
- scalebool, optional
Set to True to scale the coefficients so that the frequency response is exactly unity at a certain frequency. That frequency is either:
0 (DC) if the first passband starts at 0 (i.e. pass_zero is True)
fs/2 (the Nyquist frequency) if the first passband ends at fs/2 (i.e the filter is a single band highpass filter); center of first passband otherwise
- nyqfloat, optional
Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in cutoff must be between 0 and nyq. Default is 1.
- fsfloat, optional
The sampling frequency of the signal. Each frequency in cutoff must be between 0 and
fs/2
. Default is 2.- gpupathbool, Optional
Optional path for filter design. gpupath == False may be desirable if filter sizes are small.
- Returns:
- h(numtaps,) ndarray
Coefficients of length numtaps FIR filter.
- Raises:
- ValueError
If any value in cutoff is less than or equal to 0 or greater than or equal to
fs/2
, if the values in cutoff are not strictly monotonically increasing, or if numtaps is even but a passband includes the Nyquist frequency.
See also
firwin2
firls
minimum_phase
remez
Examples
Low-pass from 0 to f:
>>> import cusignal >>> numtaps = 3 >>> f = 0.1 >>> cusignal.firwin(numtaps, f) array([ 0.06799017, 0.86401967, 0.06799017])
Use a specific window function:
>>> cusignal.firwin(numtaps, f, window='nuttall') array([ 3.56607041e-04, 9.99286786e-01, 3.56607041e-04])
High-pass (‘stop’ from 0 to f):
>>> cusignal.firwin(numtaps, f, pass_zero=False) array([-0.00859313, 0.98281375, -0.00859313])
Band-pass:
>>> f1, f2 = 0.1, 0.2 >>> cusignal.firwin(numtaps, [f1, f2], pass_zero=False) array([ 0.06301614, 0.88770441, 0.06301614])
Band-stop:
>>> cusignal.firwin(numtaps, [f1, f2]) array([-0.00801395, 1.0160279 , -0.00801395])
Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1]):
>>> f3, f4 = 0.3, 0.4 >>> cusignal.firwin(numtaps, [f1, f2, f3, f4]) array([-0.01376344, 1.02752689, -0.01376344])
Multi-band (passbands are [f1, f2] and [f3,f4]):
>>> cusignal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False) array([ 0.04890915, 0.91284326, 0.04890915])
- cusignal.filter_design.fir_filter_design.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None, antisymmetric=False, fs=None, gpupath=True)#
FIR filter design using the window method. From the given frequencies freq and corresponding gains gain, this function constructs an FIR filter with linear phase and (approximately) the given frequency response. Parameters ———- numtaps : int
The number of taps in the FIR filter. numtaps must be less than nfreqs.
- freqarray_like, 1-D
The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being Nyquist. The Nyquist frequency is half fs. The values in freq must be nondecreasing. A value can be repeated once to implement a discontinuity. The first value in freq must be 0, and the last value must be
fs/2
. Values 0 andfs/2
must not be repeated.- gainarray_like
The filter gains at the frequency sampling points. Certain constraints to gain values, depending on the filter type, are applied, see Notes for details.
- nfreqsint, optional
The size of the interpolation mesh used to construct the filter. For most efficient behavior, this should be a power of 2 plus 1 (e.g, 129, 257, etc). The default is one more than the smallest power of 2 that is not less than numtaps. nfreqs must be greater than numtaps.
- windowstring or (string, float) or float, or None, optional
Window function to use. Default is “hamming”. See scipy.signal.get_window for the complete list of possible values. If None, no window function is applied.
- nyqfloat, optional
Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in freq must be between 0 and nyq. Default is 1.
- antisymmetricbool, optional
Whether resulting impulse response is symmetric/antisymmetric. See Notes for more details.
- fsfloat, optional
The sampling frequency of the signal. Each frequency in cutoff must be between 0 and
fs/2
. Default is 2.
Returns#
- tapsndarray
The filter coefficients of the FIR filter, as a 1-D array of length numtaps.
See also#
firls firwin minimum_phase remez Notes —– From the given set of frequencies and gains, the desired response is constructed in the frequency domain. The inverse FFT is applied to the desired response to create the associated convolution kernel, and the first numtaps coefficients of this kernel, scaled by window, are returned. The FIR filter will have linear phase. The type of filter is determined by the value of ‘numtaps` and antisymmetric flag. There are four possible combinations:
odd numtaps, antisymmetric is False, type I filter is produced
even numtaps, antisymmetric is False, type II filter is produced
odd numtaps, antisymmetric is True, type III filter is produced
even numtaps, antisymmetric is True, type IV filter is produced
Magnitude response of all but type I filters are subjects to following constraints:
type II – zero at the Nyquist frequency
type III – zero at zero and Nyquist frequencies
type IV – zero at zero frequency
New in version 0.9.0.
References#
[1]Oppenheim, A. V. and Schafer, R. W., “Discrete-Time Signal Processing”, Prentice-Hall, Englewood Cliffs, New Jersey (1989). (See, for example, Section 7.4.)
[2]Smith, Steven W., “The Scientist and Engineer’s Guide to Digital Signal Processing”, Ch. 17. http://www.dspguide.com/ch17/1.htm
- cusignal.filter_design.fir_filter_design.kaiser_atten(numtaps, width)#
Compute the attenuation of a Kaiser FIR filter. Given the number of taps N and the transition width width, compute the attenuation a in dB, given by Kaiser’s formula:
a = 2.285 * (N - 1) * pi * width + 7.95
Parameters#
- numtapsint
The number of taps in the FIR filter.
- widthfloat
The desired width of the transition region between passband and stopband (or, in general, at any discontinuity) for the filter.
Returns#
- afloat
The attenuation of the ripple, in dB.
- cusignal.filter_design.fir_filter_design.kaiser_beta(a)#
Compute the Kaiser parameter beta, given the attenuation a. Parameters ———- a : float
The desired attenuation in the stopband and maximum ripple in the passband, in dB. This should be a positive number.
Returns#
- betafloat
The beta parameter to be used in the formula for a Kaiser window.
References#
Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.
Channelizer#
Polyphase channelize signal into n channels
Parameters#
- xarray_like
The input data to be channelized
- harray_like
The 1-D input filter; will be split into n channels of int number of taps
- n_chansint
Number of channels for channelizer
Returns#
yy : channelized output matrix
Notes#
Currently only supports simple channelizer where channel spacing is equivalent to the number of channels used (zero overlap). Number of filter taps (len of filter / n_chans) must be <=32.
Filter Design#
Resample#
- cusignal.filter_design.fir_filter_design.cmplx_sort(p)#
Sort roots based on magnitude.
- Parameters:
- parray_like
The roots to sort, as a 1-D array.
- Returns:
- p_sortedndarray
Sorted roots.
- indxndarray
Array of indices needed to sort the input p.
Examples
>>> import cusignal >>> vals = [1, 4, 1+1.j, 3] >>> p_sorted, indx = cusignal.cmplx_sort(vals) >>> p_sorted array([1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j]) >>> indx array([0, 2, 3, 1])
- cusignal.filter_design.fir_filter_design.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None, gpupath=True)#
FIR filter design using the window method.
This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if numtaps is odd and Type II if numtaps is even.
Type II filters always have zero response at the Nyquist frequency, so a ValueError exception is raised if firwin is called with numtaps even and having a passband whose right end is at the Nyquist frequency.
- Parameters:
- numtapsint
Length of the filter (number of coefficients, i.e. the filter order + 1). numtaps must be odd if a passband includes the Nyquist frequency.
- cutofffloat or 1D array_like
Cutoff frequency of filter (expressed in the same units as fs) OR an array of cutoff frequencies (that is, band edges). In the latter case, the frequencies in cutoff should be positive and monotonically increasing between 0 and fs/2. The values 0 and fs/2 must not be included in cutoff.
- widthfloat or None, optional
If width is not None, then assume it is the approximate width of the transition region (expressed in the same units as fs) for use in Kaiser FIR filter design. In this case, the window argument is ignored.
- windowstring or tuple of string and parameter values, optional
Desired window to use. See cusignal.get_window for a list of windows and required parameters.
- pass_zero{True, False, ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’},
optional If True, the gain at the frequency 0 (i.e. the “DC gain”) is 1. If False, the DC gain is 0. Can also be a string argument for the desired filter type (equivalent to
btype
in IIR design functions).New in version 1.3.0: Support for string arguments.
- scalebool, optional
Set to True to scale the coefficients so that the frequency response is exactly unity at a certain frequency. That frequency is either:
0 (DC) if the first passband starts at 0 (i.e. pass_zero is True)
fs/2 (the Nyquist frequency) if the first passband ends at fs/2 (i.e the filter is a single band highpass filter); center of first passband otherwise
- nyqfloat, optional
Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in cutoff must be between 0 and nyq. Default is 1.
- fsfloat, optional
The sampling frequency of the signal. Each frequency in cutoff must be between 0 and
fs/2
. Default is 2.- gpupathbool, Optional
Optional path for filter design. gpupath == False may be desirable if filter sizes are small.
- Returns:
- h(numtaps,) ndarray
Coefficients of length numtaps FIR filter.
- Raises:
- ValueError
If any value in cutoff is less than or equal to 0 or greater than or equal to
fs/2
, if the values in cutoff are not strictly monotonically increasing, or if numtaps is even but a passband includes the Nyquist frequency.
See also
firwin2
firls
minimum_phase
remez
Examples
Low-pass from 0 to f:
>>> import cusignal >>> numtaps = 3 >>> f = 0.1 >>> cusignal.firwin(numtaps, f) array([ 0.06799017, 0.86401967, 0.06799017])
Use a specific window function:
>>> cusignal.firwin(numtaps, f, window='nuttall') array([ 3.56607041e-04, 9.99286786e-01, 3.56607041e-04])
High-pass (‘stop’ from 0 to f):
>>> cusignal.firwin(numtaps, f, pass_zero=False) array([-0.00859313, 0.98281375, -0.00859313])
Band-pass:
>>> f1, f2 = 0.1, 0.2 >>> cusignal.firwin(numtaps, [f1, f2], pass_zero=False) array([ 0.06301614, 0.88770441, 0.06301614])
Band-stop:
>>> cusignal.firwin(numtaps, [f1, f2]) array([-0.00801395, 1.0160279 , -0.00801395])
Multi-band (passbands are [0, f1], [f2, f3] and [f4, 1]):
>>> f3, f4 = 0.3, 0.4 >>> cusignal.firwin(numtaps, [f1, f2, f3, f4]) array([-0.01376344, 1.02752689, -0.01376344])
Multi-band (passbands are [f1, f2] and [f3,f4]):
>>> cusignal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False) array([ 0.04890915, 0.91284326, 0.04890915])
- cusignal.filter_design.fir_filter_design.firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None, antisymmetric=False, fs=None, gpupath=True)#
FIR filter design using the window method. From the given frequencies freq and corresponding gains gain, this function constructs an FIR filter with linear phase and (approximately) the given frequency response. Parameters ———- numtaps : int
The number of taps in the FIR filter. numtaps must be less than nfreqs.
- freqarray_like, 1-D
The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being Nyquist. The Nyquist frequency is half fs. The values in freq must be nondecreasing. A value can be repeated once to implement a discontinuity. The first value in freq must be 0, and the last value must be
fs/2
. Values 0 andfs/2
must not be repeated.- gainarray_like
The filter gains at the frequency sampling points. Certain constraints to gain values, depending on the filter type, are applied, see Notes for details.
- nfreqsint, optional
The size of the interpolation mesh used to construct the filter. For most efficient behavior, this should be a power of 2 plus 1 (e.g, 129, 257, etc). The default is one more than the smallest power of 2 that is not less than numtaps. nfreqs must be greater than numtaps.
- windowstring or (string, float) or float, or None, optional
Window function to use. Default is “hamming”. See scipy.signal.get_window for the complete list of possible values. If None, no window function is applied.
- nyqfloat, optional
Deprecated. Use `fs` instead. This is the Nyquist frequency. Each frequency in freq must be between 0 and nyq. Default is 1.
- antisymmetricbool, optional
Whether resulting impulse response is symmetric/antisymmetric. See Notes for more details.
- fsfloat, optional
The sampling frequency of the signal. Each frequency in cutoff must be between 0 and
fs/2
. Default is 2.
Returns#
- tapsndarray
The filter coefficients of the FIR filter, as a 1-D array of length numtaps.
See also#
firls firwin minimum_phase remez Notes —– From the given set of frequencies and gains, the desired response is constructed in the frequency domain. The inverse FFT is applied to the desired response to create the associated convolution kernel, and the first numtaps coefficients of this kernel, scaled by window, are returned. The FIR filter will have linear phase. The type of filter is determined by the value of ‘numtaps` and antisymmetric flag. There are four possible combinations:
odd numtaps, antisymmetric is False, type I filter is produced
even numtaps, antisymmetric is False, type II filter is produced
odd numtaps, antisymmetric is True, type III filter is produced
even numtaps, antisymmetric is True, type IV filter is produced
Magnitude response of all but type I filters are subjects to following constraints:
type II – zero at the Nyquist frequency
type III – zero at zero and Nyquist frequencies
type IV – zero at zero frequency
New in version 0.9.0.
References#
[1]Oppenheim, A. V. and Schafer, R. W., “Discrete-Time Signal Processing”, Prentice-Hall, Englewood Cliffs, New Jersey (1989). (See, for example, Section 7.4.)
[2]Smith, Steven W., “The Scientist and Engineer’s Guide to Digital Signal Processing”, Ch. 17. http://www.dspguide.com/ch17/1.htm
- cusignal.filter_design.fir_filter_design.kaiser_atten(numtaps, width)#
Compute the attenuation of a Kaiser FIR filter. Given the number of taps N and the transition width width, compute the attenuation a in dB, given by Kaiser’s formula:
a = 2.285 * (N - 1) * pi * width + 7.95
Parameters#
- numtapsint
The number of taps in the FIR filter.
- widthfloat
The desired width of the transition region between passband and stopband (or, in general, at any discontinuity) for the filter.
Returns#
- afloat
The attenuation of the ripple, in dB.
- cusignal.filter_design.fir_filter_design.kaiser_beta(a)#
Compute the Kaiser parameter beta, given the attenuation a. Parameters ———- a : float
The desired attenuation in the stopband and maximum ripple in the passband, in dB. This should be a positive number.
Returns#
- betafloat
The beta parameter to be used in the formula for a Kaiser window.
References#
Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.
Peak Finding#
Peak Finding#
- cusignal.peak_finding.peak_finding.argrelextrema(data, comparator, axis=0, order=1, mode='clip')#
Calculate the relative extrema of data.
- Parameters:
- datandarray
Array in which to find the relative extrema.
- comparatorcallable
Function to use to compare two data points. Should take two arrays as arguments.
- axisint, optional
Axis over which to select from data. Default is 0.
- orderint, optional
How many points on each side to use for the comparison to consider
comparator(n, n+x)
to be True.
- Returns:
- extrematuple of ndarrays
Indices of the maxima in arrays of integers.
extrema[k]
is the array of indices of axis k of data. Note that the return value is a tuple even when data is one-dimensional.
Examples
>>> from cusignal import argrelextrema >>> import cupy as cp >>> x = cp.array([2, 1, 2, 3, 2, 0, 1, 0]) >>> argrelextrema(x, cp.greater) (array([3, 6]),) >>> y = cp.array([[1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]]) ... >>> argrelextrema(y, cp.less, axis=1) (array([0, 2]), array([2, 1]))
- cusignal.peak_finding.peak_finding.argrelmax(data, axis=0, order=1, mode='clip')#
Calculate the relative maxima of data.
- Parameters:
- datandarray
Array in which to find the relative maxima.
- axisint, optional
Axis over which to select from data. Default is 0.
- orderint, optional
How many points on each side to use for the comparison to consider
comparator(n, n+x)
to be True.
- Returns:
- extrematuple of ndarrays
Indices of the maxima in arrays of integers.
extrema[k]
is the array of indices of axis k of data. Note that the return value is a tuple even when data is one-dimensional.
See also
argrelextrema
,argrelmin
,find_peaks
Notes
This function uses argrelextrema with cp.greater as comparator. Therefore it requires a strict inequality on both sides of a value to consider it a maximum. This means flat maxima (more than one sample wide) are not detected. In case of one-dimensional data find_peaks can be used to detect all local maxima, including flat ones.
Examples
>>> from cusignal import argrelmax >>> import cupy as cp >>> x = cp.array([2, 1, 2, 3, 2, 0, 1, 0]) >>> argrelmax(x) (array([3, 6]),) >>> y = cp.array([[1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]]) ... >>> argrelmax(y, axis=1) (array([0]), array([1]))
- cusignal.peak_finding.peak_finding.argrelmin(data, axis=0, order=1, mode='clip')#
Calculate the relative minima of data.
- Parameters:
- datandarray
Array in which to find the relative minima.
- axisint, optional
Axis over which to select from data. Default is 0.
- orderint, optional
How many points on each side to use for the comparison to consider
comparator(n, n+x)
to be True.
- Returns:
- extrematuple of ndarrays
Indices of the minima in arrays of integers.
extrema[k]
is the array of indices of axis k of data. Note that the return value is a tuple even when data is one-dimensional.
See also
argrelextrema
,argrelmax
,find_peaks
Notes
This function uses argrelextrema with cp.less as comparator. Therefore it requires a strict inequality on both sides of a value to consider it a minimum. This means flat minima (more than one sample wide) are not detected. In case of one-dimensional data find_peaks can be used to detect all local minima, including flat ones, by calling it with negated data.
Examples
>>> from cusignal import argrelmin >>> import cupy as cp >>> x = cp.array([2, 1, 2, 3, 2, 0, 1, 0]) >>> argrelmin(x) (array([1, 5]),) >>> y = cp.array([[1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]]) ... >>> argrelmin(y, axis=1) (array([0, 2]), array([2, 1]))
Window Functions#
Windows#
- cusignal.windows.windows.barthann(M, sym=True)#
Return a modified Bartlett-Hann window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.barthann(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Bartlett-Hann window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Bartlett-Hann window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.bartlett(M, sym=True)#
Return a Bartlett window.
The Bartlett window is very similar to a triangular window, except that the end points are at zero. It is often used in signal processing for tapering a signal, without generating too much ripple in the frequency domain.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The triangular window, with the first and last samples equal to zero and the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
See also
triang
A triangular window that does not touch zero at the ends
Notes
The Bartlett window is defined as
\[w(n) = \frac{2}{M-1} \left( \frac{M-1}{2} - \left|n - \frac{M-1}{2}\right| \right)\]Most references to the Bartlett window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. Note that convolution with this window produces linear interpolation. It is also known as an apodization (which means”removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function. The Fourier transform of the Bartlett is the product of two sinc functions. Note the excellent discussion in Kanasewich. [2]
References
[1]M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika 37, 1-16, 1950.
[2]E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 109-110.
[3]A.V. Oppenheim and R.W. Schafer, “Discrete-Time Signal Processing”, Prentice-Hall, 1999, pp. 468-471.
[4]Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function
[5]W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, “Numerical Recipes”, Cambridge University Press, 1986, page 429.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.bartlett(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Bartlett window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Bartlett window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.blackman(M, sym=True)#
Return a Blackman window.
The Blackman window is a taper formed by using the first three terms of a summation of cosines. It was designed to have close to the minimal leakage possible. It is close to optimal, only slightly worse than a Kaiser window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The Blackman window is defined as
\[w(n) = 0.42 - 0.5 \cos(2\pi n/M) + 0.08 \cos(4\pi n/M)\]The “exact Blackman” window was designed to null out the third and fourth sidelobes, but has discontinuities at the boundaries, resulting in a 6 dB/oct fall-off. This window is an approximation of the “exact” window, which does not null the sidelobes as well, but is smooth at the edges, improving the fall-off rate to 18 dB/oct. [3]
Most references to the Blackman window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function. It is known as a “near optimal” tapering function, almost as good (by some measures) as the Kaiser window.
References
[1]Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.
[2]Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
[3]Harris, Fredric J. (Jan 1978). “On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform”. Proceedings of the IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.blackman(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Blackman window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = cp.abs(fftshift(A / cp.abs(A).max())) >>> response = 20 * cp.log10(cp.maximum(response, 1e-10)) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Blackman window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.blackmanharris(M, sym=True)#
Return a minimum 4-term Blackman-Harris window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.blackmanharris(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Blackman-Harris window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Blackman-Harris window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.bohman(M, sym=True)#
Return a Bohman window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.bohman(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Bohman window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Bohman window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.boxcar(M, sym=True)#
Return a boxcar or rectangular window.
Also known as a rectangular window or Dirichlet window, this is equivalent to no window at all.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
Whether the window is symmetric. (Has no effect for boxcar.)
- Returns:
- wndarray
The window, with the maximum value normalized to 1.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.boxcar(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Boxcar window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the boxcar window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.chebwin(M, at, sym=True)#
Return a Dolph-Chebyshev window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- atfloat
Attenuation (in dB).
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value always normalized to 1
Notes
This window optimizes for the narrowest main lobe width for a given order M and sidelobe equiripple attenuation at, using Chebyshev polynomials. It was originally developed by Dolph to optimize the directionality of radio antenna arrays.
Unlike most windows, the Dolph-Chebyshev is defined in terms of its frequency response:
\[W(k) = \frac {\cos\{M \cos^{-1}[\beta \cos(\frac{\pi k}{M})]\}} {\cosh[M \cosh^{-1}(\beta)]}\]where
\[\beta = \cosh \left [\frac{1}{M} \cosh^{-1}(10^\frac{A}{20}) \right ]\]and 0 <= abs(k) <= M-1. A is the attenuation in decibels (at).
The time domain window is then generated using the IFFT, so power-of-two M are the fastest to generate, and prime number M are the slowest.
The equiripple condition in the frequency domain creates impulses in the time domain, which appear at the ends of the window.
References
[1]C. Dolph, “A current distribution for broadside arrays which optimizes the relationship between beam width and side-lobe level”, Proceedings of the IEEE, Vol. 34, Issue 6
[2]Peter Lynch, “The Dolph-Chebyshev Window: A Simple Optimal Filter”, American Meteorological Society (April 1997) http://mathsci.ucd.ie/~plynch/Publications/Dolph.pdf
[3]F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transforms”, Proceedings of the IEEE, Vol. 66, No. 1, January 1978
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.chebwin(51, at=100) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Dolph-Chebyshev window (100 dB)") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Dolph-Chebyshev window (100 dB)") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.cosine(M, sym=True)#
Return a window with a simple cosine shape.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
New in version 0.13.0.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.cosine(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Cosine window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the cosine window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]") >>> plt.show()
- cusignal.windows.windows.exponential(M, center=None, tau=1.0, sym=True)#
Return an exponential (or Poisson) window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- centerfloat, optional
Parameter defining the center location of the window function. The default value if not given is
center = (M-1) / 2
. This parameter must take its default value for symmetric windows.- taufloat, optional
Parameter defining the decay. For
center = 0
usetau = -(M-1) / ln(x)
ifx
is the fraction of the window remaining at the end.- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The Exponential window is defined as
\[w(n) = e^{-|n-center| / \tau}\]References
S. Gade and H. Herlufsen, “Windows to FFT analysis (Part I)”, Technical Review 3, Bruel & Kjaer, 1987.
Examples
Plot the symmetric window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> M = 51 >>> tau = 3.0 >>> window = cusignal.exponential(M, tau=tau) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Exponential Window (tau=3.0)") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -35, 0]) >>> plt.title("Frequency response of the Exponential window (tau=3.0)") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
This function can also generate non-symmetric windows:
>>> tau2 = -(M-1) / np.log(0.01) >>> window2 = cusignal.exponential(M, 0, tau2, False) >>> plt.figure() >>> plt.plot(cp.asnumpy(window2)) >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
- cusignal.windows.windows.flattop(M, sym=True)#
Return a flat top window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
Flat top windows are used for taking accurate measurements of signal amplitude in the frequency domain, with minimal scalloping error from the center of a frequency bin to its edges, compared to others. This is a 5th-order cosine window, with the 5 terms optimized to make the main lobe maximally flat. [1]
References
[1]D’Antona, Gabriele, and A. Ferrero, “Digital Signal Processing for Measurement Systems”, Springer Media, 2006, p. 70 :doi:`10.1007/0-387-28666-7`.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.flattop(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Flat top window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the flat top window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.gaussian(M, std, sym=True)#
Return a Gaussian window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- stdfloat
The standard deviation, sigma.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The Gaussian window is defined as
\[w(n) = e^{ -\frac{1}{2}\left(\frac{n}{\sigma}\right)^2 }\]Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.gaussian(51, std=7) >>> plt.plot(cp.asnumpy(window)) >>> plt.title(r"Gaussian window ($\sigma$=7)") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title(r"Frequency response of the Gaussian window ($\sigma$=7)") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.general_cosine(M, a, sym=True)#
Generic weighted sum of cosine terms window
- Parameters:
- Mint
Number of points in the output window
- aarray_like
Sequence of weighting coefficients. This uses the convention of being centered on the origin, so these will typically all be positive numbers, not alternating sign.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
References
[1]A. Nuttall, “Some windows with very good sidelobe behavior,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.
[2]Heinzel G. et al., “Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows”, February 15, 2002 https://holometer.fnal.gov/GH_FFT.pdf
Examples
Heinzel describes a flat-top window named “HFT90D” with formula: [2]
\[w_j = 1 - 1.942604 \cos(z) + 1.340318 \cos(2z) - 0.440811 \cos(3z) + 0.043097 \cos(4z)\]where
\[z = \frac{2 \pi j}{N}, j = 0...N - 1\]Since this uses the convention of starting at the origin, to reproduce the window, we need to convert every other coefficient to a positive number:
>>> HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097]
The paper states that the highest sidelobe is at -90.2 dB. Reproduce Figure 42 by plotting the window and its frequency response, and confirm the sidelobe level in red:
>>> from cusignal.windows import general_cosine >>> from cupy.fft import fft, fftshift >>> import cupy as cp >>> import matplotlib.pyplot as plt
>>> window = general_cosine(1000, HFT90D, sym=False) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("HFT90D window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 10000) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = cp.abs(fftshift(A / cp.abs(A).max())) >>> response = 20 * cp.log10(cp.maximum(response, 1e-10)) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-50/1000, 50/1000, -140, 0]) >>> plt.title("Frequency response of the HFT90D window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]") >>> plt.axhline(-90.2, color='red') >>> plt.show()
- cusignal.windows.windows.general_gaussian(M, p, sig, sym=True)#
Return a window with a generalized Gaussian shape.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- pfloat
Shape parameter. p = 1 is identical to gaussian, p = 0.5 is the same shape as the Laplace distribution.
- sigfloat
The standard deviation, sigma.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The generalized Gaussian window is defined as
\[w(n) = e^{ -\frac{1}{2}\left|\frac{n}{\sigma}\right|^{2p} }\]the half-power point is at
\[(2 \log(2))^{1/(2 p)} \sigma\]Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.general_gaussian(51, p=1.5, sig=7) >>> plt.plot(cp.asnumpy(window)) >>> plt.title(r"Generalized Gaussian window (p=1.5, $\sigma$=7)") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title(r"Freq. resp. of the gen. Gaussian " ... r"window (p=1.5, $\sigma$=7)") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.general_hamming(M, alpha, sym=True)#
Return a generalized Hamming window.
The generalized Hamming window is constructed by multiplying a rectangular window by one period of a cosine function [1].
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- alphafloat
The window coefficient, \(\alpha\)
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The generalized Hamming window is defined as
\[w(n) = \alpha - \left(1 - \alpha\right) \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1\]Both the common Hamming window and Hann window are special cases of the generalized Hamming window with \(\alpha\) = 0.54 and \(\alpha\) = 0.5, respectively [2].
References
[1]DSPRelated, “Generalized Hamming Window Family”, https://www.dsprelated.com/freebooks/sasp/Generalized_Hamming_Window_Family.html
[2]Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function
[3]Riccardo Piantanida ESA, “Sentinel-1 Level 1 Detailed Algorithm Definition”, https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Level-1-Detailed-Algorithm-Definition
[4]Matthieu Bourbigot ESA, “Sentinel-1 Product Definition”, https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Product-Definition
Examples
The Sentinel-1A/B Instrument Processing Facility uses generalized Hamming windows in the processing of spaceborne Synthetic Aperture Radar (SAR) data [3]. The facility uses various values for the \(\alpha\) parameter based on operating mode of the SAR instrument. Some common \(\alpha\) values include 0.75, 0.7 and 0.52 [4]. As an example, we plot these different windows.
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> fig1, spatial_plot = plt.subplots() >>> spatial_plot.set_title("Generalized Hamming Windows") >>> spatial_plot.set_ylabel("Amplitude") >>> spatial_plot.set_xlabel("Sample")
>>> fig2, freq_plot = plt.subplots() >>> freq_plot.set_title("Frequency Responses") >>> freq_plot.set_ylabel("Normalized magnitude [dB]") >>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]")
>>> for alpha in [0.75, 0.7, 0.52]: ... window = cusignal.general_hamming(41, alpha) ... spatial_plot.plot(cp.asnumpy(window), label="{:.2f}".format(alpha)) ... A = fft(window, 2048) / (len(window)/2.0) ... freq = cp.linspace(-0.5, 0.5, len(A)) ... response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) ... freq_plot.plot( ... cp.asnumpy(freq), cp.asnumpy(response), ... label="{:.2f}".format(alpha) ... ) >>> freq_plot.legend(loc="upper right") >>> spatial_plot.legend(loc="upper right")
- cusignal.windows.windows.get_window(window, Nx, fftbins=True)#
Return a window of a given length and type.
- Parameters:
- windowstring, float, or tuple
The type of window to create. See below for more details.
- Nxint
The number of samples in the window.
- fftbinsbool, optional
If True (default), create a “periodic” window, ready to use with ifftshift and be multiplied by the result of an FFT (see also fftpack.fftfreq). If False, create a “symmetric” window, for use in filter design.
- Returns:
- get_windowndarray
Returns a window of length Nx and type window
Notes
Window types:
~cusignal.windows.windows.boxcar
~cusignal.windows.windows.triang
~cusignal.windows.windows.blackman
~cusignal.windows.windows.hamming
~cusignal.windows.windows.hann
~cusignal.windows.windows.bartlett
~cusignal.windows.windows.flattop
~cusignal.windows.windows.parzen
~cusignal.windows.windows.bohman
~cusignal.windows.windows.blackmanharris
~cusignal.windows.windows.nuttall
~cusignal.windows.windows.barthann
~cusignal.windows.windows.kaiser (needs beta)
~cusignal.windows.windows.gaussian (needs standard deviation)
- ~cusignal.windows.windows.general_gaussian
(needs power, width)
~cusignal.windows.windows.slepian (needs width)
- ~cusignal.windows.windows.dpss
(needs normalized half-bandwidth)
~cusignal.windows.windows.chebwin (needs attenuation)
~cusignal.windows.windows.exponential (needs decay scale)
~cusignal.windows.windows.tukey (needs taper fraction)
If the window requires no parameters, then window can be a string.
If the window requires parameters, then window must be a tuple with the first argument the string name of the window, and the next arguments the needed parameters.
If window is a floating point number, it is interpreted as the beta parameter of the ~cusignal.windows.windows.kaiser window.
Each of the window types listed above is also the name of a function that can be called directly to create a window of that type.
Examples
>>> import cusignal >>> cusignal.get_window('triang', 7) array([ 0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375]) >>> cusignal.get_window(('kaiser', 4.0), 9) array([0.08848053, 0.32578323, 0.63343178, 0.89640418, 1., 0.89640418, 0.63343178, 0.32578323, 0.08848053]) >>> cusignal.get_window(4.0, 9) array([0.08848053, 0.32578323, 0.63343178, 0.89640418, 1., 0.89640418, 0.63343178, 0.32578323, 0.08848053])
- cusignal.windows.windows.hamming(M, sym=True)#
Return a Hamming window.
The Hamming window is a taper formed by using a raised cosine with non-zero endpoints, optimized to minimize the nearest side lobe.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The Hamming window is defined as
\[w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1\]The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and is described in Blackman and Tukey. It was recommended for smoothing the truncated autocovariance function in the time domain. Most references to the Hamming window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function.
References
[1]Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.
[2]E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 109-110.
[3]Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function
[4]W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, “Numerical Recipes”, Cambridge University Press, 1986, page 425.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.hamming(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Hamming window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Hamming window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.hann(M, sym=True)#
Return a Hann window.
The Hann window is a taper formed by using a raised cosine or sine-squared with ends that touch zero.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The Hann window is defined as
\[w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1\]The window was named for Julius von Hann, an Austrian meteorologist. It is also known as the Cosine Bell. It is sometimes erroneously referred to as the “Hanning” window, from the use of “hann” as a verb in the original paper and confusion with the very similar Hamming window.
Most references to the Hann window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function.
References
[1]Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.
[2]E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 106-108.
[3]Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function
[4]W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, “Numerical Recipes”, Cambridge University Press, 1986, page 425.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.hann(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Hann window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = cp.abs(fftshift(A / cp.abs(A).max())) >>> response = 20 * cp.log10(np.maximum(response, 1e-10)) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Hann window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.kaiser(M, beta, sym=True)#
Return a Kaiser window.
The Kaiser window is a taper formed by using a Bessel function.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- betafloat
Shape parameter, determines trade-off between main-lobe width and side lobe level. As beta gets large, the window narrows.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
Notes
The Kaiser window is defined as
\[w(n) = I_0\left( \beta \sqrt{1-\frac{4n^2}{(M-1)^2}} \right)/I_0(\beta)\]with
\[\quad -\frac{M-1}{2} \leq n \leq \frac{M-1}{2},\]where \(I_0\) is the modified zeroth-order Bessel function.
The Kaiser was named for Jim Kaiser, who discovered a simple approximation to the DPSS window based on Bessel functions. The Kaiser window is a very good approximation to the Digital Prolate Spheroidal Sequence, or Slepian window, which is the transform which maximizes the energy in the main lobe of the window relative to total energy.
The Kaiser can approximate other windows by varying the beta parameter. (Some literature uses alpha = beta/pi.) [4]
beta
Window shape
0
Rectangular
5
Similar to a Hamming
6
Similar to a Hann
8.6
Similar to a Blackman
A beta value of 14 is probably a good starting point. Note that as beta gets large, the window narrows, and so the number of samples needs to be large enough to sample the increasingly narrow spike, otherwise NaNs will be returned.
Most references to the Kaiser window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means “removing the foot”, i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function.
References
[1]J. F. Kaiser, “Digital Filters” - Ch 7 in “Systems analysis by digital computer”, Editors: F.F. Kuo and J.F. Kaiser, p 218-285. John Wiley and Sons, New York, (1966).
[2]E.R. Kanasewich, “Time Sequence Analysis in Geophysics”, The University of Alberta Press, 1975, pp. 177-178.
[3]Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function
[4]F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978. :doi:`10.1109/PROC.1978.10837`.
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.kaiser(51, beta=14) >>> plt.plot(cp.asnumpy(window)) >>> plt.title(r"Kaiser window ($\beta$=14)") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title(r"Frequency response of the Kaiser window ($\beta$=14)") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.nuttall(M, sym=True)#
Return a minimum 4-term Blackman-Harris window according to Nuttall.
This variation is called “Nuttall4c” by Heinzel. [2]
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
References
[1]A. Nuttall, “Some windows with very good sidelobe behavior,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.
[2]Heinzel G. et al., “Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows”, February 15, 2002 https://holometer.fnal.gov/GH_FFT.pdf
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.nuttall(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Nuttall window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Nuttall window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.parzen(M, sym=True)#
- cusignal.windows.windows.taylor(M, nbar=4, sll=30, norm=True, sym=True)#
Return a Taylor window. The Taylor window taper function approximates the Dolph-Chebyshev window’s constant sidelobe level for a parameterized number of near-in sidelobes, but then allows a taper beyond [2]. The SAR (synthetic aperature radar) community commonly uses Taylor weighting for image formation processing because it provides strong, selectable sidelobe suppression with minimum broadening of the mainlobe [Rc5e098a1a171-1]. Parameters ———- M : int
Number of points in the output window. If zero or less, an empty array is returned.
- nbarint, optional
Number of nearly constant level sidelobes adjacent to the mainlobe.
- sllfloat, optional
Desired suppression of sidelobe level in decibels (dB) relative to the DC gain of the mainlobe. This should be a positive number.
- normbool, optional
When True (default), divides the window by the largest (middle) value for odd-length windows or the value that would occur between the two repeated middle values for even-length windows such that all values are less than or equal to 1. When False the DC gain will remain at 1 (0 dB) and the sidelobes will be sll dB down.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
Returns#
- outarray
The window. When norm is True (default), the maximum value is normalized to 1 (though the value 1 does not appear if M is even and sym is True).
See Also#
chebwin, kaiser, bartlett, blackman, hamming, hanning References ———- .. [Rc5e098a1a171-1] W. Carrara, R. Goodman, and R. Majewski, “Spotlight Synthetic
Aperture Radar: Signal Processing Algorithms” Pages 512-513, July 1995.
[2]Armin Doerry, “Catalog of Window Taper Functions for Sidelobe Control”, 2017. https://www.researchgate.net/profile/Armin_Doerry/publication/316281181_Catalog_of_Window_Taper_Functions_for_Sidelobe_Control/links/58f92cb2a6fdccb121c9d54d/Catalog-of-Window-Taper-Functions-for-Sidelobe-Control.pdf
Examples#
Plot the window and its frequency response: >>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt >>> window = signal.windows.taylor(51, nbar=20, sll=100, norm=False) >>> plt.plot(window) >>> plt.title(“Taylor window (100 dB)”) >>> plt.ylabel(“Amplitude”) >>> plt.xlabel(“Sample”) >>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title(“Frequency response of the Taylor window (100 dB)”) >>> plt.ylabel(“Normalized magnitude [dB]”) >>> plt.xlabel(“Normalized frequency [cycles per sample]”)
- cusignal.windows.windows.triang(M, sym=True)#
Return a triangular window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
See also
bartlett
A triangular window that touches zero
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.triang(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Triangular window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample")
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = cp.abs(fftshift(A / cp.abs(A).max())) >>> response = 20 * cp.log10(cp.maximum(response, 1e-10)) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the triangular window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
- cusignal.windows.windows.tukey(M, alpha=0.5, sym=True)#
Return a Tukey window, also known as a tapered cosine window.
- Parameters:
- Mint
Number of points in the output window. If zero or less, an empty array is returned.
- alphafloat, optional
Shape parameter of the Tukey window, representing the fraction of the window inside the cosine tapered region. If zero, the Tukey window is equivalent to a rectangular window. If one, the Tukey window is equivalent to a Hann window.
- symbool, optional
When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.
- Returns:
- wndarray
The window, with the maximum value normalized to 1 (though the value 1 does not appear if M is even and sym is True).
References
[1]Harris, Fredric J. (Jan 1978). “On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform”. Proceedings of the IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`
[2]Wikipedia, “Window function”, https://en.wikipedia.org/wiki/Window_function#Tukey_window
Examples
Plot the window and its frequency response:
>>> import cusignal >>> import cupy as cp >>> from cupy.fft import fft, fftshift >>> import matplotlib.pyplot as plt
>>> window = cusignal.tukey(51) >>> plt.plot(cp.asnumpy(window)) >>> plt.title("Tukey window") >>> plt.ylabel("Amplitude") >>> plt.xlabel("Sample") >>> plt.ylim([0, 1.1])
>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = cp.linspace(-0.5, 0.5, len(A)) >>> response = 20 * cp.log10(cp.abs(fftshift(A / cp.abs(A).max()))) >>> plt.plot(cp.asnumpy(freq), cp.asnumpy(response)) >>> plt.axis([-0.5, 0.5, -120, 0]) >>> plt.title("Frequency response of the Tukey window") >>> plt.ylabel("Normalized magnitude [dB]") >>> plt.xlabel("Normalized frequency [cycles per sample]")
Waveforms#
Waveform Generation#
- cusignal.waveforms.waveforms.chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True, type='real')#
Frequency-swept cosine generator.
In the following, ‘Hz’ should be interpreted as ‘cycles per unit’; there is no requirement here that the unit is one second. The important distinction is that the units of rotation are cycles, not radians. Likewise, t could be a measurement of space instead of time.
- Parameters:
- tarray_like
Times at which to evaluate the waveform.
- f0float
Frequency (e.g. Hz) at time t=0.
- t1float
Time at which f1 is specified.
- f1float
Frequency (e.g. Hz) of the waveform at time t1.
- method{‘linear’, ‘quadratic’, ‘logarithmic’, ‘hyperbolic’}, optional
Kind of frequency sweep. If not given, linear is assumed. See Notes below for more details.
- phifloat, optional
Phase offset, in degrees. Default is 0.
- vertex_zerobool, optional
This parameter is only used when method is ‘quadratic’. It determines whether the vertex of the parabola that is the graph of the frequency is at t=0 or t=t1.
- type{‘real’, ‘complex’}, optional
Specify output chirp type, only applicable when method is ‘linear’.
- Returns:
- yndarray
A numpy array containing the signal evaluated at t with the requested time-varying frequency. More precisely, the function returns
cos(phase + (pi/180)*phi)
where phase is the integral (from 0 to t) of2*pi*f(t)
.f(t)
is defined below.
Examples
The following will be used in the examples:
>>> from cusignal import chirp, spectrogram >>> import matplotlib.pyplot as plt >>> import cupy as cp
For the first example, we’ll plot the waveform for a linear chirp from 6 Hz to 1 Hz over 10 seconds:
>>> t = cp.linspace(0, 10, 5001) >>> w = chirp(t, f0=6, f1=1, t1=10, method='linear') >>> plt.plot(cp.asnumpy(t), cp.asnumpy(w)) >>> plt.title("Linear Chirp, f(0)=6, f(10)=1") >>> plt.xlabel('t (sec)') >>> plt.show()
For the remaining examples, we’ll use higher frequency ranges, and demonstrate the result using cusignal.spectrogram. We’ll use a 10 second interval sampled at 8000 Hz.
>>> fs = 8000 >>> T = 10 >>> t = cp.linspace(0, T, T*fs, endpoint=False)
Quadratic chirp from 1500 Hz to 250 Hz over 10 seconds (vertex of the parabolic curve of the frequency is at t=0):
>>> w = chirp(t, f0=1500, f1=250, t1=10, method='quadratic') >>> ff, tt, Sxx = spectrogram(w, fs=fs, noverlap=256, nperseg=512, ... nfft=2048) >>> plt.pcolormesh(cp.asnumpy(tt), cp.asnumpy(ff[:513]), cp.asnumpy(Sxx[:513]), cmap='gray_r') >>> plt.title('Quadratic Chirp, f(0)=1500, f(10)=250') >>> plt.xlabel('t (sec)') >>> plt.ylabel('Frequency (Hz)') >>> plt.grid() >>> plt.show()
- cusignal.waveforms.waveforms.gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False, retenv=False)#
Return a Gaussian modulated sinusoid:
exp(-a t^2) exp(1j*2*pi*fc*t).
If retquad is True, then return the real and imaginary parts (in-phase and quadrature). If retenv is True, then return the envelope (unmodulated signal). Otherwise, return the real part of the modulated sinusoid.
- Parameters:
- tndarray or the string ‘cutoff’
Input array.
- fcint, optional
Center frequency (e.g. Hz). Default is 1000.
- bwfloat, optional
Fractional bandwidth in frequency domain of pulse (e.g. Hz). Default is 0.5.
- bwrfloat, optional
Reference level at which fractional bandwidth is calculated (dB). Default is -6.
- tprfloat, optional
If t is ‘cutoff’, then the function returns the cutoff time for when the pulse amplitude falls below tpr (in dB). Default is -60.
- retquadbool, optional
If True, return the quadrature (imaginary) as well as the real part of the signal. Default is False.
- retenvbool, optional
If True, return the envelope of the signal. Default is False.
- Returns:
- yIndarray
Real part of signal. Always returned.
- yQndarray
Imaginary part of signal. Only returned if retquad is True.
- yenvndarray
Envelope of signal. Only returned if retenv is True.
See also
cusignal.morlet
Examples
Plot real component, imaginary component, and envelope for a 5 Hz pulse, sampled at 100 Hz for 2 seconds:
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt >>> t = cp.linspace(-1, 1, 2 * 100, endpoint=False) >>> i, q, e = cusignal.gausspulse(t, fc=5, retquad=True, retenv=True) >>> plt.plot(cp.asnumpy(t), cp.asnumpy(i), cp.asnumpy(t), cp.asnumpy(q), cp.asnumpy(t), cp.asnumpy(e), '--')
- cusignal.waveforms.waveforms.sawtooth(t, width=1.0)#
Return a periodic sawtooth or triangle waveform. The sawtooth waveform has a period
2*pi
, rises from -1 to 1 on the interval 0 towidth*2*pi
, then drops from 1 to -1 on the intervalwidth*2*pi
to2*pi
. width must be in the interval [0, 1]. Note that this is not band-limited. It produces an infinite number of harmonics, which are aliased back and forth across the frequency spectrum. Parameters ———- t : array_likeTime.
- widtharray_like, optional
Width of the rising ramp as a proportion of the total cycle. Default is 1, producing a rising ramp, while 0 produces a falling ramp. width = 0.5 produces a triangle wave. If an array, causes wave shape to change over time, and must be the same length as t.
Returns#
- yndarray
Output array containing the sawtooth waveform.
Examples#
A 5 Hz waveform sampled at 500 Hz for 1 second: >>> from scipy import signal >>> import matplotlib.pyplot as plt >>> t = np.linspace(0, 1, 500) >>> plt.plot(t, signal.sawtooth(2 * np.pi * 5 * t))
- cusignal.waveforms.waveforms.square(t, duty=0.5)#
Return a periodic square-wave waveform.
The square wave has a period
2*pi
, has value +1 from 0 to2*pi*duty
and -1 from2*pi*duty
to2*pi
. duty must be in the interval [0,1].Note that this is not band-limited. It produces an infinite number of harmonics, which are aliased back and forth across the frequency spectrum.
- Parameters:
- tarray_like
The input time array.
- dutyarray_like, optional
Duty cycle. Default is 0.5 (50% duty cycle). If an array, causes wave shape to change over time, and must be the same length as t.
- Returns:
- yndarray
Output array containing the square waveform.
Examples
A 5 Hz waveform sampled at 500 Hz for 1 second:
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt >>> t = cp.linspace(0, 1, 500, endpoint=False) >>> plt.plot(cp.asnumpy(t), cp.asnumpy(cusignal.square(2 * cp.pi * 5 * t))) >>> plt.ylim(-2, 2)
A pulse-width modulated sine wave:
>>> plt.figure() >>> sig = cp.sin(2 * cp.pi * t) >>> pwm = cusignal.square(2 * cp.pi * 30 * t, duty=(sig + 1)/2) >>> plt.subplot(2, 1, 1) >>> plt.plot(cp.asnumpy(t), cp.asnumpy(sig)) >>> plt.subplot(2, 1, 2) >>> plt.plot(cp.asnumpy(t), cp.asnumpy(pwm)) >>> plt.ylim(-1.5, 1.5)
- cusignal.waveforms.waveforms.unit_impulse(shape, idx=None, dtype=<class 'float'>)#
Unit impulse signal (discrete delta function) or unit basis vector.
- Parameters:
- shapeint or tuple of int
Number of samples in the output (1-D), or a tuple that represents the shape of the output (N-D).
- idxNone or int or tuple of int or ‘mid’, optional
Index at which the value is 1. If None, defaults to the 0th element. If
idx='mid'
, the impulse will be centered atshape // 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 isnumpy.float64
.
- Returns:
- yndarray
Output array containing an impulse signal.
Notes
The 1D case is also known as the Kronecker delta.
Examples
An impulse at the 0th element (\(\delta[n]\)):
>>> import cusignal >>> import cupy as cp >>> cusignal.unit_impulse(8) array([ 1., 0., 0., 0., 0., 0., 0., 0.])
Impulse offset by 2 samples (\(\delta[n-2]\)):
>>> cusignal.unit_impulse(7, 2) array([ 0., 0., 1., 0., 0., 0., 0.])
2-dimensional impulse, centered:
>>> cusignal.unit_impulse((3, 3), 'mid') array([[ 0., 0., 0.], [ 0., 1., 0.], [ 0., 0., 0.]])
Impulse at (2, 2), using broadcasting:
>>> cusignal.unit_impulse((4, 4), 2) array([[ 0., 0., 0., 0.], [ 0., 0., 0., 0.], [ 0., 0., 1., 0.], [ 0., 0., 0., 0.]])
Spectrum Analysis#
Spectral#
- cusignal.spectral_analysis.spectral.coherence(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', axis=-1)#
Estimate the magnitude squared coherence estimate, Cxy, of discrete-time signals X and Y using Welch’s method.
Cxy = abs(Pxy)**2/(Pxx*Pyy)
, where Pxx and Pyy are power spectral density estimates of X and Y, and Pxy is the cross spectral density estimate of X and Y.- Parameters:
- xarray_like
Time series of measurement values
- yarray_like
Time series of measurement values
- fsfloat, optional
Sampling frequency of the x and y time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
- npersegint, optional
Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
- noverlap: int, optional
Number of points to overlap between segments. If None,
noverlap = nperseg // 2
. Defaults to None.- nfftint, optional
Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
- detrendstr or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
- axisint, optional
Axis along which the coherence is computed for both inputs; the default is over the last axis (i.e.
axis=-1
).
- Returns:
- fndarray
Array of sample frequencies.
- Cxyndarray
Magnitude squared coherence of x and y.
See also
periodogram
Simple, optionally modified periodogram
lombscargle
Lomb-Scargle periodogram for unevenly sampled data
welch
Power spectral density by Welch’s method.
csd
Cross spectral density by Welch’s method.
Notes
An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.
References
[1]P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967.
[2]Stoica, Petre, and Randolph Moses, “Spectral Analysis of Signals” Prentice Hall, 2005
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt
Generate two test signals with some common features.
>>> fs = 10e3 >>> N = 1e5 >>> amp = 20 >>> freq = 1234.0 >>> noise_power = 0.001 * fs / 2 >>> time = cp.arange(N) / fs >>> b, a = cusignal.butter(2, 0.25, 'low') >>> x = cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape) >>> # lfilter not implemented in cuSignal >>> y = cusignal.lfilter(b, a, x) >>> x += amp*cp.sin(2*cp.pi*freq*time) >>> y += cp.random.normal(scale=0.1*cp.sqrt(noise_power), size=time.shape)
Compute and plot the coherence.
>>> f, Cxy = cusignal.coherence(x, y, fs, nperseg=1024) >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Cxy)) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('Coherence') >>> plt.show()
- cusignal.spectral_analysis.spectral.csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')#
Estimate the cross power spectral density, Pxy, using Welch’s method.
- Parameters:
- xarray_like
Time series of measurement values
- yarray_like
Time series of measurement values
- fsfloat, optional
Sampling frequency of the x and y time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
- npersegint, optional
Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
- noverlap: int, optional
Number of points to overlap between segments. If None,
noverlap = nperseg // 2
. Defaults to None.- nfftint, optional
Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
- detrendstr or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
- return_onesidedbool, optional
If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
- scaling{ ‘density’, ‘spectrum’ }, optional
Selects between computing the cross spectral density (‘density’) where Pxy has units of V**2/Hz and computing the cross spectrum (‘spectrum’) where Pxy has units of V**2, if x and y are measured in V and fs is measured in Hz. Defaults to ‘density’
- axisint, optional
Axis along which the CSD is computed for both inputs; the default is over the last axis (i.e.
axis=-1
).- average{ ‘mean’, ‘median’ }, optional
Method to use when averaging periodograms. Defaults to ‘mean’.
- Returns:
- fndarray
Array of sample frequencies.
- Pxyndarray
Cross spectral density or cross power spectrum of x,y.
See also
periodogram
Simple, optionally modified periodogram
lombscargle
Lomb-Scargle periodogram for unevenly sampled data
welch
Power spectral density by Welch’s method. [Equivalent to csd(x,x)]
coherence
Magnitude squared coherence by Welch’s method.
Notes
By convention, Pxy is computed with the conjugate FFT of X multiplied by the FFT of Y.
If the input series differ in length, the shorter series will be zero-padded to match.
An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.
References
[1]P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967.
[2]Rabiner, Lawrence R., and B. Gold. “Theory and Application of Digital Signal Processing” Prentice-Hall, pp. 414-419, 1975
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt
Generate two test signals with some common features.
>>> fs = 10e3 >>> N = 1e5 >>> amp = 20 >>> freq = 1234.0 >>> noise_power = 0.001 * fs / 2 >>> time = cp.arange(N) / fs >>> b, a = cusignal.butter(2, 0.25, 'low') >>> x = cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape) >>> # lfilter not currently implemented in cuSignal >>> y = signal.lfilter(b, a, x) >>> x += amp*cp.sin(2*cp.pi*freq*time) >>> y += cp.random.normal(scale=0.1*cp.sqrt(noise_power), size=time.shape)
Compute and plot the magnitude of the cross spectral density.
>>> f, Pxy = cusignal.csd(x, y, fs, nperseg=1024) >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.abs(Pxy))) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('CSD [V**2/Hz]') >>> plt.show()
- cusignal.spectral_analysis.spectral.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2)#
Perform the inverse Short Time Fourier transform (iSTFT). Parameters ———- Zxx : array_like
STFT of the signal to be reconstructed. If a purely real array is passed, it will be cast to a complex data type.
- fsfloat, optional
Sampling frequency of the time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window. Must match the window used to generate the STFT for faithful inversion.
- npersegint, optional
Number of data points corresponding to each STFT segment. This parameter must be specified if the number of data points per segment is odd, or if the STFT was padded via
nfft > nperseg
. If None, the value depends on the shape of Zxx and input_onesided. If input_onesided is True,nperseg=2*(Zxx.shape[freq_axis] - 1)
. Otherwise,nperseg=Zxx.shape[freq_axis]
. Defaults to None.- noverlapint, optional
Number of points to overlap between segments. If None, half of the segment length. Defaults to None. When specified, the COLA constraint must be met (see Notes below), and should match the parameter used to generate the STFT. Defaults to None.
- nfftint, optional
Number of FFT points corresponding to each STFT segment. This parameter must be specified if the STFT was padded via
nfft > nperseg
. If None, the default values are the same as for nperseg, detailed above, with one exception: if input_onesided is True andnperseg==2*Zxx.shape[freq_axis] - 1
, nfft also takes on that value. This case allows the proper inversion of an odd-length unpadded STFT usingnfft=None
. Defaults to None.- input_onesidedbool, optional
If True, interpret the input array as one-sided FFTs, such as is returned by stft with
return_onesided=True
and numpy.fft.rfft. If False, interpret the input as a a two-sided FFT. Defaults to True.- boundarybool, optional
Specifies whether the input signal was extended at its boundaries by supplying a non-None
boundary
argument to stft. Defaults to True.- time_axisint, optional
Where the time segments of the STFT is located; the default is the last axis (i.e.
axis=-1
).- freq_axisint, optional
Where the frequency axis of the STFT is located; the default is the penultimate axis (i.e.
axis=-2
).
Returns#
- tndarray
Array of output data times.
- xndarray
iSTFT of Zxx.
See Also#
stft: Short Time Fourier Transform check_COLA: Check whether the Constant OverLap Add (COLA) constraint
is met
check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met Notes —– In order to enable inversion of an STFT via the inverse STFT with istft, the signal windowing must obey the constraint of “nonzero overlap add” (NOLA): .. math:: sum_{t}w^{2}[n-tH] ne 0 This ensures that the normalization factors that appear in the denominator of the overlap-add reconstruction equation .. math:: x[n]=frac{sum_{t}x_{t}[n]w[n-tH]}{sum_{t}w^{2}[n-tH]} are not zero. The NOLA constraint can be checked with the check_NOLA function. An STFT which has been modified (via masking or otherwise) is not guaranteed to correspond to a exactly realizible signal. This function implements the iSTFT via the least-squares estimation algorithm detailed in [2], which produces a signal that minimizes the mean squared error between the STFT of the returned signal and the modified STFT. .. versionadded:: 0.19.0 References ———- .. [Rb890c6d0cb06-1] Oppenheim, Alan V., Ronald W. Schafer, John R. Buck
“Discrete-Time Signal Processing”, Prentice Hall, 1999.
[2]Daniel W. Griffin, Jae S. Lim “Signal Estimation from Modified Short-Time Fourier Transform”, IEEE 1984, 10.1109/TASSP.1984.1164317
Examples#
>>> from scipy import signal >>> import matplotlib.pyplot as plt Generate a test signal, a 2 Vrms sine wave at 50Hz corrupted by 0.001 V**2/Hz of white noise sampled at 1024 Hz. >>> fs = 1024 >>> N = 10*fs >>> nperseg = 512 >>> amp = 2 * np.sqrt(2) >>> noise_power = 0.001 * fs / 2 >>> time = cp.arange(N) / float(fs) >>> carrier = amp * cp.sin(2*cp.pi*50*time) >>> noise = cp.random.normal(scale=cp.sqrt(noise_power), ... size=time.shape) >>> x = carrier + noise Compute the STFT, and plot its magnitude >>> f, t, Zxx = cusignal.stft(x, fs=fs, nperseg=nperseg) >>> f = cp.asnumpy(f) >>> t = cp.asnumpy(t) >>> Zxx = cp.asnumpy(Zxx) >>> plt.figure() >>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud') >>> plt.ylim([f[1], f[-1]]) >>> plt.title('STFT Magnitude') >>> plt.ylabel('Frequency [Hz]') >>> plt.xlabel('Time [sec]') >>> plt.yscale('log') >>> plt.show() Zero the components that are 10% or less of the carrier magnitude, then convert back to a time series via inverse STFT >>> Zxx = cp.where(cp.abs(Zxx) >= amp/10, Zxx, 0) >>> _, xrec = cusignal.istft(Zxx, fs) >>> xrec = cp.asnumpy(xrec) >>> x = cp.asnumpy(x) >>> time = cp.asnumpy(time) >>> carrier = cp.asnumpy(carrier) Compare the cleaned signal with the original and true carrier signals. >>> plt.figure() >>> plt.plot(time, x, time, xrec, time, carrier) >>> plt.xlim([2, 2.1])*+ >>> plt.xlabel('Time [sec]') >>> plt.ylabel('Signal') >>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier']) >>> plt.show() Note that the cleaned signal does not start as abruptly as the original, since some of the coefficients of the transient were also removed: >>> plt.figure() >>> plt.plot(time, x, time, xrec, time, carrier) >>> plt.xlim([0, 0.1]) >>> plt.xlabel('Time [sec]') >>> plt.ylabel('Signal') >>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier']) >>> plt.show()
- cusignal.spectral_analysis.spectral.lombscargle(x, y, freqs)#
Computes the Lomb-Scargle periodogram. The Lomb-Scargle periodogram was developed by Lomb [1] and further extended by Scargle [2] to find, and test the significance of weak periodic signals with uneven temporal sampling. When normalize is False (default) the computed periodogram is unnormalized, it takes the value
(A**2) * N/4
for a harmonic signal with amplitude A for sufficiently large N. When normalize is True the computed periodogram is normalized by the residuals of the data around a constant reference model (at zero). Input arrays should be one-dimensional and will be cast to float64.- Parameters:
- xarray_like
Sample times.
- yarray_like
Measurement values.
- freqsarray_like
Angular frequencies for output periodogram.
- precenterbool, optional
Pre-center amplitudes by subtracting the mean.
- normalizebool, optional
Compute normalized periodogram.
- Returns:
- pgramarray_like
Lomb-Scargle periodogram.
- Raises:
- ValueError
If the input arrays x and y do not have the same shape.
See also
istft
Inverse Short Time Fourier Transform
check_COLA
Check whether the Constant OverLap Add (COLA) constraint is met
welch
Power spectral density by Welch’s method
spectrogram
Spectrogram by Welch’s method
csd
Cross spectral density by Welch’s method
Notes
This subroutine calculates the periodogram using a slightly modified algorithm due to Townsend [3] which allows the periodogram to be calculated using only a single pass through the input arrays for each frequency. The algorithm running time scales roughly as O(x * freqs) or O(N^2) for a large number of samples and frequencies.
References
[1]N.R. Lomb “Least-squares frequency analysis of unequally spaced data”, Astrophysics and Space Science, vol 39, pp. 447-462, 1976
[2]J.D. Scargle “Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data”, The Astrophysical Journal, vol 263, pp. 835-853, 1982
[3]R.H.D. Townsend, “Fast calculation of the Lomb-Scargle periodogram using graphics processing units.”, The Astrophysical Journal Supplement Series, vol 191, pp. 247-253, 2010
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt First define some input parameters for the signal: >>> A = 2. >>> w = 1. >>> phi = 0.5 * cp.pi >>> nin = 1000 >>> nout = 100000 >>> frac_points = 0.9 # Fraction of points to select Randomly select a fraction of an array with timesteps: >>> r = cp.random.rand(nin) >>> x = cp.linspace(0.01, 10*cp.pi, nin) >>> x = x[r >= frac_points] Plot a sine wave for the selected times: >>> y = A * cp.sin(w*x+phi) Define the array of frequencies for which to compute the periodogram: >>> f = cp.linspace(0.01, 10, nout) Calculate Lomb-Scargle periodogram: >>> pgram = cusignal.lombscargle(x, y, f, normalize=True) Now make a plot of the input data: >>> plt.subplot(2, 1, 1) >>> plt.plot(cp.asnumpy(x), cp.asnumpy(y), 'b+') Then plot the normalized periodogram: >>> plt.subplot(2, 1, 2) >>> plt.plot(cp.asnumpy(f), cp.asnumpy(pgram)) >>> plt.show()
- cusignal.spectral_analysis.spectral.periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1)#
Estimate power spectral density using a periodogram.
- Parameters:
- xarray_like
Time series of measurement values
- fsfloat, optional
Sampling frequency of the x time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to ‘boxcar’.
- nfftint, optional
Length of the FFT used. If None the length of x will be used.
- detrendstr or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
- return_onesidedbool, optional
If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
- scaling{ ‘density’, ‘spectrum’ }, optional
Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’
- axisint, optional
Axis along which the periodogram is computed; the default is over the last axis (i.e.
axis=-1
).
- Returns:
- fndarray
Array of sample frequencies.
- Pxxndarray
Power spectral density or power spectrum of x.
See also
welch
Estimate power spectral density using Welch’s method
lombscargle
Lomb-Scargle periodogram for unevenly sampled data
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt >>> cp.random.seed(1234)
Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.
>>> fs = 10e3 >>> N = 1e5 >>> amp = 2*cp.sqrt(2) >>> freq = 1234.0 >>> noise_power = 0.001 * fs / 2 >>> time = cp.arange(N) / fs >>> x = amp*cp.sin(2*cp.pi*freq*time) >>> x += cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)
Compute and plot the power spectral density.
>>> f, Pxx_den = cusignal.periodogram(x, fs) >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Pxx_den)) >>> plt.ylim([1e-7, 1e2]) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('PSD [V**2/Hz]') >>> plt.show()
If we average the last half of the spectral density, to exclude the peak, we can recover the noise power on the signal.
>>> cp.mean(Pxx_den[25000:]) 0.00099728892368242854
Now compute and plot the power spectrum.
>>> f, Pxx_spec = cusignal.periodogram(x, fs, 'flattop', scaling='spectrum') >>> plt.figure() >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.sqrt(Pxx_spec))) >>> plt.ylim([1e-4, 1e1]) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('Linear spectrum [V RMS]') >>> plt.show()
The peak height in the power spectrum is an estimate of the RMS amplitude.
>>> cp.sqrt(Pxx_spec.max()) 2.0077340678640727
- cusignal.spectral_analysis.spectral.spectrogram(x, fs=1.0, window=('tukey', 0.25), nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, mode='psd')#
Compute a spectrogram with consecutive Fourier transforms.
Spectrograms can be used as a way of visualizing the change of a nonstationary signal’s frequency content over time.
- Parameters:
- xarray_like
Time series of measurement values
- fsfloat, optional
Sampling frequency of the x time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Tukey window with shape parameter of 0.25.
- npersegint, optional
Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
- noverlapint, optional
Number of points to overlap between segments. If None,
noverlap = nperseg // 8
. Defaults to None.- nfftint, optional
Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
- detrendstr or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
- return_onesidedbool, optional
If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
- scaling{ ‘density’, ‘spectrum’ }, optional
Selects between computing the power spectral density (‘density’) where Sxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Sxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’.
- axisint, optional
Axis along which the spectrogram is computed; the default is over the last axis (i.e.
axis=-1
).- modestr, optional
Defines what kind of return values are expected. Options are [‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’]. ‘complex’ is equivalent to the output of stft with no padding or boundary extension. ‘magnitude’ returns the absolute magnitude of the STFT. ‘angle’ and ‘phase’ return the complex angle of the STFT, with and without unwrapping, respectively.
- Returns:
- fndarray
Array of sample frequencies.
- tndarray
Array of segment times.
- Sxxndarray
Spectrogram of x. By default, the last axis of Sxx corresponds to the segment times.
See also
periodogram
Simple, optionally modified periodogram
lombscargle
Lomb-Scargle periodogram for unevenly sampled data
welch
Power spectral density by Welch’s method.
csd
Cross spectral density by Welch’s method.
Notes
An appropriate amount of overlap will depend on the choice of window and on your requirements. In contrast to welch’s method, where the entire data stream is averaged over, one may wish to use a smaller overlap (or perhaps none at all) when computing a spectrogram, to maintain some statistical independence between individual segments. It is for this reason that the default window is a Tukey window with 1/8th of a window’s length overlap at each end.
New in version 0.16.0.
References
[1]Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”, Prentice Hall, 1999.
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt
Generate a test signal, a 2 Vrms sine wave whose frequency is slowly modulated around 3kHz, corrupted by white noise of exponentially decreasing magnitude sampled at 10 kHz.
>>> fs = 10e3 >>> N = 1e5 >>> amp = 2 * cp.sqrt(2) >>> noise_power = 0.01 * fs / 2 >>> time = cp.arange(N) / float(fs) >>> mod = 500*cp.cos(2*cp.pi*0.25*time) >>> carrier = amp * cp.sin(2*cp.pi*3e3*time + mod) >>> noise = cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape) >>> noise *= cp.exp(-time/5) >>> x = carrier + noise
Compute and plot the spectrogram.
>>> f, t, Sxx = cusignal.spectrogram(x, fs) >>> plt.pcolormesh(cp.asnumpy(t), cp.asnumpy(f), cp.asnumpy(Sxx)) >>> plt.ylabel('Frequency [Hz]') >>> plt.xlabel('Time [sec]') >>> plt.show()
Note, if using output that is not one sided, then use the following:
>>> f, t, Sxx = cusignal.spectrogram(x, fs, return_onesided=False) >>> plt.pcolormesh(cp.asnumpy(t), cp.fft.fftshift(f), cp.fft.fftshift(Sxx, axes=0)) >>> plt.ylabel('Frequency [Hz]') >>> plt.xlabel('Time [sec]') >>> plt.show()
- cusignal.spectral_analysis.spectral.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend=False, return_onesided=True, boundary='zeros', padded=True, axis=-1)#
Compute the Short Time Fourier Transform (STFT).
STFTs can be used as a way of quantifying the change of a nonstationary signal’s frequency and phase content over time.
- Parameters:
- xarray_like
Time series of measurement values
- fsfloat, optional
Sampling frequency of the x time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
- npersegint, optional
Length of each segment. Defaults to 256.
- noverlapint, optional
Number of points to overlap between segments. If None,
noverlap = nperseg // 2
. Defaults to None. When specified, the COLA constraint must be met (see Notes below).- nfftint, optional
Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
- detrendstr or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to False.
- return_onesidedbool, optional
If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
- boundarystr or None, optional
Specifies whether the input signal is extended at both ends, and how to generate the new values, in order to center the first windowed segment on the first input point. This has the benefit of enabling reconstruction of the first input point when the employed window function starts at zero. Valid options are
['even', 'odd', 'constant', 'zeros', None]
. Defaults to ‘zeros’, for zero padding extension. I.e.[1, 2, 3, 4]
is extended to[0, 1, 2, 3, 4, 0]
fornperseg=3
.- paddedbool, optional
Specifies whether the input signal is zero-padded at the end to make the signal fit exactly into an integer number of window segments, so that all of the signal is included in the output. Defaults to True. Padding occurs after boundary extension, if boundary is not None, and padded is True, as is the default.
- axisint, optional
Axis along which the STFT is computed; the default is over the last axis (i.e.
axis=-1
).
- Returns:
- fndarray
Array of sample frequencies.
- tndarray
Array of segment times.
- Zxxndarray
STFT of x. By default, the last axis of Zxx corresponds to the segment times.
See also
welch
Power spectral density by Welch’s method.
spectrogram
Spectrogram by Welch’s method.
csd
Cross spectral density by Welch’s method.
lombscargle
Lomb-Scargle periodogram for unevenly sampled data
Notes
In order to enable inversion of an STFT via the inverse STFT in istft, the signal windowing must obey the constraint of “Nonzero OverLap Add” (NOLA), and the input signal must have complete windowing coverage (i.e.
(x.shape[axis] - nperseg) % (nperseg-noverlap) == 0
). The padded argument may be used to accomplish this.Given a time-domain signal \(x[n]\), a window \(w[n]\), and a hop size \(H\) = nperseg - noverlap, the windowed frame at time index \(t\) is given by
\[x_{t}[n]=x[n]w[n-tH]\]The overlap-add (OLA) reconstruction equation is given by
\[x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]}\]The NOLA constraint ensures that every normalization term that appears in the denomimator of the OLA reconstruction equation is nonzero. Whether a choice of window, nperseg, and noverlap satisfy this constraint can be tested with check_NOLA.
References
[1]Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”, Prentice Hall, 1999.
[2]Daniel W. Griffin, Jae S. Lim “Signal Estimation from Modified Short-Time Fourier Transform”, IEEE 1984, 10.1109/TASSP.1984.1164317
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt
Generate a test signal, a 2 Vrms sine wave whose frequency is slowly modulated around 3kHz, corrupted by white noise of exponentially decreasing magnitude sampled at 10 kHz.
>>> fs = 10e3 >>> N = 1e5 >>> amp = 2 * cp.sqrt(2) >>> noise_power = 0.01 * fs / 2 >>> time = cp.arange(N) / float(fs) >>> mod = 500*cp.cos(2*cp.pi*0.25*time) >>> carrier = amp * cp.sin(2*cp.pi*3e3*time + mod) >>> noise = cp.random.normal(scale=cp.sqrt(noise_power), ... size=time.shape) >>> noise *= cp.exp(-time/5) >>> x = carrier + noise
Compute and plot the STFT’s magnitude.
>>> f, t, Zxx = cusignal.stft(x, fs, nperseg=1000) >>> plt.pcolormesh(cp.asnumpy(t), cp.asnumpy(f), cp.asnumpy(cp.abs(Zxx)), \ vmin=0, vmax=amp) >>> plt.title('STFT Magnitude') >>> plt.ylabel('Frequency [Hz]') >>> plt.xlabel('Time [sec]') >>> plt.show()
- cusignal.spectral_analysis.spectral.vectorstrength(events, period)#
Determine the vector strength of the events corresponding to the given period.
The vector strength is a measure of phase synchrony, how well the timing of the events is synchronized to a single period of a periodic signal.
If multiple periods are used, calculate the vector strength of each. This is called the “resonating vector strength”.
- Parameters:
- events1D array_like
An array of time points containing the timing of the events.
- periodfloat or array_like
The period of the signal that the events should synchronize to. The period is in the same units as events. It can also be an array of periods, in which case the outputs are arrays of the same length.
- Returns:
- strengthfloat or 1D array
The strength of the synchronization. 1.0 is perfect synchronization and 0.0 is no synchronization. If period is an array, this is also an array with each element containing the vector strength at the corresponding period.
- phasefloat or array
The phase that the events are most strongly synchronized to in radians. If period is an array, this is also an array with each element containing the phase for the corresponding period.
References
[1]van Hemmen, JL, Longtin, A, and Vollmayr, AN. Testing resonating vector strength: Auditory system, electric fish, and noise. Chaos 21, 047508 (2011); :doi:`10.1063/1.3670512`.
[2]van Hemmen, JL. Vector strength after Goldberg, Brown, and von Mises: biological and mathematical perspectives. Biol Cybern. 2013 Aug;107(4):385-96. :doi:`10.1007/s00422-013-0561-7`.
[3]van Hemmen, JL and Vollmayr, AN. Resonating vector strength: what happens when we vary the “probing” frequency while keeping the spike times fixed. Biol Cybern. 2013 Aug;107(4):491-94. :doi:`10.1007/s00422-013-0560-8`.
- cusignal.spectral_analysis.spectral.welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')#
Estimate power spectral density using Welch’s method.
Welch’s method [1] computes an estimate of the power spectral density by dividing the data into overlapping segments, computing a modified periodogram for each segment and averaging the periodograms.
- Parameters:
- xarray_like
Time series of measurement values
- fsfloat, optional
Sampling frequency of the x time series. Defaults to 1.0.
- windowstr or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even by default. See get_window for a list of windows and required parameters. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window.
- npersegint, optional
Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window.
- noverlapint, optional
Number of points to overlap between segments. If None,
noverlap = nperseg // 2
. Defaults to None.- nfftint, optional
Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
- detrendstr or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to ‘constant’.
- return_onesidedbool, optional
If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.
- scaling{ ‘density’, ‘spectrum’ }, optional
Selects between computing the power spectral density (‘density’) where Pxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’
- axisint, optional
Axis along which the periodogram is computed; the default is over the last axis (i.e.
axis=-1
).- average{ ‘mean’, ‘median’ }, optional
Method to use when averaging periodograms. Defaults to ‘mean’.
- Returns:
- fndarray
Array of sample frequencies.
- Pxxndarray
Power spectral density or power spectrum of x.
See also
periodogram
Simple, optionally modified periodogram
lombscargle
Lomb-Scargle periodogram for unevenly sampled data
Notes
An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.
If noverlap is 0, this method is equivalent to Bartlett’s method [2].
New in version 0.12.0.
References
[1]P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967.
[2]M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika, vol. 37, pp. 1-16, 1950.
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt >>> cp.random.seed(1234)
Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.
>>> fs = 10e3 >>> N = 1e5 >>> amp = 2*cp.sqrt(2) >>> freq = 1234.0 >>> noise_power = 0.001 * fs / 2 >>> time = cp.arange(N) / fs >>> x = amp*cp.sin(2*cp.pi*freq*time) >>> x += cp.random.normal(scale=cp.sqrt(noise_power), size=time.shape)
Compute and plot the power spectral density.
>>> f, Pxx_den = cusignal.welch(x, fs, nperseg=1024) >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Pxx_den)) >>> plt.ylim([0.5e-3, 1]) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('PSD [V**2/Hz]') >>> plt.show()
If we average the last half of the spectral density, to exclude the peak, we can recover the noise power on the signal.
>>> cp.mean(Pxx_den[256:]) 0.0009924865443739191
Now compute and plot the power spectrum.
>>> f, Pxx_spec = cusignal.welch(x, fs, 'flattop', 1024, \ scaling='spectrum') >>> plt.figure() >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(cp.sqrt(Pxx_spec))) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('Linear spectrum [V RMS]') >>> plt.show()
The peak height in the power spectrum is an estimate of the RMS amplitude.
>>> cp.sqrt(Pxx_spec.max()) 2.0077340678640727
If we now introduce a discontinuity in the signal, by increasing the amplitude of a small portion of the signal by 50, we can see the corruption of the mean average power spectral density, but using a median average better estimates the normal behaviour.
>>> x[int(N//2):int(N//2)+10] *= 50. >>> f, Pxx_den = cusignal.welch(x, fs, nperseg=1024) >>> f_med, Pxx_den_med = cusignal.welch(x, fs, nperseg=1024, average='median') >>> plt.semilogy(cp.asnumpy(f), cp.asnumpy(Pxx_den), label='mean') >>> plt.semilogy(cp.asnumpy(f_med), cp.asnumpy(Pxx_den_med), \ label='median') >>> plt.ylim([0.5e-3, 1]) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('PSD [V**2/Hz]') >>> plt.legend() >>> plt.show()
Acoustics#
- cusignal.acoustics.cepstrum.complex_cepstrum(x, n=None, axis=-1)#
Calculates the complex cepstrum of a real valued input sequence x where the cepstrum is defined as the inverse Fourier transform of the log magnitude DFT (spectrum) of a signal. It’s primarily used for source/speaker separation in speech signal processing.
The input is altered to have zero-phase at pi radians (180 degrees) Parameters ———- x : ndarray
Input sequence, if x is a matrix, return cepstrum in direction of axis
- nint
Size of Fourier Transform; If none, will use length of input array
- axis: int
Direction for cepstrum calculation
Returns#
- cepsndarray
Complex cepstrum result
- cusignal.acoustics.cepstrum.inverse_complex_cepstrum(ceps, ndelay)#
Compute the inverse complex cepstrum of a real sequence. ceps : ndarray
Real sequence to compute inverse complex cepstrum of.
- ndelay: int
The amount of samples of circular delay added to x.
Returns#
- xndarray
The inverse complex cepstrum of the real sequence ceps.
The inverse complex cepstrum is given by .. math:: x[n] = F^{-1}left{exp(F(c[n]))right} where \(c_[n]\) is the input signal and \(F\) and :math:`F_{-1} are respectively the forward and backward Fourier transform.
- cusignal.acoustics.cepstrum.minimum_phase(x, n=None)#
Compute the minimum phase reconstruction of a real sequence. x : ndarray
Real sequence to compute the minimum phase reconstruction of.
- n{None, int}, optional
Length of the Fourier transform.
Compute the minimum phase reconstruction of a real sequence using the real cepstrum. Returns ——- m : ndarray
The minimum phase reconstruction of the real sequence x.
- cusignal.acoustics.cepstrum.real_cepstrum(x, n=None, axis=-1)#
Calculates the real cepstrum of an input sequence x where the cepstrum is defined as the inverse Fourier transform of the log magnitude DFT (spectrum) of a signal. It’s primarily used for source/speaker separation in speech signal processing
- Parameters:
- xndarray
Input sequence, if x is a matrix, return cepstrum in direction of axis
- nint
Size of Fourier Transform; If none, will use length of input array
- axis: int
Direction for cepstrum calculation
- Returns:
- cepsndarray
Complex cepstrum result
Wavelets#
Wavelets#
- cusignal.wavelets.wavelets.cwt(data, wavelet, widths)#
Continuous wavelet transform.
Performs a continuous wavelet transform on data, using the wavelet function. A CWT performs a convolution with data using the wavelet function, which is characterized by a width parameter and length parameter.
- Parameters:
- data(N,) ndarray
data on which to perform the transform.
- waveletfunction
Wavelet function, which should take 2 arguments. The first argument is the number of points that the returned vector will have (len(wavelet(length,width)) == length). The second is a width parameter, defining the size of the wavelet (e.g. standard deviation of a gaussian). See ricker, which satisfies these requirements.
- widths(M,) sequence
Widths to use for transform.
- Returns:
- cwt: (M, N) ndarray
Will have shape of (len(widths), len(data)).
Notes
length = min(10 * width[ii], len(data)) cwt[ii,:] = cusignal.convolve(data, wavelet(length, width[ii]), mode='same')
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt >>> t = cp.linspace(-1, 1, 200, endpoint=False) >>> sig = cp.cos(2 * cp.pi * 7 * t) + cusignal.gausspulse(t - 0.4, fc=2) >>> widths = cp.arange(1, 31) >>> cwtmatr = cusignal.cwt(sig, cusignal.ricker, widths) >>> plt.imshow(abs(cp.asnumpy(cwtmatr)), extent=[-1, 1, 31, 1], cmap='PRGn', aspect='auto', vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max()) >>> plt.show()
- cusignal.wavelets.wavelets.morlet(M, w=5.0, s=1.0, complete=True)#
Complex Morlet wavelet.
- Parameters:
- Mint
Length of the wavelet.
- wfloat, optional
Omega0. Default is 5
- sfloat, optional
Scaling factor, windowed from
-s*2*pi
to+s*2*pi
. Default is 1.- completebool, optional
Whether to use the complete or the standard version.
- Returns:
- morlet(M,) ndarray
See also
cusignal.gausspulse
Notes
The standard version:
pi**-0.25 * exp(1j*w*x) * exp(-0.5*(x**2))
This commonly used wavelet is often referred to simply as the Morlet wavelet. Note that this simplified version can cause admissibility problems at low values of w.
The complete version:
pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2))
This version has a correction term to improve admissibility. For w greater than 5, the correction term is negligible.
Note that the energy of the return wavelet is not normalised according to s.
The fundamental frequency of this wavelet in Hz is given by
f = 2*s*w*r / M
where r is the sampling rate.Note: This function was created before cwt and is not compatible with it.
- cusignal.wavelets.wavelets.morlet2(M, s, w=5)#
Complex Morlet wavelet, designed to work with cwt. Returns the complete version of morlet wavelet, normalised according to s:
exp(1j*w*x/s) * exp(-0.5*(x/s)**2) * pi**(-0.25) * sqrt(1/s)
Parameters#
- Mint
Length of the wavelet.
- sfloat
Width parameter of the wavelet.
- wfloat, optional
Omega0. Default is 5
Returns#
morlet : (M,) ndarray See Also ——– morlet : Implementation of Morlet wavelet, incompatible with cwt Notes —– .. versionadded:: 1.4.0 This function was designed to work with cwt. Because morlet2 returns an array of complex numbers, the dtype argument of cwt should be set to complex128 for best results. Note the difference in implementation with morlet. The fundamental frequency of this wavelet in Hz is given by:
f = w*fs / (2*s*np.pi)
where
fs
is the sampling rate and s is the wavelet width parameter. Similarly we can get the wavelet width parameter atf
:s = w*fs / (2*f*np.pi)
Examples#
>>> from scipy import signal >>> import matplotlib.pyplot as plt >>> M = 100 >>> s = 4.0 >>> w = 2.0 >>> wavelet = signal.morlet2(M, s, w) >>> plt.plot(abs(wavelet)) >>> plt.show() This example shows basic use of `morlet2` with `cwt` in time-frequency analysis: >>> from scipy import signal >>> import matplotlib.pyplot as plt >>> t, dt = np.linspace(0, 1, 200, retstep=True) >>> fs = 1/dt >>> w = 6. >>> sig = np.cos(2*np.pi*(50 + 10*t)*t) + np.sin(40*np.pi*t) >>> freq = np.linspace(1, fs/2, 100) >>> widths = w*fs / (2*freq*np.pi) >>> cwtm = signal.cwt(sig, signal.morlet2, widths, w=w) >>> plt.pcolormesh(t, freq, np.abs(cwtm), cmap='viridis', shading='gouraud') >>> plt.show()
- cusignal.wavelets.wavelets.qmf(hk)#
Return high-pass qmf filter from low-pass
- Parameters:
- hkarray_like
Coefficients of high-pass filter.
- cusignal.wavelets.wavelets.ricker(points, a)#
Return a Ricker wavelet, also known as the “Mexican hat wavelet”.
It models the function:
A (1 - x^2/a^2) exp(-x^2/2 a^2)
,where
A = 2/sqrt(3a)pi^1/4
.- Parameters:
- pointsint
Number of points in vector. Will be centered around 0.
- ascalar
Width parameter of the wavelet.
- Returns:
- vector(N,) ndarray
Array of length points in shape of ricker curve.
Examples
>>> import cusignal >>> import cupy as cp >>> import matplotlib.pyplot as plt
>>> points = 100 >>> a = 4.0 >>> vec2 = cusignal.ricker(points, a) >>> print(len(vec2)) 100 >>> plt.plot(cp.asnumpy(vec2)) >>> plt.show()
B-splines#
B-splines#
- cusignal.bsplines.bsplines.cubic(x)#
A cubic B-spline.
This is a special case of bspline, and equivalent to
bspline(x, 3)
.
- cusignal.bsplines.bsplines.gauss_spline(x, n)#
Gaussian approximation to B-spline basis function of order n.
- Parameters:
- nint
The order of the spline. Must be nonnegative, i.e. n >= 0
References
[1]Bouma H., Vilanova A., Bescos J.O., ter Haar Romeny B.M., Gerritsen F.A. (2007) Fast and Accurate Gaussian Derivatives Based on B-Splines. In: Sgallari F., Murli A., Paragios N. (eds) Scale Space and Variational Methods in Computer Vision. SSVM 2007. Lecture Notes in Computer Science, vol 4485. Springer, Berlin, Heidelberg
- cusignal.bsplines.bsplines.quadratic(x)#
A quadratic B-spline.
This is a special case of bspline, and equivalent to
bspline(x, 2)
.
Utilities#
Array Tools#
- cusignal.utils.arraytools.from_pycuda(pycuda_arr, device=0)#
Read in gpuarray from PyCUDA and output CuPy array
- Parameters:
- pycuda_arrPyCUDA gpuarray
- deviceint
GPU Device ID
- Returns:
- cupy_arrCuPy ndarray
- cusignal.utils.arraytools.get_pinned_array(data)#
Return populated pinned memory.
- Parameters:
- datacupy.ndarray or numpy.ndarray
The array to be copied to shared buffer
- strides: int or None
- order: char
- cusignal.utils.arraytools.get_pinned_mem(shape, dtype)#
Create a pinned memory allocation.
- Parameters:
- shapeint or tuple of ints
Output shape.
- dtypedata-type
Output data type.
- Returns:
- retndarray
Pinned memory numpy array.
Return populated shared memory between GPU and CPU.
- Parameters:
- datacupy.ndarray or numpy.ndarray
The array to be copied to shared buffer
- strides: int or None
- order: char
- streamint
Stream number (0 for default)
- portablebool
- wcbool
Return shared memory between GPU and CPU. Similar to numpy.zeros
- Parameters:
- shapendarray.shape
Size of shared memory allocation
- dtypecupy.dtype or numpy.dtype
Data type of allocation
- strides: int or None
- order: char
- streamint
Stream number (0 for default)
- portablebool
- wcbool
FFTPack Helper#
- cusignal.utils.fftpack_helper.next_fast_len(target)#
Find the next fast size of input data to fft, for zero-padding, etc.
SciPy’s FFTPACK has efficient functions for radix {2, 3, 4, 5}, so this returns the next composite of the prime factors 2, 3, and 5 which is greater than or equal to target. (These are also known as 5-smooth numbers, regular numbers, or Hamming numbers.)
- Parameters:
- targetint
Length to start searching from. Must be a positive integer.
- Returns:
- outint
The first 5-smooth number greater than or equal to target.
Notes
New in version 0.18.0.
Examples
On a particular machine, an FFT of prime length takes 133 ms:
>>> from scipy import fftpack >>> min_len = 10007 # prime length is worst case for speed >>> a = np.random.randn(min_len) >>> b = fftpack.fft(a)
Zero-padding to the next 5-smooth length reduces computation time to 211 us, a speedup of 630 times:
>>> fftpack.helper.next_fast_len(min_len) 10125 >>> b = fftpack.fft(a, 10125)
Rounding up to the next power of 2 is not optimal, taking 367 us to compute, 1.7 times as long as the 5-smooth size:
>>> b = fftpack.fft(a, 16384)
IO#
Reader#
- cusignal.io.reader.read_bin(file, buffer=None, dtype=<class 'numpy.uint8'>, num_samples=None, offset=0)#
Reads binary file into GPU memory. Can be used as a building blocks for custom unpack/pack data readers/writers.
- Parameters:
- filestr
A string of filename to be read to GPU.
- bufferndarray, optional
Pinned memory buffer to use when copying data from GPU.
- dtypedata-type, optional
Any object that can be interpreted as a numpy data type.
- num_samplesint, optional
Number of samples to be loaded to GPU. If set to 0, read in all samples.
- offsetint, optional
In the file, array data starts at this offset. Since offset is measured in bytes, it should normally be a multiple of the byte-size of dtype.
- Returns
- ——-
- outndarray
An 1-dimensional array containing binary data.
- cusignal.io.reader.read_sigmf(data_file, meta_file=None, buffer=None, num_samples=None, offset=0)#
Read and unpack binary file, with SigMF spec, to GPU memory.
- Parameters:
- data_filestr
File contain sigmf data.
- meta_filestr, optional
File contain sigmf meta.
- bufferndarray, optional
Pinned memory buffer to use when copying data from GPU.
- num_samplesint, optional
Number of samples to be loaded to GPU. If set to 0, read in all samples.
- offsetint, optional
May be specified as a non-negative integer offset. It is the number of samples before loading ‘num_samples’. ‘offset’ must be a multiple of ALLOCATIONGRANULARITY which is equal to PAGESIZE on Unix systems.
- Returns:
- outndarray
An 1-dimensional array containing unpacked binary data.
- cusignal.io.reader.unpack_bin(binary, dtype, endianness='L')#
Unpack binary file. If endianness is big-endian, it my be converted to little endian for NVIDIA GPU compatibility.
- Parameters:
- binaryndarray
The binary array to be unpack.
- dtypedata-type, optional
Any object that can be interpreted as a numpy data type.
- endianness{‘L’, ‘B’}, optional
Data set byte order
- Returns:
- outndarray
An 1-dimensional array containing unpacked binary data.
Writer#
- cusignal.io.writer.pack_bin(in1)#
Pack binary arrary. Data will be packed with little endian for NVIDIA GPU compatibility.
- Parameters:
- in1ndarray
The ndarray to be pack at binary.
- Returns:
- outndarray
An 1-dimensional array containing packed binary data.
- cusignal.io.writer.write_bin(file, binary, buffer=None, append=True)#
Writes binary array to file.
- Parameters:
- filestr
A string of filename to store output.
- binaryndarray
Binary array to be written to file.
- bufferndarray, optional
Pinned memory buffer to use when copying data from GPU.
- appendbool, optional
Append to file if created.
- Returns:
- outndarray
An 1-dimensional array containing binary data.
- cusignal.io.writer.write_sigmf(data_file, data, buffer=None, append=True)#
Pack and write binary array to file, with SigMF spec.
- Parameters:
- filestr
A string of filename to be read/unpacked to GPU.
- binaryndarray
Binary array to be written to file.
- bufferndarray, optional
Pinned memory buffer to use when copying data from GPU.
- appendbool, optional
Append to file if created.
- Returns:
Radar/Phased Array#
Radar Tools#
- cusignal.radartools.radartools.ambgfun(x, fs, prf, y=None, cut='2d', cutValue=0)#
Calculates the normalized ambiguity function for the vector x
- Parameters:
- xndarray
Input pulse waveform
- fs: int, float
Sampling rate in Hz
- prf: int, float
Pulse repetition frequency in Hz
- yndarray
Second input pulse waveform. If not given, y = x
- cutstring
Direction of one-dimensional cut through ambiguity function
- cutValueint, float
Time delay or doppler shift at which one-dimensional cut through ambiguity function is taken
- Returns:
- amfunndarray
Normalized magnitude of the ambiguity function
- cusignal.radartools.radartools.ca_cfar(array, guard_cells, reference_cells, pfa=0.001)#
Computes the cell-averaged constant false alarm rate (CA CFAR) detector threshold and returns for a given array.
- Parameters:
- arrayndarray
Array containing data to be processed.
- guard_cells_xint
One-sided guard cell count in the first dimension.
- guard_cells_yint
One-sided guard cell count in the second dimension.
- reference_cells_xint
one-sided reference cell count in the first dimension.
- reference_cells_yint
one-sided referenc cell count in the second dimension.
- pfafloat
Probability of false alarm.
- Returns:
- thresholdndarray
CFAR threshold
- returnndarray
CFAR detections
- cusignal.radartools.radartools.cfar_alpha(pfa, N)#
Computes the value of alpha corresponding to a given probability of false alarm and number of reference cells N.
- Parameters:
- pfafloat
Probability of false alarm.
- Nint
Number of reference cells.
- Returns:
- alphafloat
Alpha value.
- cusignal.radartools.radartools.pulse_compression(x, template, normalize=False, window=None, nfft=None)#
Pulse Compression is used to increase the range resolution and SNR by performing matched filtering of the transmitted pulse (template) with the received signal (x)
- Parameters:
- xndarray
Received signal, assume 2D array with [num_pulses, sample_per_pulse]
- templatendarray
Transmitted signal, assume 1D array
- normalizebool
Normalize transmitted signal
- windowarray_like, callable, string, float, or tuple, optional
Specifies the window applied to the signal in the Fourier domain.
- nfftint, size of FFT for pulse compression. Default is number of
samples per pulse
- Returns:
- compressedIQndarray
Pulse compressed output
- cusignal.radartools.radartools.pulse_doppler(x, window=None, nfft=None)#
Pulse doppler processing yields a range/doppler data matrix that represents moving target data that’s separated from clutter. An estimation of the doppler shift can also be obtained from pulse doppler processing. FFT taken across slow-time (pulse) dimension.
- Parameters:
- xndarray
Received signal, assume 2D array with [num_pulses, sample_per_pulse]
- windowarray_like, callable, string, float, or tuple, optional
Specifies the window applied to the signal in the Fourier domain.
- nfftint, size of FFT for pulse compression. Default is number of
samples per pulse
- Returns:
- pd_dataMatrixndarray
Pulse-doppler output (range/doppler matrix)
- cusignal.radartools.beamformers.mvdr(x, sv, calc_cov=True)#
Minimum variance distortionless response (MVDR) beamformer weights
- Parameters:
- xndarray
Received signal or input covariance matrix, assume 2D array with size [num_sensors, num_samples]
- sv: ndarray
Steering vector, assume 1D array with size [num_sensors, 1]
- calc_covbool
Determine whether to calculate covariance matrix. Simply put, calc_cov defines whether x input is made of sensor/observation data or is a precalculated covariance matrix
- Note: Unlike MATLAB where input matrix x is of size MxN where N represents
- the number of array elements, we assume row-major formatted data where each
- row is assumed to be complex-valued data from a given sensor (i.e. NxM)