Signal module

lazylinop.wip.signal.anti_eye(m, n=None, k=0, dtype='float')

Constructs a LazyLinOp whose equivalent array is filled with ones on the k-th antidiagonal and zeros everywhere else.

L = anti_eye(m, n, k) is such that L.toarray() == numpy.flip(numpy.eye(m, n, k), axis=1).

Args:

m: int

Number of rows.

n: int, optional

Number of columns (default is m).

k: int, optional

Anti-diagonal to place ones on.

zero for the main antidiagonal (default),
positive integer for an upper antidiagonal.
negative integer for a lower antidiagonal,

if k >= n or k <= - m then anti_eye(m, n, k) is zeros((m, n)) (k-th antidiagonal is out of operator shape).

dtype: str or numpy.dtype, optional

Defaultly 'float'.

Returns:

The anti-eye LazyLinOp.

Raises:

ValueError: “m and n must be >= 1.”

Examples:

>>> import lazylinop as lz
>>> import numpy as np
>>> x = np.arange(3)
>>> L = lz.wip.signal.anti_eye(3)
>>> np.allclose(L @ x, np.flip(x))
True
>>> L = lz.wip.signal.anti_eye(3, n=3, k=0)
>>> L.toarray()
array([[0., 0., 1.],
       [0., 1., 0.],
       [1., 0., 0.]])
>>> L = lz.wip.signal.anti_eye(3, n=3, k=1)
>>> L.toarray()
array([[0., 1., 0.],
       [1., 0., 0.],
       [0., 0., 0.]])
>>> L = lz.wip.signal.anti_eye(3, n=3, k=-1)
>>> L.toarray()
array([[0., 0., 0.],
       [0., 0., 1.],
       [0., 1., 0.]])
>>> L = lz.wip.signal.anti_eye(3, n=4, k=0)
>>> L.toarray()
array([[0., 0., 0., 1.],
       [0., 0., 1., 0.],
       [0., 1., 0., 0.]])
>>> L = lz.wip.signal.anti_eye(3, n=4, k=1)
>>> L.toarray()
array([[0., 0., 1., 0.],
       [0., 1., 0., 0.],
       [1., 0., 0., 0.]])
>>> L = lz.wip.signal.anti_eye(3, n=4, k=-1)
>>> L.toarray()
array([[0., 0., 0., 0.],
       [0., 0., 0., 1.],
       [0., 0., 1., 0.]])

lazylinop.wip.signal.bc(L, n=1, bn=0, an=0, boundary='periodic')

Builds a periodic or symmetric boundary condition LazyLinOp.

For an input \(x_1, x_2, \ldots, x_N\):

A symmetric boundary condition is such that:

\(x_N, ..., x_2, x_1 | x_1, x_2, ..., x_N | x_N, ..., x_2, x_1\)
A periodic boundary condition is such that:

\(x_1, x_2, ..., x_N | x_1, x_2, ..., x_N | x_1, x_2, ..., x_N\)

If the input is a 2d array, it works on each column and returns the resulting horizontal concatenation (a 2d-array).

Args:

L: int: Size of the input.
n: int, optional: Duplicate the signal this number of times on both side.
bn: int, optional: Add this number of elements before. It first adds n times the input to the left of the original input and then adds bn elements.
an: int, optional: Add this number of elements after. It first adds n times the input to the right of the original input and then adds an elements.
boundary: str, optional: 'wrap'/'periodic' (default) or 'symm'/'symmetric' boundary condition.

Returns:

LazyLinOp

Raises:

ValueError: L must be strictly positive.
ValueError: n must be >= 0.
ValueError: bn and an must be >= 0.
ValueError: bn and an must be <= L.
ValueError: boundary is either ‘periodic’ (‘wrap’) or ‘symmetric’ (‘symm’).

Examples:

>>> import lazylinop as lz
>>> import numpy as np
>>> N = 3
>>> x = np.arange(N).astype(np.float_)
>>> L = lz.wip.signal.bc(N, n=1, boundary='periodic')
>>> L @ x
array([0., 1., 2., 0., 1., 2., 0., 1., 2.])
>>> L = lz.wip.signal.bc(N, n=1, boundary='symmetric')
>>> L @ x
array([2., 1., 0., 0., 1., 2., 2., 1., 0.])
>>> X = np.array([[0., 0.], [1., 1.], [2., 2.]])
>>> X
array([[0., 0.],
       [1., 1.],
       [2., 2.]])
>>> L @ X
array([[2., 2.],
       [1., 1.],
       [0., 0.],
       [0., 0.],
       [1., 1.],
       [2., 2.],
       [2., 2.],
       [1., 1.],
       [0., 0.]])
>>> L = lz.wip.signal.bc(N, n=1, bn=1, boundary='periodic')
>>> L @ x
array([2., 0., 1., 2., 0., 1., 2., 0., 1., 2.])
>>> L = lz.wip.signal.bc(N, n=1, bn=2, an=1, boundary='symmetric')
>>> L @ x
array([1., 2., 2., 1., 0., 0., 1., 2., 2., 1., 0., 0.])

See also

vstack(), hstack(), bc2d()

lazylinop.wip.signal.bc2d(shape, x=1, y=1, b0=0, a0=0, b1=0, a1=0, boundary='periodic')

Constructs a periodic or symmetric boundary condition lazy linear operator. It will be applied to a flattened image. It basically add image on bottom, left, top and right side. Symmetric boundary condition is something like (on both axis): xN, …, x2, x1 | x1, x2, …, xN | xN, …, x2, x1 while a periodic boundary condition is something like (on both axis): x1, x2, …, xN | x1, x2, …, xN | x1, x2, …, xN

Args:

shape: tuple Shape of the image x: int, optional 2 * x + 1 signals along the first axis y: int, optional 2 * y + 1 signals along the second axis b0: int, optional Add b0 lines before (along x). a0: int, optional Add a0 lines after (along x). b1: int, optional Add b1 lines before (along y). a1: int, optional Add a1 lines after (along y). boundary: str, optional wrap/periodic (default) or symm/symmetric boundary condition

Returns:

LazyLinOp

Raises:

ValueError: x and y must be >= 0.
ValueError: shape expects tuple (R, C).
ValueError: b0, a0, b1 and a1 must be >= 0.
ValueError: b0 (resp. b1) must be < shape[0] (resp. shape[1]).
ValueError: a0 (resp. a1) must be < shape[0] (resp. shape[1]).
ValueError: boundary excepts ‘wrap’, ‘periodic’, ‘symm’ or ‘symmetric’.

Examples:

>>> from lazylinop.wip.signal import bc2d
>>> import numpy as np
>>> X = np.array([[0., 1.], [2., 3.]])
>>> X
array([[0., 1.],
       [2., 3.]])
>>> Op = bc2d(X.shape, x=1, y=1, boundary='periodic')
>>> (Op @ X.ravel()).reshape(2 * (2 * 1 + 1), 2 * (2 * 1 + 1))
array([[0., 1., 0., 1., 0., 1.],
       [2., 3., 2., 3., 2., 3.],
       [0., 1., 0., 1., 0., 1.],
       [2., 3., 2., 3., 2., 3.],
       [0., 1., 0., 1., 0., 1.],
       [2., 3., 2., 3., 2., 3.]])
>>> X = np.array([[0., 1.], [2., 3.], [4., 5.]])
>>> X
array([[0., 1.],
       [2., 3.],
       [4., 5.]])
>>> Op = bc2d(X.shape, x=1, y=0, a0=2, a1=1, boundary='symmetric')
>>> (Op @ X.ravel()).reshape(3 * (2 * 1 + 1) + 2, 2 * (2 * 0 + 1) + 1)
array([[4., 5., 5.],
       [2., 3., 3.],
       [0., 1., 1.],
       [0., 1., 1.],
       [2., 3., 3.],
       [4., 5., 5.],
       [4., 5., 5.],
       [2., 3., 3.],
       [0., 1., 1.],
       [0., 1., 1.],
       [2., 3., 3.]])

lazylinop.wip.signal.convolve(in1, in2, mode='full', method='scipy encapsulation', disable_jit=0)

Creates a lazy linear operator that corresponds to the convolution of a signal of size in1 with a kernel in2. If signal is a 2d array (in1, batch), return convolution per column. Kernel argument in2 must be 1d array even if input signal is a 2d array.

Args:

in1: int: Length of the input.
in2: np.ndarray: 1d kernel to convolve with the signal, shape is (K, ).

mode: str, optional

‘full’ computes convolution (input + padding).

‘valid’ computes ‘full’ mode and extract centered output that does not depend on the padding.

‘same’ computes ‘full’ mode and extract centered output that has the same shape that the input.

‘circ’ computes circular convolution

method: str, optional

‘auto’ use best implementation.

‘direct’ direct computation using nested for loops. It uses Numba jit. Larger the signal is better the performances are. Larger the batch size is better the performances are.

‘scipy encapsulation’ (default) to use lazy encapsulation of Scipy.signal convolve function. For large batch and small signal sizes, it uses encapsulation of sp.linalg.toeplitz.

‘pyfaust toeplitz’ to use pyfaust implementation of Toeplitz matrix.

‘oa’ to use lazylinop implementation of overlap-add method.

The following methods only work with mode='circ'.

‘scipy fft’ use Scipy implementation of FFT to compute circular convolution.

‘pyfaust circ’ use pyfaust implementation of circulant matrix.

‘pyfaust dft’ use pyfaust implementation of DFT.

disable_jit: int, optional: If 0 (default) disable Numba jit.

Returns:

LazyLinOp

Raises:

Exception Negative dimension of the signal and/or the kernel is not equal to 1. ValueError mode is either ‘full’ (default), ‘valid’, ‘same’ or ‘circ’. Exception in1 expects int. ValueError Size of the kernel is greater than the size of signal and mode is valid. ValueError method is not in: ‘auto’, ‘direct’, ‘scipy encapsulation’, ‘pyfaust toeplitz’, ‘oa’, ‘scipy fft’, ‘pyfaust circ’, ‘pyfaust dft’ ValueError method=’pyfaust circ’, ‘scipy fft’ or ‘pyfaust dft’ works only with mode=’circ’.

Examples:

>>> import numpy as np
>>> from lazylinop.wip.signal import convolve
>>> import scipy as sp
>>> x = np.random.rand(1024)
>>> kernel = np.random.rand(32)
>>> c1 = convolve(x.shape[0], kernel, mode='same', method='direct') @ x
>>> c2 = sp.signal.convolve(x, kernel, mode='same', method='auto')
>>> np.allclose(c1, c2)
True
>>> N = 32768
>>> x = np.random.rand(N)
>>> kernel = np.random.rand(48)
>>> c1 = convolve(N, kernel, mode='circ', method='scipy fft') @ x
>>> c2 = convolve(N, kernel, mode='circ', method='pyfaust dft') @ x
>>> np.allclose(c1, c2)
True

lazylinop.wip.signal.convolve2d(in1, in2, mode='full', boundary='fill', method='auto', **kwargs)

Constructs a 2d convolution lazy linear operator C. C @ X.flatten() is the result of the convolution. Of note, 2d-array X has been flattened. The output is also flattened and you have to reshape it according to the convolution mode. Toeplitz based method use the fact that convolution of a kernel with an image can be written as a sum of Kronecker product between eye and Toeplitz matrices.

Args:

in1: tuple, Shape (X, Y) of the signal to convolve with kernel. in2: np.ndarray Kernel to use for the convolution, shape is (K, L) mode: str, optional ‘full’ computes convolution (input + padding) ‘valid’ computes ‘full’ mode and extract centered output that does not depend on the padding. ‘same’ computes ‘full’ mode and extract centered output that has the same shape that the input. See also Scipy documentation of scipy.signal.convolve function for more details boundary: str, optional ‘fill’ pads input array with zeros (default) ‘wrap’ periodic boundary conditions ‘symm’ symmetrical boundary conditions method: str, optional ‘auto’ to use the best method according to the kernel and input array dimensions. ‘scipy encapsulation’ use scipy.signal.convolve2d as a lazy linear operator ‘scipy toeplitz’ to use lazy encapsulation of Scipy implementation of Toeplitz matrix. ‘pyfaust toeplitz’ to use FAuST implementation of Toeplitz matrix. ‘scipy fft’ to use Fast-Fourier-Transform to compute convolution.

Returns:

LazyLinOp or np.ndarray

Raises:

ValueError: mode is either ‘full’ (default), ‘valid’ or ‘same’.
ValueError: boundary is either ‘fill’ (default), ‘wrap’ or ‘symm’
Exception: Size of the kernel is greater than the size of signal (mode is valid).
ValueError: Unknown method.
Exception: in1 expects tuple as (X, Y).
Exception: in1 expects array with shape (X, Y).
ValueError: Negative dimension value is not allowed.

Examples:

>>> from lazylinop.wip.signal import convolve2d
>>> import scipy as sp
>>> X = np.random.randn(6, 6)
>>> kernel = np.random.randn(3, 3)
>>> c1 = convolve2d(X.shape, kernel, mode='same') @ X.flatten()
>>> c2 = sp.signal.convolve2d(X, kernel, mode='same').flatten()
>>> np.allclose(c1, c2)
True

lazylinop.wip.signal.dsconvolve(in1, in2, mode='full', offset=0, every=2, disable_jit=0)

Creates convolution plus down-sampling lazy linear operator. If input is a 2d array shape=(in1, batch), return convolution per column. offset (0 or 1) argument determines the first element to compute while every argument determines distance between two elements (1 or 2). The ouput of convolution followed by down-sampling C @ x is equivalent to scipy.signal.convolve(x, in2, mode)[offset::every].

Args:

in1: int: Length of the input.
in2: np.ndarray: 1d kernel to convolve with the signal, shape is (K, ).

mode: str, optional

‘full’ computes convolution (input + padding).

‘valid’ computes ‘full’ mode and extract centered output that does not depend on the padding.

‘same’ computes ‘full’ mode and extract centered output that has the same shape that the input.

offset: int, optional: First element to keep (default is 0).
every: int, optional: Keep element every this number (default is 2).
disable_jit: int, optional: If 0 (default) enable Numba jit.

Returns:

LazyLinOp

Raises:

ValueError: Number of dimensions of the kernel must be 1.
ValueError: Negative input length.
ValueError: mode is either ‘full’ (default), ‘valid’ or ‘same’
Exception: Length of the input are expected (int).
ValueError: Size of the kernel is greater than the size of signal and mode is valid.
Exception: mode, offset and every are incompatibles with kernel and signal sizes.
ValueError: offset must be either 0 or 1.
ValueError: every must be either 1 or 2.

Examples:

>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop.wip.signal import dsconvolve
>>> x = np.random.rand(1024)
>>> kernel = np.random.rand(32)
>>> Op = dsconvolve(x.shape, kernel, mode='same', offset=0, every=2)
>>> c1 = Op @ x
>>> c2 = sp.signal.convolve(x, kernel, mode='same', method='auto')
>>> np.allclose(c1, c2[0::2])
True

lazylinop.wip.signal.ds_mconv(in1, in2, in3, mode='full', offset=0, disable_jit=0)

Creates convolution plus down-sampling lazy linear operator. It first computes convolution with in2 and in3 filters. Then, it performs down-sampling (keep 1 every 2) on both convolution results (it is useful for Discrete-Wavelet-Transform). If input x is a 1d array, C @ x return concatenation of both convolution. If input X is a 2d array, C @ X return concatenation per column. offset (0 or 1) argument determines the first element to compute. The ouput C @ x is equivalent to the concatenation of scipy.signal.convolve(x, in2, mode)[offset::2] and scipy.signal.convolve(x, in3, mode)[offset::2].

Args:

in1: int: Length of the input.
in2: np.ndarray: First 1d kernel to convolve with the signal, shape is (K, ).
in3: np.ndarray: Second 1d kernel to convolve with the signal, shape is (K, ).

mode: str, optional

‘full’ computes convolution (input + padding).

‘valid’ computes ‘full’ mode and extract centered output that does not depend on the padding.

‘same’ computes ‘full’ mode and extract centered output that has the same shape that the input.

offset: int, optional: First element to keep (default is 0).
disable_jit: int, optional: If 0 (default) enable Numba jit.

Returns:

LazyLinOp

Raises:

ValueError: Number of dimensions of the kernel must be 1.
ValueError: Negative input length.
ValueError: mode is either ‘full’ (default), ‘valid’ or ‘same’
Exception: Length of the input are expected (int).
Exception: in2 and in3 must have the same length.
ValueError: Size of the kernel is greater than the size of signal and mode is valid.
Exception: mode and offset values are incompatibles with kernel and signal sizes.
ValueError: offset must be either 0 or 1.

Examples:

>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop.wip.signal import ds_conv_lh
>>> x = np.random.rand(1024)
>>> l = np.random.rand(32)
>>> h = np.random.rand(32)
>>> Op = ds_conv_lh(x.shape, l, h, mode='same', offset=0)
>>> c1 = Op @ x
>>> c2 = sp.signal.convolve(x, l, mode='same', method='auto')
>>> c3 = sp.signal.convolve(x, h, mode='same', method='auto')
>>> np.allclose(c1, np.hstack((c2[0::2], c3[0::2])))
True

lazylinop.wip.signal.dct(N, type=2, n=None, norm='backward', workers=None, orthogonalize=None, backend='scipy')

Returns a LazyLinOp for the DCT of size N.

If the input array is a batch of vectors, apply DCT per column. The function provides two backends. One is based on the DCT from SciPy while the other one builds DCT from lazylinop building blocks.

Args:

N: int

Size of the input (N > 0).

type: int, optional

Defaut is 2. See SciPy DCT for more details.

n: int, optional

Default is None. See SciPy DCT for more details.

norm: str, optional

Default is ‘backward’. See SciPy DCT for more details.

workers: int, optional

Default is None See SciPy DCT for more details.

orthogonalize: bool, optional

Default is None. See SciPy DCT for more details.

backend: str, optional

'scipy' (default) uses scipy.fft.dct to compute DCT.
'lazylinop' uses building-blocks to compute DCT.

lazylinop backend computation of L @ x is equivalent to sp.fft.dct(x, type, n, 0, norm, False, 1, orthogonalize).

Returns:

LazyLinOp

Raises:

Exception: DCT I: size of the input must be >= 2.
ValueError: norm must be either ‘backward’, ‘ortho’ or ‘forward’.
ValueError: type must be either 1, 2, 3 or 4.

Example:

>>> import lazylinop as lz
>>> import numpy as np
>>> import scipy as sp
>>> F = lz.wip.signal.dct(32)
>>> x = np.random.rand(32)
>>> y = sp.fft.dct(x)
>>> np.allclose(F @ x, y)
True

lazylinop.wip.signal.dst(N, type=2, n=None, norm='backward', workers=None, orthogonalize=None, backend='scipy')

Returns a LazyLinOp for the DST of size N.

If the input array is a batch of vectors, apply DST per column. The function provides two backends. One is based on the DST from SciPy while the other one builds DST from lazylinop building blocks.

Args:

N: int

Size of the input (N > 0).

type: int, optional

Defaut is 2. See SciPy DST for more details.

n: int, optional

Default is None. See SciPy DST for more details.

norm: str, optional

Default is ‘backward’. See SciPy DST for more details.

workers: int, optional

Default is None. It only works when backend='lazylinop'. See SciPy DST for more details.

orthogonalize: bool, optional

Default is None. See SciPy DST for more details.

backend: str, optional

'scipy' (default) uses scipy.fft.dst to compute DST.
'lazylinop' uses building-blocks to compute DST.

lazylinop backend computation of L @ x is equivalent to sp.fft.dst(x, type, n, 0, norm, False, 1, orthogonalize).

Returns:

LazyLinOp

Raises:

Exception: DST I: N and n must be > 1.
Exception: N and n must be >= 1.
ValueError: norm must be either ‘backward’, ‘ortho’ or ‘forward’.
ValueError: type must be either 1, 2, 3 or 4.

Example:

>>> import lazylinop as lz
>>> import numpy as np
>>> import scipy as sp
>>> F = lz.wip.signal.dst(32)
>>> x = np.random.rand(32)
>>> y = sp.fft.dst(x)
>>> np.allclose(F @ x, y)
True

See also

Wikipedia, SciPy DST, SciPy DSTn.

lazylinop.wip.signal.dwt1d(N, wavelet=pywt._extensions._pywt.Wavelet(name='haar', filter_bank=([0.7071067811865476, 0.7071067811865476], [-0.7071067811865476, 0.7071067811865476], [0.7071067811865476, 0.7071067811865476], [0.7071067811865476, -0.7071067811865476])), mode='zero', level=None, backend='pywavelets', **kwargs)

Constructs a Discrete Wavelet Transform (DWT) lazy linear operator. For the first level of decomposition Op @ x returns the array [cA, cD] while for the nth level of decomposition Op @ x returns the array [cAn, cDn, cDn-1, …, cD2, cD1]. The 1d array of coefficients corresponds to the concatenation of the list of arrays Pywavelets returns. Of note, the function follows the format returned by Pywavelets module.

Args:

N: int: Size of the input array.
wavelet: pywt.Wavelet: Wavelet from Pywavelets module.

mode: str, optional

‘periodic’, signal is treated as periodic signal.

‘symmetric’, use mirroring to pad the signal.

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

level: int, optional

Decomposition level, by default (None) return all.

kwargs:

lfilter: np.ndarray: Quadrature mirror low-pass filter.
hfilter: np.ndarray: Quadrature mirror high-pass filter. If lfilter and hfilter, ignore wavelet argument.

backend: str, optional

‘pywavelets’ (default) or ‘pyfaust’ for the underlying computation of the DWT.

Returns:

The DWT LazyLinOp.

Raises:

ValueError: Decomposition level must be >= 0.
ValueError: mode is either ‘periodic’, ‘symmetric’ or ‘zero’.
ValueError: level is greater than the maximum decomposition level.
ValueError: backend must be either ‘pywavelets’ or ‘lazylinop’.
ValueError: backend ‘pywavelets’ works only for mode=’zero’. We are currently implementing adjoint for others modes.

Examples:

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

lazylinop.wip.signal.dwt2d(shape, wavelet=pywt._extensions._pywt.Wavelet(name='haar', filter_bank=([0.7071067811865476, 0.7071067811865476], [-0.7071067811865476, 0.7071067811865476], [0.7071067811865476, 0.7071067811865476], [0.7071067811865476, -0.7071067811865476])), mode='zero', level=None, backend='pywavelets', **kwargs)

Constructs a multiple levels 2d DWT lazy linear operator. If the lazy linear operator is applied to a 1d array it returns the array [cA, cH, cV, cD] for the first decomposition level. For the nth level of decomposition it returns the array [cAn, cHn, cVn, cDn, …, cH1, cV1, cD1]. cAn matrix has been flattened as-well-as the cHs, cVs and cDs. The 1d array of coefficients corresponds to the concatenation of the list of arrays Pywavelets returns. Of note, you must apply DWT LazyLinOp to a flattened 2d array.

Args:

shape: tuple: Shape of the 2d input array (M, N).
wavelet: pywt.Wavelet: Wavelet from Pywavelets module.

mode: str, optional

‘periodic’, image is treated as periodic image

‘symmetric’, use mirroring to pad the signal

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

level: int, optional: If level is None compute full decomposition (default).
backend: str, optional: ‘pywavelets’ (default) or ‘pyfaust’ for the underlying computation of the DWT.
kwargs:: lfilter: np.ndarray Quadrature mirror low-pass filter. hfilter: np.ndarray Quadrature mirror high-pass filter. If lfilter and hfilter and backend=’lazylinop’, ignore wavelet argument.

Returns:

LazyLinOp

Raises:

Exception: Shape expects tuple (M, N).
ValueError: Decomposition level must be >= 0.
Exception: Decomposition level is greater than the maximum decomposition level.
ValueError: mode is either ‘zero’, ‘periodic’ or ‘symmetric’.
ValueError: backend must be either ‘pywavelets’ or ‘lazylinop’.
ValueError: backend ‘pywavelets’ works only for mode=’zero’. We are currently implementing adjoint for others modes.
Exception: Length of the wavelet must be > 1.

Examples:

>>> from lazylinop.wip.signal import dwt2d
>>> import numpy as np
>>> import pywt
>>> X = np.array([[1., 2.], [3., 4.]])
>>> Op = dwt2d(X.shape, wavelet=pywt.Wavelet('haar'), level=1)
>>> y = Op @ X.flatten()
>>> cA, (cH, cV, cD) = pywt.wavedec2(X, wavelet='haar', level=1)
>>> z = np.concatenate([cA, cH, cV, cD], axis=1)
>>> np.allclose(y, z)
True

lazylinop.wip.signal.fft(N, backend='scipy', disable_jit=0, **kwargs)

Returns a LazyLinOp for the DFT of size N.

Args:

N: int: Size of the input (N > 0).
backend:: ‘scipy’ (default) or ‘pyfaust’ for the underlying computation of the DFT. If we denote F the LazyLinOp DFT, X a batch of vectors and F @ X the DFT of X, backend ‘scipy’ uses fft (resp. fftn) encapsulation when batch size is 1 (resp. > 1). ‘iterative’ backend use Numba prange function and shows good performances when the size of the signal is greater than 1e6 and the size of the batch of vectors is greater than 1000. We do not recommend to use it for small size.
kwargs:: Any key-value pair arguments to pass to the SciPy or pyfaust dft backend.

Returns:

LazyLinOp

Raises:

ValueError: Invalid norm value for ‘scipy’ backend.
ValueError: backend must be either ‘scipy’ or ‘pyfaust’.

Example:

>>> from lazylinop.wip.signal import fft
>>> import numpy as np
>>> F1 = fft(32, norm='ortho')
>>> F2 = fft(32, backend='pyfaust')
>>> x = np.random.rand(32)
>>> np.allclose(F1 @ x, F2 @ x)
True
>>> y = F1 @ x
>>> np.allclose(F1.H @ y, x)
True
>>> np.allclose(F2.H @ y, x)
True

lazylinop.wip.signal.fft2(shape, backend='scipy', **kwargs)

Returns a LazyLinOp for the 2D DFT of size n.

Args:

shape: tuple: The signal shape to apply the fft2 to.
backend: str, optional: ‘scipy’ (default) or ‘pyfaust’ for the underlying computation of the 2D DFT. pyfaust allows power of two only.
kwargs:: Any key-value pair arguments to pass to the scipy or pyfaust dft backend.

Returns:

LazyLinOp

Raises:

ValueError: backend must be either scipy or pyfaust.

Example:

>>> from lazylinop.wip.signal import fft2
>>> import numpy as np
>>> F_scipy = fft2((32, 32), norm='ortho')
>>> F_pyfaust = fft2((32, 32), backend='pyfaust')
>>> x = np.random.rand(32, 32)
>>> np.allclose(F_scipy @ x.ravel(), F_pyfaust @ x.ravel())
True
>>> y = F_scipy @ x.ravel()
>>> np.allclose(F_scipy.H @ y, x.ravel())
True
>>> np.allclose(F_pyfaust.H @ y, x.ravel())
True

See also

SciPy fft2, pyfaust dft

lazylinop.wip.signal.flip(shape, start=0, end=None, axis=0)

Constructs a flip lazy linear operator.

Args:

shape: tuple shape of the input start: int, optional flip from start (default is 0) end: int, optional stop flip (not included, default is None) axis: int, optional if axis=0 (default) flip per column, if axis=1 flip per row it does not apply if shape[1] is None.

Returns:

The flip LazyLinOp

Raises:

ValueError: start is < 0.
ValueError: start is >= number of elements along axis.
ValueError: end is < 1.
ValueError: end is > number of elements along axis.
Exception: end is <= start.

Examples:

>>> import numpy as np
>>> from lazylinop.wip.signal import flip
>>> x = np.arange(6)
>>> x
array([0, 1, 2, 3, 4, 5])
>>> y = flip(x.shape, 0, 5) @ x
>>> y
array([4, 3, 2, 1, 0, 5])
>>> z = flip(x.shape, 2, 4) @ x
>>> z
array([0, 1, 3, 2, 4, 5])
>>> X = np.eye(5, M=5, k=0)
>>> X
array([[1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 1., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 1.]])
>>> flip(X.shape, 1, 4) @ X
array([[1., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 1., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 1.]])

lazylinop.wip.signal.fwht(N, backend='pyfaust', disable_jit=0)

Creates a Fast-Walsh-Hadamard-Transform lazy linear operator. The size of the signal has to be a power of two greater than 1.

Args:

N: int

Size the input array.

backend: str, optional

It can be ‘direct’, ‘kronecker’, ‘pyfaust’ (default), ‘pytorch’ (wip) or ‘scipy’.

‘kronecker’ uses kron to compute FWHT.
‘direct’ uses N * log(N) algorithm and Numba jit to compute FWHT.
‘pytorch’ does not work yet.
‘scipy’ uses scipy.linalg.hadamard.
‘pyfaust’ uses wht from pyfaust.
‘auto’ uses ‘pyfaust’ for \(N<=32768\), ‘direct’ otherwise. Of note, this condition depends on your hardware. It may change in the near future.

disable_jit: int, optional

If 0 (default) enable Numba jit.

Returns:

LazyLinOp

Raises:

ValueError: The size of the signal must be a power of two, greater or equal to two.
ValueError: backend argument expects either ‘direct’, ‘kronecker’, ‘pytorch’ (do not work yet), pyfaust or ‘scipy’.

Examples:

>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop.wip.signal import fwht
>>> x = np.random.randn(16)
>>> y = fwht(x.shape[0]) @ x
>>> np.allclose(y, sp.linalg.hadamard(x.shape[0]) @ x)
True
>>> X = np.random.randn(8, 4)
>>> Y = fwht(X.shape[0]) @ X
>>> np.allclose(Y, sp.linalg.hadamard(X.shape[0]) @ X)
True

lazylinop.wip.signal.slices(N, start, end, window=None, nhop=None, offset=(0, None))

Constructs a multiple slices lazy linear operator.

If window and hop are None (default) extract slices given by the intervals [start[i], end[i]] for all i < len(start). Element start[i] must be lesser than element end[i]. Element end[i] must be greater or equal than element start[i - 1]. If start[i] = end[i], extract only one element.
If window and nhop are not None extract slices given by window and nhop arguments. For a given signal, window size (window) and number of elements (nhop) between two consecutive windows the lazy linear operator concatenates the windows into a signal that could be larger than the original one. The number of windows is given by 1 + (signal length - window) // nhop. Therefore, the length of the output signal is nwindows * window.

If batch greater than 1, apply to each column.

Args:

N: int: Length of the input.
start: np.ndarray: List of first element to keep.
end: np.ndarray: List of last element to keep.
window: int, optional: Size of the window.
nhop: int, optional: Number of elements between two windows.
offset: tuple, optional: First element to keep, default is 0. Last element to keep, default is None.

Returns:

The multiple slices LazyLinOp

Raises:

Exception: start and end do not have the same length.
ValueError: start must be positive.
ValueError: end must be strictly positive.
Exception: end must be >= start.
Exception: start must be < N.
Exception: end must be < N.
Exception: start is empty.
ValueError: window argument expects a value > 0 and <= signal length.
ValueError: nhop argument expects a value > 0.
ValueError: offset[0] is < 0.
ValueError: offset[0] >= number of elements.
ValueError: offset[1] > number of elements.
Exception: offset[1] is <= offset[0].

Examples:

>>> import numpy as np
>>> import lazylinop as lz
>>> # use list to extract multiple slices
>>> x = np.arange(10)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> lz.wip.signal.slices(x.shape[0], [0, 5], [2, 8]) @ x
array([0, 1, 2, 5, 6, 7, 8])
>>> # use window = 1 and nhop
>>> x = np.arange(10)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> lz.wip.signal.slices(x.shape[0], [], [], 1, 2, (0, 10)) @ x
array([0, 2, 4, 6, 8])
>>> X = np.arange(30).reshape((10, 3))
>>> X
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14],
       [15, 16, 17],
       [18, 19, 20],
       [21, 22, 23],
       [24, 25, 26],
       [27, 28, 29]])
>>> lz.wip.signal.slices(x.shape[0], [], [], 1, 2, (0, 10)) @ X
array([[ 0,  1,  2],
       [ 6,  7,  8],
       [12, 13, 14],
       [18, 19, 20],
       [24, 25, 26]])
>>> # use window > 1 and nhop
>>> x = np.array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])
>>> lz.wip.signal.slices(x.shape[0], [], [], 5, 2) @ x
array([0., 1., 2., 3., 4., 2., 3., 4., 5., 6., 4., 5., 6., 7., 8.])
>>> X = x.reshape(5, 2)
>>> X
array([[0., 1.],
       [2., 3.],
       [4., 5.],
       [6., 7.],
       [8., 9.]])
>>> lz.wip.signal.slices(X.shape[0], [], [], 3, 2) @ X
array([[0., 1.],
       [2., 3.],
       [4., 5.],
       [4., 5.],
       [6., 7.],
       [8., 9.]])

lazylinop.wip.signal.oa(L, X, overlap=1)

return overlap-add linear operator. The overlap-add linear operator adds last overlap of block i > 0 with first overlap of block i + 1. Of note, block i = 0 (of size L - overlap) does not change.

Args:

L: int Block size. X: int Number of blocks. overlap: int Size of the overlap < L (strictly positive).

Returns:

LazyLinOp or np.ndarray

Raises:

ValueError: L is strictly positive.
ValueError: X is strictly positive.
ValueError: overlap must be > 0 and <= L

Examples:

>>> from lazylinop.wip.signal import oa
>>> import numpy as np
>>> signal = np.full(5, 1.0)
>>> oa(1, 5, overlap=1) @ signal
array([5.])
>>> signal = np.full(10, 1.0)
>>> oa(2, 5, overlap=1) @ signal
array([1., 2., 2., 2., 2., 1.])

lazylinop.wip.signal.stft(N, fs=1.0, window='hann', nperseg=256, noverlap=None, boundary='zeros', padded=True, scaling='spectrum')

Constructs a Short-Time-Fourier-Transform lazy linear operator. The function uses basic operators to build STFT operator.

Args:

N: int: Length of the input array.
fs: int, optional: Sampling frequency (1 is default).
window: str, optional: Window name to use to avoid discontinuity if the segment was looped indefinitly (‘hann’ is default). See scipy.signal.get_window for a list of available windows.
nperseg: int, optional: Number of samples in a frame (256 is default).
noverlap: int, optional: Number of samples to overlap between two consecutive segments (None is default correspoding to nperseg // 2).
boundary: str or None, optional: How to extend signal at both ends (‘zeros’ is default). Only option is ‘zeros’, others are WiP.
padded: bool, optional: Zero-pad the signal such that new length fits exactly into an integer number of window segments (True is default).
scaling: str, optional: Scaling mode (‘spectrum’ is default) follows scipy.signal.stft function, other possible choice is ‘psd’.

Returns:

LazyLinOp

Raises:

ValueError: window argument expects ‘hann’.
ValueError: scaling argument expects ‘spectrum’ or ‘psd’.
ValueError: noverlap expects value less than nperseg.
ValueError: nperseg expects value greater than 0.

Examples: