Signal module

Signal module#

This module provides (mostly) orthonormal lazy linear operators associated to fast linear transforms commonly used in signal processing.

List of transforms#

lazylinop.signal.convolve() The lazy convolution operator is obtained as L = convolve(N, filter, ...) with N the dimension of the input signal and filter the filter (given as a vector), so that y = L @ x yields the prescribed convolution. Whenever needed z = L.H @ y computes the adjoint of the convolution operator, appropriately taking into account all boundary conditions that you may have specified using the “mode” option of convolve().
lazylinop.signal.dct() Discrete Cosine Transform of types I to IV
lazylinop.signal.dst() Discrete Sine Transform of types I to IV
lazylinop.signal.mdct() Modified DCT
lazylinop.signal.imdct() Inverse of the Modified DCT
lazylinop.signal.dft() Discrete Fourier Transform
lazylinop.signal.wht() Walsh-Hadamard Transform
lazylinop.signal.dwt() Discrete Wavelet Transform
lazylinop.signal.idwt() Inverse of the Discrete Wavelet Transform
lazylinop.signal.dzt() Zak Transform
lazylinop.signal.stft() Short-Time-Fourier-Transform with a given window
lazylinop.signal.istft() Inverse of Short-Time-Fourier-Transform
lazylinop.signal.nufft() Nonuniform Fast-Fourier-Transform

Padding or cropping#

Traditional implementations of signal processing transforms offer the possibility to either crop or zero-pad the analyzed signal. We chose not to include this in our implementation, as mimicing this feature is simple with the generic lazylinop interface. Consider for the example F = dft(N), the operator associated to the NxN DFT matrix. To apply it to a signal x of length n, we can define G = eye(N, n), and observe that F @ G does exactly what we need:

padding: if \(n<N\), G @ x is a zero padded version of x, of size \(N\), so that H = F @ G is the lazy linear operator that computes the DFT after zero padding. An even simpler implementation is H = F[:,:n].
cropping: if \(n>N\), G @ x is a cropped version of x of size \(N\), and H = F @ G is again exactly what you need.

Inverse transforms#

The lazy linear operator associated to a given transform (e.g., the DFT) is obtained as F = dft(N) with \(N\) the dimension of the input signal, and y = F @ x yields the DFT of x, a vector of size \(N\).

You may be surprised that we generally do not provide any implementation for inverse transforms (i.e., matching pairs of functions fft/ifft, dct/idct etc.). This is simply due to the fact that since most of the transforms we provide (to the exception of the STFT, MDCT, NUFFT and DWT when non-orthogonal wavelet is used) are orthonormal, their inverse is simply their adjoint, e.g. F.H is the lazy linear operator associated to the inverse DFT: z = F.H @ y recovers x.

Orthonormality corresponds to the fact that both F.H @ F and F @ F.H are lazy linear operators that act as the identity matrix (up to numerical precision). When F is a DCT/DST/WHT operator, its array version F.toarray() is real-valued, so you can also use F.T instead of F.H, and both F @ F.T and F.T @ F act as the identity.

The main exception to orthonormality are the MDCT, the NUFFT and the STFT: for example, the operator M = mdct(N) acts on signals of size N, and is only defined when N is even, with y = M @ x yielding a signal of size N/2. The matrix version of M is (real-valued and) rectangular of size \(\left[N/2,N\right]\) so that M.T @ M cannot correspond to an identity, however we chose to normalize M in such a way that both M @ M.T and M @ M.H correspond to the identity.

The other exception is the DWT: see its documentation (dwt()) for details on parameters ensuring that L = dwt(...) satisfies L @ L.T = L.T @ L = Id (or only L.T @ L = Id when the shape of L is rectangular).

Other normalizations#

Traditional implemementations of signal processing transforms offer various normalizations (e.g., with division by N, sqrt(N), or no division, or slightly more subtle operations for the DCT/DST). They can all be mimicked (if really needed) by pre- or post-composing the transform F. As an illustration, SciPy’s fft and ifft with the default normalization is mimicked as follows.

>>> import numpy as np
>>> from scipy.fft import fft as sp_fft
>>> from scipy.fft import ifft as sp_ifft
>>> from lazylinop.signal import dft as lz_dft
>>> N = 32
>>> x = np.random.randn(N)
>>> F = lz_dft(N)
>>> scale = sqrt(N)
>>> y = scale * F @ x
>>> np.allclose(y, sp_fft(x))
True
>>> x_ = F.H @ y / scale
>>> np.allclose(x_, sp_ifft(y))
True

To mimick SciPy’s DCT/DST called with orthogonalize=True, the same trick holds where scale depends on the transform’s type (I,II,III,IV) and the choice of norm ('backward' or 'forward'; NB: scale = 1 if norm = 'ortho').

Mimicking the DCT/DST with orthogonalize=False requires pre- and/or post-composing by diagonal operators that depend on the type. For example, the default DCT-II behavior (norm = 'ortho'):

>>> from lazylinop.signal import dct as lz_dct
>>> from scipy.fft import dct as sp_dct
>>> from lazylinop.basicops import diag
>>> import numpy as np
>>> N = 32
>>> x = np.random.randn(N)
>>> F = lz_dct(N)
>>> v = np.full(N, 1.0)
>>> v[0] = np.sqrt(2.0)
>>> y = diag(v) @ F @ x
>>> z = sp_dct(x, 2, N, 0, 'ortho', False, 1, orthogonalize=False)
>>> np.allclose(y, z)

Note: to be coherent with other transforms, our implementation of the Walsh-Hadamard Transform is orthonormal, unlike the Hadamard function of Scipy.

Transforms#

lazylinop.signal.convolve(N, filter, mode='full', backend='scipy_convolve')#

Returns a LazyLinOp L for the 1d convolution of signal(s) of size N with a kernel filter.

Shape of L is \((M,~N)\). See mode for output size \(M\).

Args:

N: int

Size of the signal to be convolved.

filter: np.ndarray, cp.ndarray or torch.Tensor

Filter to be applied to the signal. filter must be a 1d-array.

mode: str, optional

'full': compute full convolution including at points of non-complete overlapping of inputs. Output size M is N + filter.size - 1.
'valid': compute only fully overlapping part of the convolution. This is the ‘full’ center part of (output) size M = max(N, K) - min(N, K) + 1 where K = filter.size.
'same': compute only the center part of 'full' to obtain an output size M equal to the input size N.
'circ': compute circular convolution (filter.size <= N must be satisfied). Output size M is N.

backend: str, optional

'auto': use best backend available according to the namespace of filter.
- for mode='circ' use backend='fft'.
- for any other mode use backend='direct' if filter is a NumPy array and filter.size < log2(M), backend='scipy_convolve' if filter is a NumPy array and filter.size >= log2(M), backend='scipy_convolve' if filter is a CuPy or NumPy array, backend='torch' if filter is a torch tensor.
'scipy_convolve': (default) encapsulates (cupyx).scipy.signal.convolve.

It uses internally the best SciPy backend between 'fft' and 'direct' (see scipy.signal.choose_conv_method). If the filter and the input are CuPy arrays the function uses cupyx.scipy.signal.convolve.
'direct': direct computation using nested for-loops with Numba and parallelization. It does not work if the input is CuPy array or torch tensor.
'fft': (only for mode='circ') compute circular convolution using (cupyx).scipy.fft.fft.
'toeplitz': encapsulate (cupyx).scipy.linalg.toeplitz if N < 2048, (cupyx).scipy.linalg.matmul_toeplitz otherwise. Of note, there is no torch implementation for Toeplitz matrix.
'oa': use Lazylinop implementation of overlap-add backend.
'torch': use torchaudio.functional.convolve.

Returns:

LazyLinOp

Examples:

>>> import numpy as np
>>> import scipy.signal as sps
>>> import lazylinop.signal as lzs
>>> N = 1024
>>> x = np.random.randn(N)
>>> kernel = np.random.randn(32)
>>> L = lzs.convolve(N, kernel)
>>> c1 = L @ x
>>> c2 = sps.convolve(x, kernel)
>>> np.allclose(c1, c2)
True
>>> N = 32768
>>> x = np.random.randn(N)
>>> kernel = np.random.randn(48)
>>> L = lzs.convolve(N, kernel, mode='circ', backend='fft')
>>> c1 = L @ x
>>> L = lzs.convolve(N, kernel, mode='circ', backend='direct')
>>> c2 = L @ x
>>> np.allclose(c1, c2)
True

See also

lazylinop.signal.dct(N, type=2, backend='scipy')#

Returns a LazyLinOp `L for the Direct Cosine Transform (DCT).

Shape of L is \((N,~N)\).

L is orthonormal, and the LazyLinOp for the inverse DCT is L.T.

The function provides two backends: SciPy and Lazylinop.

Args:

N: int

Size of the input (N > 0).

type: int, optional

1, 2, 3, 4 (I, II, III, IV). Defaut is 2. See SciPy DCT and CuPy DCT for more details.

backend: str, optional

'scipy' (default) uses (cupyx).scipy.fft.dct: encapsulation for the underlying computation of the DCT.

'lazylinop' uses a composition of basic Lazylinop operators to compute the DCT (fft(), vstack() etc.).

Returns:

LazyLinOp

Example:

>>> from lazylinop.signal import dct as lz_dct
>>> from scipy.fft import dct as sp_dct
>>> import numpy as np
>>> N = 32
>>> x = np.random.randn(N)
>>> F = lz_dct(N)
>>> y = F @ x
>>> np.allclose(y, sp_dct(x, norm='ortho'))
True
>>> # compute the inverse DCT
>>> x_ = F.T @ y
>>> np.allclose(x_, x)
True
>>> # To mimick SciPy DCT II norm='ortho' and orthogonalize=False
>>> from lazylinop.basicops import diag
>>> v = np.full(N, 1.0)
>>> v[0] = np.sqrt(2.0)
>>> y = diag(v) @ F @ x
>>> z = sp_dct(x, 2, N, 0, 'ortho', False, 1, orthogonalize=False)
>>> np.allclose(y, z)
True

References:

[1] A Fast Cosine Transform in One and Two Dimensions,: by J. Makhoul, IEEE Transactions on acoustics, speech and signal processing vol. 28(1), pp. 27-34, :doi:`10.1109/TASSP.1980.1163351` (1980).

See also

lazylinop.signal.dst(N, type=2, backend='scipy')#

Returns a LazyLinOp L for the Direct Sine Transform (DST).

Shape of L is \((N,~N)\).

L is orthonormal, and the LazyLinOp for the inverse DST is L.T.

The function provides two backends: SciPy and Lazylinop.

Args:

N: int

Size of the input (N > 0).

type: int, optional

1, 2, 3, 4 (I, II, III, IV). Defaut is 2. See SciPy DST and CuPy DST for more details.

backend: str, optional

'scipy' (default) uses (cupyx).scipy.fft.dst to compute the DST depending on the input array.
'lazylinop' Uses a composition of basic Lazylinop operators to compute the DST (fft(), vstack() etc.).

Returns:

LazyLinOp

Example:

>>> from lazylinop.signal import dst as lz_dst
>>> from scipy.fft import dst as sp_dst
>>> import numpy as np
>>> x = np.random.randn(32)
>>> F = lz_dst(32)
>>> y = F @ x
>>> np.allclose(y, sp_dst(x, norm='ortho'))
True
>>> # compute the inverse DST
>>> x_ = F.T @ y
>>> np.allclose(x_, x)
True

See also

lazylinop.signal.mdct(N, window=('vorbis', 128), backend='scipy')#

Returns a LazyLinOp L for the Modified Direct Cosine Transform (MDCT).

Shape of L is \((M,~N)\) with \(M=n\frac{W}{2}\).

\(n\) is the number of chunks and \(W\) is the window size.

Our implementation matches the one from TensorFlow, L @ x where L = mdct(N, window=('vorbis', 128)) is equivalent to:

import tensorflow as tf
tf.signal.mdct(x,
               frame_length=128,
               window_fn=tf.signal.vorbis_window,
               norm='ortho')

If we consider MDCT using a rectangular window with a scale \(\frac{1}{\sqrt{2}}\), L = mdct(N, window=('None', N)), X = L @ x yiels a vector X of size \(H=\frac{N}{2}\) such that:

\[\begin{equation} X_k=\frac{1}{\sqrt{H}}\sum_{n=0}^{N-1}x_n\cos\left(\frac{\pi}{H}\left(n+\frac{1}{2}+\frac{H}{2}\right)\left(k+\frac{1}{2}\right)\right) \end{equation}\]

y = L.T @ X yields a vector y of size \(N\) such that:

\[\begin{equation} y_n=\frac{1}{\sqrt{H}}\sum_{k=0}^{H-1}X_k\cos\left(\frac{\pi}{H}\left(n+\frac{1}{2}+\frac{H}{2}\right)\left(k+\frac{1}{2}\right)\right) \end{equation}\]

The operator L is rectangular and is not left invertible. It is however right-invertible as L @ L.T = Id. Thus, L.T can be used as a right-inverse.

After removing some details the code looks like:

# Consecutive windows are slided by hop samples.
hop = win_size // 2
# Number of chunks.
n = 1 + (N - win_size) // hop
# Lazy linear operator that extracts and
# concatenates chunks from a signal of length N.
G = chunk(N, win_size, hop)

and using the mixed Kronecker product property \((I^T\otimes A)\mathtt{vec}(X)=\mathtt{vec}(AXI)=\mathtt{vec}(AX)\):

K = kron(eye(n), _helper(win_size, backend) @ diag(win))
L = K @ G

where _helper(...) encapsulates underlying implementation using DCT of type IV (see Ref [1] for more details).

The function provides two backends: SciPy and Lazylinop for the underlying computation of the DCT of type IV.

Args:

N: int

Length of the input array.

window: (str, int) or (str, int, float), optional

Window, a tuple (name: str, win_size: int) or (name: str, win_size: int, beta: float). Window size must be a mutliple of 4. Default is ('vorbis', 128). beta has no effect excepts for 'kaiser_bessel_derived' window. Possible windows are:

'None' corresponds to a rectangular window with a scale \(\frac{1}{\sqrt{2}}\).
'kaiser_bessel_derived' see scipy.signal.window.kaiser_bessel_derived for more details.
'vorbis' (default) or 'sin' see https://en.wikipedia.org/wiki/Modified_discrete_cosine_transform for more details.

backend: str, optional

'scipy' (default) uses scipy.fft.dct encapsulation for the underlying computation of the DCT of type IV. backend='scipy' does not work with CuPy input array because cupyx.signal.fft.dct only implements type II and III.
'lazylinop' uses pre-built Lazylinop operators (Lazylinop dct(), eye(), kron(), vstack() etc.) to build the pipeline that will compute the MDCT and the underlying DCT of type IV.

Returns:

LazyLinOp

Example:

>>> from lazylinop.signal import mdct
>>> import numpy as np
>>> x = np.random.randn(64)
>>> # MDCT with a rectangular window
>>> # of size equal to the size of the input.
>>> L = mdct(64, ('None', 64))
>>> y = L @ x
>>> y.shape[0] == 32
True
>>> x_ = L.T @ y
>>> x_.shape[0] == 64
True

References:

[1] Xuancheng Shao, Steven G. Johnson, Type-IV DCT, DST, and MDCT algorithms with reduced numbers of arithmetic operations, Signal Processing, Volume 88, Issue 6, 2008, Pages 1313-1326, ISSN 0165-1684, https://doi.org/10.1016/j.sigpro.2007.11.024.

See also

lazylinop.signal.imdct(N, window=('vorbis', 128), backend='scipy')#

Returns a LazyLinOp iL for the inverse Modified Direct Cosine Transform (iMDCT).

Shape of iL is \((P,~M)\) with \(M=n\frac{W}{2}\) and \(P=(n-1)\frac{W}{2}+W\).

\(n\) is the number of chunks and \(W\) is the window size.

N is the size of the input signal of the MDCT LazyLinop, not of the iMDCT.

Our implementation matches the one from TensorFlow, iL @ y where iL = imdct(N, window=('vorbis', N)) is equivalent to:

import tensorflow as tf
tf.signal.inverse_mdct(y,
                       window_fn=tf.signal.vorbis_window,
                       norm='ortho')

The operator iL is rectangular and is not right invertible. It is however left-invertible as iL.T @ iL. Thus, iL.T can be used as a left-inverse.

After removing some details the code looks like:

# Consecutive windows are slided by hop samples.
hop = win_size // 2
# Number of chunks.
n = 1 + (N - win_size) // hop

if N == win_size:
    return diag(win) @ _helper(N, backend)
else:
    # Apply inverse MDCT per chunk followed by windowing.
    K = kron(eye(n), diag(win) @ _helper(win_size, backend))
    # Overlap and add.
    A = overlap_add(K.shape[0], size=win_size, overlap=win_size - hop)
    # Restrict to get the output of length N.
    iL = eye(N, A.shape[0]) @ A @ K

where _helper(...) encapsulates underlying implementation using DCT of type IV (see Ref [1] for more details).

The function provides two backends: SciPy and Lazylinop for the underlying computation of the DCT of type IV.

Args:

N: int

Length of the output array (see above).

window: (str, int) or (str, int, float), optional

'None' corresponds to a rectangular window with a scale \(\frac{1}{\sqrt{2}}\).
'kaiser_bessel_derived' see scipy.signal.window.kaiser_bessel_derived for more details.
'vorbis' (default) or 'sin' see https://en.wikipedia.org/wiki/Modified_discrete_cosine_transform for more details.

backend: str, optional

'scipy' (default) uses scipy.fft.dct encapsulation for the underlying computation of the DCT of type IV. backend='scipy' does not work with CuPy input array because cupyx.signal.fft.dct only implements type II and III.
'lazylinop' uses pre-built Lazylinop operators (Lazylinop dct(), eye(), kron(), vstack() etc.) to build the pipeline that will compute the MDCT and the underlying DCT of type IV.

Returns:

LazyLinOp

Example:

>>> from lazylinop.signal import mdct, imdct
>>> import numpy as np
>>> N = 64
>>> x = np.random.randn(N)
>>> # MDCT with a rectangular window
>>> # of size equal to the size of the input.
>>> L = imdct(N, window=('None', N))
>>> L.shape[0] == 64
True
>>> L.shape[1] == 32
True

References:

[1] Xuancheng Shao, Steven G. Johnson, Type-IV DCT, DST, and MDCT algorithms with reduced numbers of arithmetic operations, Signal Processing, Volume 88, Issue 6, 2008, Pages 1313-1326, ISSN 0165-1684, https://doi.org/10.1016/j.sigpro.2007.11.024.

See also

lazylinop.signal.dft(N)#

Returns a LazyLinOp `L for the Discrete Fourier Transform (DFT).

Shape of L is \((N,~N)\).

L is orthonormal, and the LazyLinOp for the inverse DFT is L.H.

SciPy FFT is used as underlying implementation.

Of note, we provide an alias fft of the dft function.

Args:

N: int: Size of the input (\(N > 0\)).

Returns:

LazyLinOp DFT

Examples:

>>> from lazylinop.signal import dft as lz_dft
>>> from scipy.fft import fft as sp_fft
>>> from scipy.fft import ifft as sp_ifft
>>> import numpy as np
>>> N = 32
>>> x = np.random.randn(N)
>>> F = lz_dft(N)
>>> y = F @ x
>>> np.allclose(y, sp_fft(x, norm='ortho'))
True
>>> # easy inverse
>>> x_ = F.H @ y
>>> np.allclose(x_, x)
True
>>> # To mimick SciPy FFT norm='backward'
>>> scale = np.sqrt(N)
>>> y = scale * F @ x
>>> np.allclose(y, sp_fft(x))
True
>>> x_ = F.H @ y / scale
>>> np.allclose(x_, sp_ifft(y))
True

See also

scipy.fft.fft, scipy.fft.ifft, rfft()

lazylinop.signal.wht(N, backend='auto')#

Returns a LazyLinOp `L for the Fast-Walsh-Hadamard-Transform (WHT).

Shape of L is \((N,~N)\). N must be a power of 2.

L is orthonormal and symmetric and the inverse WHT operator is L.T = L.

y = L @ x returns a vector of size \(N\) that satisfies y = scipy.linalg.hadamard(N) @ x / np.sqrt(N). The implementation is however more efficient for large values of \(N\).

Of note, we provide an alias fwht of the wht function.

Args:

N: int

Size of the input (N > 0).

backend: str, optional

Either ‘auto’ (default), ‘direct’, ‘kronecker’, ‘pyfaust’ or ‘scipy’.

‘auto’ uses ‘pyfaust’ for N<=32768, ‘direct’ otherwise. Of note, this condition depends on your hardware. It may change in the near future.
‘direct’ uses N * log(N) algorithm and Numba jit to compute FWHT. Because of Numba jit, computation of the first F @ x will be longer than the next one.
- Larger the size N of the input is better the performances are.
- Larger the batch is better the performances are.
‘kronecker’ uses kron from lazylinop to compute FWHT.
‘pyfaust’ (default) uses wht from pyfaust.
‘scipy’ uses scipy.linalg.hadamard matrix. It could be memory consuming for large N.
‘cupy’ uses cupyx.scipy.linalg.hadamard matrix. It could be memory consuming for large N.

Returns:

LazyLinOp

Examples:

>>> import numpy as np
>>> import scipy.linalg as spl
>>> import lazylinop.signal as lzs
>>> N = 16
>>> x = np.random.randn(N)
>>> H = lzs.fwht(N)
>>> y = np.sqrt(N) * (H @ x)
>>> np.allclose(y, spl.hadamard(N) @ x)
True
>>> np.allclose(H.T @ (H @ x), x)
True
>>> X = np.random.randn(N, 3)
>>> Y = np.sqrt(N) * (H @ X)
>>> np.allclose(Y, spl.hadamard(N) @ X)
True

See also

lazylinop.signal.dwt(N, wavelet='haar', mode='zero', level=None, backend='pywavelets')#

Returns a LazyLinOp L for the Discrete Wavelet Transform (DWT). L @ x will return a 1D numpy array as the concatenation of the DWT coefficients in the form [cAₙ, cDₙ, cDₙ₋₁, …, cD₂, cD₁], where n is the decomposition level, cAₙ is the approximation coefficients for level n and cDᵢ is the detail coefficients at decomposition level i.

Shape of L is \((M,~N)\) where \(M\) depends on the wavelet, input size N and decomposition level.

In general, L is not orthogonal. However, when wavelet is either the name of an orthogonal wavelet (see Pywavelets documentation), or a tuple of Quadrature-Mirrors-Filters, the following properties holds: - if mode='zero' then L.T is a left-inverse to L,

i.e. L.T @ L = Id (but since L.shape = (M, N) with \(M>N\), L does not satisfy L @ L.T = Id.

if mode='periodization' and if the signal size N is a power of two, then \(M=N\) and L is orthogonal so that L @ L.T = L.T @ L = Id

See technical notes for more details.

After removing some details the code looks like (for the first level of decomposition and without extension mode):

# Convolution operators with low-pass and
# high-pass decomposition filters.
G = convolve(N, dec_lo, mode='same')
H = convolve(N, dec_hi, mode='same')
# Downsampling operator.
DG = downsample(G.shape[0], 2, 1)
DH = downsample(H.shape[0], 2, 1)
# Decomposition operator.
L = vstack((DG @ G, DH @ H))

Args:

N: int: Size of the input array.

wavelet: str or tuple of (dec_lo, dec_hi), optional

If a str is provided, the wavelet name from Pywavelets library. In that case, the namespace of the low and high-pass filters depends on the namespace of the input array.

If a tuple (dec_lo, dec_hi) of two NumPy, CuPy arrays or torch tensors. is provided, the low and high-pass filters (for decomposition) used to define the wavelet.

The dwt() function does not test whether these two filters are actually Quadrature-Mirror-Filters.

mode: str, optional

'zero', signal is padded with zeros (default).

'periodic', signal is treated as periodic signal.

'symmetric', use mirroring to pad the signal.

'antisymmetric', signal is extended by mirroring and multiplying elements by minus one.

'reflect', signal is extended by reflecting elements.

'periodization', signal is extended like 'periodic' extension mode. Only the smallest possible number of coefficients is returned. Odd-length signal is extended first by replicating the last value.

level: int, optional: Decomposition level, by default (None) the maximum level is used.
backend: str, optional: 'pywavelets' (default) or 'lazylinop' for the underlying computation of the DWT. backend='pywavelets' does not work for CuPy array and torch tensor inputs.

Returns:

LazyLinOp DWT.

Examples:

>>> from lazylinop.signal import dwt
>>> import numpy as np
>>> import pywt
>>> N = 8
>>> x = np.arange(1, N + 1, 1)
>>> L = dwt(N, wavelet='haar', mode='periodic', level=1)
>>> y = L @ x
>>> z = pywt.wavedec(x, wavelet='haar', mode='periodic', level=1)
>>> np.allclose(y, np.concatenate(z))
True

See also

lazylinop.signal.idwt(N, wavelet='haar', mode='zero', level=None, backend='pywavelets')#

Returns a LazyLinOp iL for the inverse Discrete Wavelet Transform (iDWT).

If L = dwt(N, wavelet, mode, level, backend) is the DWT operator of shape \((M,~N)\) with \(M\ge~N\), then iL = idwt(N, wavelet, mode, level, backend) is the iDWT operator such that iL @ L = Id. As a result, if y = L @ x is the coefficients at level decomposition level, then the N-dimensionnal signal x can be reconstructed from the M-dimensionnal vector y as iL @ y.

Shape of iL is \((N,~M)\).

N is the size of the input signal of the DWT LazyLinop, not of the iDWT.
The order is not a typo since N is the input length of the associated DWT operator L, which is of shape \((M,~N)\), and iL is of the same shape as L.H.

In general, iL is not orthogonal. See technical notes for more details.

The implementation uses pre-built Lazylinop operators to build iDWT operator.

After removing some details the code looks like (for the first level of decomposition):

# Wavelet length.
W = pywt.Wavelet(wavelet).dec_len
# Number of detail and approximation coefficients.
m = (N + W - 1) // 2
# Upsampling operator.
U = downsample(2 * m, step=2, start=0).H
# Convolution operators with low-pass and
# high-pass reconstruction filters.
Cl = convolve(U.shape[0], rec_lo, mode='same')
Ch = convolve(U.shape[0], rec_hi, mode='same')
C = hstack((Cl @ U, Ch @ U))
# Restrict.
R = eye(N, C.shape[0], k=(C.shape[0] - N) // 2)
# Reconstruction operator.
iL = R @ C

Reconstruction using 'dmey' wavelet does not ensure np.allclose(x, L @ y) == True.

import numpy as np
import pywt
N = 133
x = np.random.randn(N)
y = pywt.wavedec(x, wavelet='dmey', level=1, mode='periodic')
z = pywt.waverec(y, wavelet='dmey', mode='periodic')
np.allclose(x, z[:N])
False
y = pywt.wavedec(x, wavelet='coif11', level=1, mode='periodic')
z = pywt.waverec(y, wavelet='coif11', mode='periodic')
np.allclose(x, z[:N])
True

Args:

N: int: Length of the output array (i.e., length of the input array of the associated DWT LazyLinOp, see above).

wavelet: str or tuple of (rec_lo, rec_hi), optional

If a str is provided, the wavelet name from Pywavelets library

If a tuple (rec_lo, rec_hi) of two NumPy, CuPy arrays or torch tensors is provided, the low and high-pass filters (for reconstruction) used to define the wavelet.

The idwt() function does not test whether these two filters are actually Quadrature-Mirror-Filters.

mode: str, optional

'zero', signal is padded with zeros (default).

'periodic', signal is treated as periodic signal.

'periodization', signal is extended like 'periodic' extension mode. Only the smallest possible number of coefficients is returned. Odd-length signal is extended first by replicating the last value.

'symmetric', use mirroring to pad the signal.

'antisymmetric', signal is extended by mirroring and multiplying elements by minus one.

'reflect', signal is extended by reflecting elements.

level: int, optional: Decomposition level, by default (None) the maximum level is used.
backend: str, optional: ‘pywavelets’ (default) or ‘lazylinop’ for the underlying computation of the iDWT. backend='pywavelets' does not work for CuPy array and torch tensor inputs.

Returns:

LazyLinOp iDWT.

Examples:

>>> from lazylinop.signal import dwt, idwt
>>> import numpy as np
>>> N = 8
>>> x = np.arange(1, N + 1, 1)
>>> W = dwt(N, mode='periodic', level=2)
>>> y = W @ x
>>> L = idwt(N, mode='periodic', level=2)
>>> x_ = L @ y
>>> np.allclose(x, x_)
True

See also

lazylinop.signal.dzt(N, a=1)#

Returns a LazyLinOp L for the Discrete Zak Transform (DZT). The parameter a must divide input length N.

\[\begin{equation} Z_x^a\left(n, k\right)=\frac{1}{\sqrt{a}}\sum_{l=0}^ax\left(n+lL\right)e^{-2j\pi\frac{kl}{a}} \end{equation}\]

Shape of L is \((N,~N)\).

L is orthonormal, and the LazyLinOp for the inverse DZT is L.H.

Args:

N: int: Size of the input (\(N > 0\)).
a: int, optional: The parameter a must divide input length N. Default is 1.

Returns:

LazyLinOp DZT

Examples:

>>> from lazylinop.signal import dzt
>>> import numpy as np
>>> Z = dzt(32, 1)
>>> x = np.random.randn(32)
>>> y = Z @ x
>>> x_ = Z.H @ y
>>> np.allclose(x_, x)
True

References:

[1] H. Bölcskei and F. Hlawatsch. Discrete Zak transforms, polyphase: transforms, and applications. IEEE Trans. Signal Process., 45(4):851–866, april 1997.

See also

Zak transform (Wikipedia), LTFAT.

lazylinop.signal.stft(N, window=('hann', 256), hop=None, fft_size=None)#

Returns a LazyLinOp for the Short-Time-Fourier-Transform (STFT), so that L @ x is equivalent to:

import scipy as sp
win_size = window[1]
STFT = sp.signal.ShortTimeFFT(
    win=sp.signal.windows.get_window(window, win_size),
    hop=hop,
    fs=1.0,
    fft_mode='twosided',
    fft_size=None,
    phase_shift=None
)
STFT.stft(x)

The STFT computes and concatenates Fourier transforms of windowed signal slices where consecutive windows are slided by hop samples.

The implementation uses pre-built Lazylinop operators to build the STFT operator.

After removing some details the code looks like:

# For a 1d real valued array
win = sp.signal.get_window(win_name, win_size)
# Number of chunks.
n = 1 + (N - win_size) // hop
# Lazy linear operator that extracts and
# concatenates chunks from a signal of length N.
G = chunk(N, win_size, hop)
# Lazy operator that multiplies a (batch of)
# chunk(s) of length win_len by the window.
W = diag(win)
# Lazy operator computing the DFT of a (batch of)
# windowed chunk(s) of length win_len with appropriate padding.
FW = dft(fft_size)[:win_size, :] @ W

and using the mixed Kronecker product property \((I^T\otimes A)\mathtt{vec}(X)=vec(AXI)=vec(AX)\):

# Final lazy operator that extract chunks, applies F on each chunk.
L = kron(eye(n), FW) @ G

Shape of L is \((M,~N)\). In general \(M\ge~N\).

Args:

N: int

Length of the input array.

window: NumPy/CuPy array, torch tensor or (str, int), optional

Window, either directly provided as an 1d (real or complex valued) array, or as a pair (name: str, win_size: int) which is equivalent to use window=scipy.signal.get_window(name, length). Default is ('hann', 256). See scipy.signal.get_window for a list of available windows.

Be aware that if window does not have the same namespace xp = array_api_compat.array_namespace(x) than x, computation of L @ x uses an extra xp.asarray(window, device=array_api_compat.device(x), copy=True). This incur a loss of performance.

hop: int, optional

The increment in sample location between two consecutives slices. None (default) corresponds to len(window) // 2.

fft_size: int, optional: Size of the FFT, if zero-padding is required. Default value is None. It corresponds to fft_size=len(window).

Returns:

LazyLinOp L for the STFT

Examples:

>>> import lazylinop as lz
>>> import numpy as np
>>> import scipy as sp
>>> x = np.random.randn(1024)
>>> L = lz.signal.stft(x.shape[0])
>>> win = sp.signal.windows.get_window('hann', 256)
>>> STFT = sp.signal.ShortTimeFFT(win=win, hop=128, fs=1.0, fft_mode="twosided", phase_shift=None)
>>> np.allclose(STFT.stft(x).ravel(order='F'), L @ x)
True

See also

lazylinop.signal.istft(N, window=('hann', 256), hop=None, fft_size=None)#

Returns a LazyLinOp LazyLinOp iL for the inverse Short-Time-Fourier-Transform (iSTFT), so that if L is the LazyLinop of the STFT with the same arguments, iL @ L is the identity.

Shape of iL is \((N,~M)\) with \(M\ge~N\).

N is the size of the input signal of the STFT LazyLinop, not of the iSTFT.
The order is not a typo since N is the input length of the associated STFT operator L, which is of shape \((M,~N)\), and iL is of the same shape as L.H.

Args:

N: int

Length of the output array (i.e., length of the input array of the associated STFT LazyLinOp, see above).

window: NumPy array, CuPy array, torch tensor or (str, int), optional

Window, either directly provided as an 1d (real or complex valued) array , or as a pair (name: str, win_size: int) which is equivalent to use window=scipy.signal.get_window(name, length). Default is ('hann', 256). See scipy.signal.get_window for a list of available windows.

hop: int, optional

The increment in sample location between two consecutives slices. None (default) corresponds to len(window) // 2.

fft_size: int, optional: See stft() for more details.

Returns:

LazyLinOp iL for the inverse STFT

Examples:

>>> import lazylinop as lz
>>> import numpy as np
>>> N = 1024
>>> x = np.random.randn(N)
>>> L = lz.signal.stft(N)
>>> y = L @ x
>>> iL = lz.signal.istft(N)
>>> np.allclose(x, iL @ y)
True

See also

lazylinop.signal.nufft(N, samples=None, frequencies=None, backend='lazylinop', eps=1e-06)#

Returns a LazyLinOp L for the Non-Uniform Discrete Fourier Transform (NUFFT).

Shape of L is \((N,~N)\).

In general, L.H is not the inverse NUFFT (as the generic NUFFT is not orthonormal).

L = nufft(N, None, None) (default) is functionally equivalent to the Fast-Fourier-Transform operator and returns F = fft(N).

The nonuniform Discrete Fourier Transform is defined (for possibly nonuniform samples \(s_n\) and possibly nonuniform frequencies \(\omega_k\)) as:

\[\begin{equation} y_k\equiv\left(\omega_k\right)~\text{where}~y\left(\omega\right)=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x_ne^{-2\pi is_n\omega}\text{,}~\forall\omega\in\left[0,~N\right] \end{equation}\]

Args:

N: int

Size of the input (\(N > 0\)).

samples: 1d array

Samples \(s_n\in\left[0,~1\right[\) must be real. The size of samples must be \(N\). Default value is None (uniform samples \(s_n=\frac{n}{N}\) with \(0\le n\lt N\)).

frequencies: 1d array

Frequencies \(\omega_k\in\left[0,~N\right[\) must be real. The size of frequencies must be \(N\). Default value is None (uniform frequencies \(\omega_k=k\) with \(0\le k\lt N\)).

backend: str, optional

Three different implementations are available:

'lazylinop' is the default value. The underlying implementation is based on [1] and LazyLinOp’s pre-built operators.
'finuftt' uses encapsulation of fiNUFFT Python package.
'direct' uses direct computation of NUFFT. In the case of small \(N\) there is no need for a fast algorithm. eps argument has no effect for backend='direct'.

eps: float, optional

Precision to achieve (see [1] and fiNUFFT Python package for more details). Default value is 1e-6. It has not effect for backend='direct'.

Returns:

LazyLinOp NUFFT

Examples:

>>> from lazylinop.signal import nufft
>>> import numpy as np
>>> N = 32
>>> x = np.random.randn(N)
>>> samples = np.random.rand(N)
>>> k = 4
>>> F = nufft(N, samples, None)
>>> y = F @ x
>>> # Naive implementation of NUFFT-II.
>>> z = np.sum(x * np.exp(-2j * np.pi * k * samples)) / np.sqrt(N)
>>> np.allclose(y[k], z)
True

References:

[1] A Nonuniform Fast Fourier Transform Based on Low Rank Approximation. Diego Ruiz-Antolin and Alex Townsend SIAM Journal on Scientific Computing 2018 40:1, A529-A547

Utility functions#

lazylinop.signal.chunk(N, size, hop=1, start=0, stop=None)#

Returns a LazyLinOp that extracts chunks from a signal of size N.

This operator returns the concatenation of regularly spaced chunks of size elements from the provided signal. The spacing between consecutive chunks is of hop elements of the original signal. If start and/or stop is provided, the first chunk starts at the start element of the signal (and subsequents chunks are also shifted) and stops at the stop element.

Args:

size: int: Size of each chunk.
hop: int, optional: Number of elements in the original signal between two chunks. Default is 1.
start: int, optional: Index of the element in the original signal where chunking starts. Default is 0.
stop: int, optional: Index of the element (excluded) in the original signal where chunking stops. Default is None, meaning that chunkins stops at the end of the signal.

Returns:

The chunk LazyLinOp of shape \((M,N)\), where \(M=size imes Nchunks\) with \(Nchunks\) the number of chunks.

Chunk description#

Examples:

>>> import numpy as np
>>> import lazylinop as lz
>>> N = 10
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> # extract chunks of 5 elements every 2 element
>>> lz.signal.chunk(N, 5, 2) @ x
array([0, 1, 2, 3, 4, 2, 3, 4, 5, 6, 4, 5, 6, 7, 8])

lazylinop.signal.decimate(N, step, start=0, stop=None)#

Returns a LazyLinOp that decimates/subsamples a signal x of size N.

This operator performs decimation on provided signal, where the result contains elements selected every step elements in the original signal, starting at start up to stop, if provided.

When L = decimate(N, step, start, stop), computing y = L@x is equivalent to y = x[start:stop:step].

Args:

N: int: Length of the input vector.
step: int: The stride to use for selecting elements in the original signal.
start: int, optional: Index of the element in the original signal where decimation starts. Default is 0.
stop: int, optional: Index of the element (excluded) in the original signal where decimation stops. Default is None, meaning that chunkins stops at the end of the signal.

Returns:

The decimate LazyLinOp.

Examples:

>>> import numpy as np
>>> import lazylinop as lz
>>> N = 10
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> # decimate every 2 element, starting at the third
>>> lz.signal.decimate(N, 2, 3) @ x
array([3, 5, 7, 9])

lazylinop.signal.downsample(N, step, start=0, stop=None)#

Returns a LazyLinOp that downsamples a signal x of size N.

This operator performs decimation on provided signal, where the result contains elements selected every step elements in the original signal, starting at start up to stop, if provided.

When L = downsample(N, step, start, stop), computing y = L @ x is equivalent to y = x[start:stop:step].

Args:

N: int: Length of the input vector.
step: int: The stride to use for selecting elements in the original signal.
start: int, optional: Index of the element in the original signal where decimation starts. Default is 0.
stop: int, optional: Index of the element (excluded) in the original signal where decimation stops. Default is None, meaning that chunkins stops at the end of the signal.

Returns:

The downsample LazyLinOp.

Examples:

>>> import numpy as np
>>> import lazylinop as lz
>>> N = 10
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> # Downsample every 2 element, starting at the third
>>> lz.signal.downsample(N, 2, 3) @ x
array([3, 5, 7, 9])

lazylinop.signal.overlap_add(N, size, overlap=1)#

Returns a LazyLinOp L for overlap-add linear operator.

Such an operator L applies to an input signal x with N elements, where N must be a multiple of size, and is considered as the concatenation of \(N/size\) chunks of size elements.

The output y = L@x results from the addition of these chunks with overlap: The first overlap elements of a block \(i+1\) with last overlap elements of block \(i\). Blocks size is given by size.

Shape of L is \((M, N)\) with \(M = N - overlap imes ( N/size - 1)\)

Args:

N: int: Size of the input.
size: int: Block size. Must divide N.
overlap: int: Size of the overlap between blocks. Must be strictly positive, and smaller than the block size size.

Returns:

LazyLinOp

Overlap-add description#

Examples:

>>> from lazylinop.signal import overlap_add
>>> import numpy as np
>>> x = np.full(5, 1.0)
>>> # Do nothing because input size is equal to window
>>> L = overlap_add(5, 5, overlap=1)
>>> np.allclose(L @ x, x)
True
>>> x = np.full(6, 1.0)
>>> L = overlap_add(6, 2, overlap=1)
>>> L @ x
array([1., 2., 2., 1.])

lazylinop.signal.dwt_coeffs_sizes(N, wavelet='haar', level=None, mode='zero')#

Return a list of int that gives the size of the coefficients [cAn, cDn, ..., cD2, cD1].

Args:

N, wavelet, level, mode:: See dwt() for more details.

Returns:

list of int.

Examples:

>>> from lazylinop.signal.dwt import dwt_coeffs_sizes
>>> import pywt
>>> W = pywt.Wavelet('haar').dec_len
>>> dwt_coeffs_sizes(5, 'haar', level=2)
[2, 2, 3]

lazylinop.signal.dwt_to_pywt_coeffs(x, N, wavelet='haar', level=None, mode='zero')#

Returns Pywavelets compatible [cAn, cDn, ..., cD2, cD1] computed from x where n is the decomposition level.

Args:

x: np.ndarray: List of coefficients [cAn, cDn, ..., cD2, cD1].
N, wavelet, level, mode:: See dwt() for more details.

Returns:

Pywavelets compatible list [cAn, cDn, ..., cD2, cD1].

Examples:

>>> import numpy as np
>>> from lazylinop.signal import dwt, dwt_to_pywt_coeffs
>>> import pywt
>>> N = 11
>>> x = np.arange(N)
>>> L = dwt(N, wavelet='haar', level=2, mode='zero')
>>> y = dwt_to_pywt_coeffs(L @ x, N, 'haar', level=2, mode='zero')
>>> z = pywt.wavedec(x, wavelet='haar', level=2, mode='zero')
>>> np.allclose(y[0], z[0])
True
>>> np.allclose(y[1], z[1])
True

Signal module

Contents

Signal module#

List of transforms#

Padding or cropping#

Inverse transforms#

Other normalizations#

Transforms#

Utility functions#

Chunk description#

Overlap-add description#