Signal module

Preambule

This module provides (mostly) orthonormal lazy linear operators associated to fast linear transforms commonly used in signal processing. This includes the Discrete Fourier Transform (dft), Discrete Wavelet Transform (dwt), Discrete Cosine/Sine Transforms (dct/dst) of types I to IV, the Modified DCT (mdct), the Walsh-Hadamard Transform (wht) and the Zak Transform (dzt). We also provide a lazy linear operator for convolution with a given filter.

Convolution:

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

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

Inverse transforms: you may be surprised that we are not providing 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 all the transforms we provide (to the exception of the MDCT 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 is the MDCT: 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.

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.

Warning

Do not confuse padder with pad (see Construction for more details).

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.

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

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 = N - filter.size + 1 (filter.size <= N must be satisfied).
'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

'scipy_convolve': (default) encapsulate scipy.signal.convolve.

It uses internally the best SciPy backend between 'fft' and 'direct' (see scipy.signal.choose_conv_method).
'auto': use best backend available.
- for mode='circ' use backend='fft'.
- for any other mode use backend='direct' if filter.size < log2(M), backend='scipy_convolve' otherwise.
'direct': direct computation using nested for-loops with Numba and parallelization.
'fft': (only for mode='circ') compute circular convolution using SciPy FFT.
'toeplitz': encapsulate scipy.linalg.toeplitz if N < 2048, scipy.linalg.matmul_toeplitz otherwise.
'oa': use Lazylinop implementation of overlap-add backend.

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 for more details.

backend: str, optional

'scipy' (default) uses scipy.fft.dct to compute 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 for more details.

backend: str, optional

'scipy' (default) uses scipy.fft.dst to compute the DST.
'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, backend='scipy')

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

Shape of L is \((H,~N)\) with \(H=N/2\). \(N\) must be multiple of \(4\).

X = L @ x yiels a vector X of size \(H\) 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}\]

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

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

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}\]

Args:

N: int

Size of the signal. N must be even.

backend: str, optional

'scipy' (default) uses scipy.fft.dct encapsulation for the underlying computation of the DCT.
'lazylinop' uses pre-built Lazylinop operators (Lazylinop fft(), eye(), vstack() etc.) to build the pipeline that will compute the MDCT.

Returns:

LazyLinOp

Example:

>>> from lazylinop.signal import mdct
>>> import numpy as np
>>> x = np.random.randn(64)
>>> L = mdct(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.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.

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 and decomposition level. In general, L is not orthogonal.

Args:

N: int: Size of the input array.

wavelet: str or tuple of (np.ndarray, np.ndarray), optional

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

If a tuple of two np.ndarray is provided, the low and high-pass filters used to define the wavelet. The dwt() function does not test Quadrature-Mirror-Filters of the custom wavelet.

mode: str, optional

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

‘periodic’, signal is treated as periodic signal. Only works with backend='lazylinop'.

‘reflect’, signal is extended by reflecting elements.

‘symmetric’, use mirroring to pad the signal. Only works with backend='lazylinop'.

‘zero’, signal is padded with zeros (default).

level: int, optional: Decomposition level, by default (None) the maximum level is used
backend: str, optional: ‘pywavelets’ (default) or ‘pyfaust’ for the underlying computation of the DWT.

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, mode='periodic', level=1, backend='lazylinop')
>>> y = L @ x
>>> z = pywt.wavedec(x, wavelet='haar', mode='periodic', level=1)
>>> np.allclose(y, np.concatenate(z))
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.

Utils

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