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(in1s, in2, mode='full', method='scipy_convolve', workers=None)

Builds a LazyLinOp for the 1d convolution of signal(s) of size in1s with a kernel in2.

See below, in1s and in2 for input sizes. See mode for output size.

Args:
in1s: int

Size of the first input signal(s).

The input shape can be (in1s,) for a single signal and (in1s, b) for a batch of b signals (to compute one convolution per column of size in1s with 1d-array in2).

See workers for parallelization of the 2d-array case.

in2: np.ndarray

Second input. Kernel to convolve with the signal. in2 is always a 1d-array.

mode: str, optional
  • 'full': compute full convolution including at points of non-complete overlapping of inputs. Output size is in1s + in2.size - 1.

  • 'valid': compute only fully overlapping part of the convolution. This is the ‘full’ center part of (output) size ins1 - in2.size + 1 (in2.size <= in1s must be satisfied).

  • 'same': compute only the center part of 'full' to obtain an output size equal to the input size in1s.

  • 'circ': compute circular convolution (in2.size <= in1s must be satisfied).

method: str, optional
  • 'scipy_convolve': (default) encapsulate scipy.signal.convolve.

    It uses internally best SciPy method between 'fft' and 'direct' (see scipy.signal.choose_conv_method).

  • 'auto': use best method available.

    • for mode='circ' use method='fft'.

    • for any other mode use method='direct' if in2.size < log2(output_size), method='scipy_convolve' otherwise.

    Note that, according to their performance, method='toeplitz' and method='oa' are never considered.

  • 'direct': direct computation using nested for-loops with Numba and parallelization (see workers).

  • 'toeplitz': encapsulate scipy.linalg.toeplitz if in1s < 2048, scipy.linalg.matmul_toeplitz otherwise.

  • 'oa': use Lazylinop implementation of overlap-add method.

  • 'fft': (only for mode='circ') compute circular convolution using SciPy FFT.

workers: int, optional

The number of threads used to parallelize method='direct', 'toeplitz', 'scipy_convolve' using respectively Numba, NumPy and SciPy capabilities. Default is os.cpu_count() (number of CPU threads available).

Environment override

workers can be overridden from the environment using NUMBA_NUM_THREADS for method='direct' and OMP_NUM_THREADS for method='toeplitz', 'scipy_convolve'.

Returns:

LazyLinOp

Raises:
TypeError

in1s must be an int.

ValueError

size in1s < 0

ValueError

mode is not valid (‘full’ (default), ‘valid’, ‘same’ or ‘circ’).

ValueError

in2.size > in1s and mode=’valid’

ValueError

in2.size > in1s and mode=’circ’

Exception

in2 must be 1d array.

ValueError

method is not in: ‘auto’, ‘direct’, ‘toeplitz’, ‘scipy_convolve’, ‘oa’, ‘fft’,

ValueError

method=’fft’ works only with mode=’circ’.

Examples:
>>> import numpy as np
>>> from lazylinop.wip.signal import convolve
>>> import scipy as sp
>>> x = np.random.randn(1024)
>>> kernel = np.random.randn(32)
>>> c1 = convolve(x.shape[0], kernel, method='direct') @ x
>>> c2 = sp.signal.convolve(x, kernel, method='auto')
>>> np.allclose(c1, c2)
True
>>> N = 32768
>>> x = np.random.randn(N)
>>> kernel = np.random.randn(48)
>>> c1 = convolve(N, kernel, mode='circ', method='fft') @ x
>>> c2 = convolve(N, kernel, mode='circ', method='direct') @ 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='ortho', workers=None, orthogonalize=None, backend='scipy')

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

Operator dimensions: \(n \times N\) (or \(N \times N\) if n=None).

If the input array is 2d, one DCT per column is applied (see workers for parallelization).

The function provides two backends: SciPy and Lazylinop.

To compute the inverse DCT, simply use dct(...).inv() (see example below). It works for any norm and orthogonalize configuration. For more details about the precise calculation you can consult Computing the inverse DCT from the transpose operator.

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.

n: int, optional

Length of the output / size of the DCT.

  • If n < N, crop the input before DCT.

  • If n > N, zero-pad the input before DCT.

  • If n=None or equivalent n=N, no cropping/padding.

norm: str, optional

Normalization mode:

  • norm='ortho' (default): normalization with scale factor \(1/\sqrt{2n}\) or \(1/\sqrt{2(n-1)}\) if type=1.

  • norm=None: no normalization.

  • norm='1/(2n)': the operator is scaled by \(1/(2n)\) or \(1/(2(n-1))\) if type=1.

See below the table for conversion to SciPy DCT norm argument.

The operator is orthogonal iff orthogonalize=True, norm='ortho' and \(N = n\) (square operator). See orthogonalize doc for more details.

workers: int, optional

Number of workers (default None means os.cpu_count()) to use for parallel computation if input has multiple columns. Works only for backend='scipy'.

See SciPy DCT for more details.

orthogonalize: bool, optional

Orthogonalization scaling.

Default is None (equivalent to True if norm='ortho', False otherwise).

  • orthogonalize=True:

    • type I: the inputs x[0], x[n-1] are scaled by \(\sqrt{2}\). Output y[0], y[n-1] are scaled by \({1 / \sqrt{2}}\).

    • type II: output y[0] is scaled by \({1 / \sqrt{2}}\).

    • type III: input x[0] is scaled by \({\sqrt{2}}\).

  • orthogonalize=False: no additional scaling of any input/output entry.

For type IV: no additional scaling whatever is orthogonalize.

See norm argument and SciPy DCT for more details.

backend: str, optional
  • 'scipy' (default) uses scipy.fft.dct to compute the DCT.

  • 'lazylinop' uses building-blocks to compute the DCT (Lazylinop fft(), vstack() etc.).

Returns:

LazyLinOp

Raises:
Exception

DCT I: size of the input must be >= 2.

ValueError

norm must be either ‘ortho’, ‘1/(2n)’ or None.

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, norm='ortho')
>>> np.allclose(F @ x, y)
True
>>> # compute the inverse DCT
>>> x_ = F.inv() @ y
>>> np.allclose(x_, x)
True

Correspondence to scipy.fft.dct/idct norm argument:

Forward operator

Lazylinop

SciPy

dct(N, norm=None) @ x

dct(x, norm='backward')

dct(N, norm='ortho') @ x

dct(x, norm='ortho')

dct(N, norm='1/(2n)') @ x

dct(x, norm='forward')

Transpose operator

Lazylinop

SciPy

dct(N, norm=None).T @ y

idct(y, norm='forward')

dct(N, norm='ortho').T @ y

idct(y, norm='ortho')

dct(N, norm='1/(2n)').T @ y

idct(y, norm='backward')

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

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

Operator dimensions: \(n \times N\) (or \(N \times N\) if n=None).

If the input array is 2d, one DST per column is applied (see workers for parallelization).

The function provides two backends: SciPy and Lazylinop.

To compute the inverse DST, simply use dst(...).inv() (see example below). It works for any norm and orthogonalize configuration. For more details about the precise calculation you can consult Computing the inverse DST from the transpose operator.

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.

n: int, optional

Length of the output / size of the DCT.

  • If n < N, crop the input before DCT.

  • If n > N, zero-pad the input before DCT.

  • If n=None or equivalent n=N, no cropping/padding.

norm: str, optional

Normalization mode:

  • norm='ortho' (default): normalization with scale factor \(1/\sqrt{2n}\) or \(1/\sqrt{2(n+1)}\) if type=1.

  • norm=None: no normalization.

  • norm='1/(2n)': the operator is scaled by \(1/(2n)\) or \(1/(2(n+1))\) if type=1.

See below the table for conversion to SciPy DST norm argument.

The operator is orthogonal iff orthogonalize=True, norm='ortho' and \(N = n\) (square operator). See orthogonalize doc for more details.

workers: int, optional

Number of workers (default None means os.cpu_count()) to use for parallel computation if input has multiple columns. Works only for backend='scipy'.

See SciPy DST for more details.

orthogonalize: bool, optional

Orthogonalization scaling.

Default is None (equivalent to True if norm='ortho', False otherwise).

  • orthogonalize=True:

    • type II: output y[n-1] is scaled by \({1 / \sqrt{2}}\).

    • type III: input x[n-1] is scaled by \({\sqrt{2} / 2}\).

  • orthogonalize=False: no additional scaling of any input/output entry.

For types I and IV: no additional scaling whatever is orthogonalize.

See norm argument and SciPy DST for more details.

backend: str, optional
  • 'scipy' (default) uses scipy.fft.dst to compute the DST.

  • 'lazylinop' uses building-blocks to compute the DST (Lazylinop fft(), vstack() etc.).

Returns:

LazyLinOp

Raises:
Exception

DST I: N and n must be > 1.

Exception

N and n must be >= 1.

ValueError

norm must be either ‘ortho’, ‘1/(2n)’ or None.

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, norm='ortho')
>>> np.allclose(F @ x, y)
True
>>> # compute the inverse DST
>>> x_ = F.inv() @ y
>>> np.allclose(x_, x)
True

Correspondence to scipy.fft.dst/idst norm argument:

Forward operator

Lazylinop

SciPy

dst(N, norm=None) @ x

dst(x, norm='backward')

dst(N, norm='ortho') @ x

dst(x, norm='ortho')

dst(N, norm='1/(2n)') @ x

dst(x, norm='forward')

Transpose operator

Lazylinop

SciPy

dst(N, norm=None).T @ y

idst(y, norm='forward')

dst(N, norm='ortho').T @ y

idst(y, norm='ortho')

dst(N, norm='1/(2n)').T @ y

idst(y, norm='backward')

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, n=None, norm='ortho', workers=None)

Builds a Discrete Fourier Transform (DFT) LazyLinOp.

Operator dimensions: \(n \times N\) (or \(N \times N\) if n=None).

SciPy FFT is used as underlying implementation.

To compute the inverse FFT, simply use fft(...).inv() (see example below). It works for any norm.

Args:
N: int

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

n: int, optional

Length of the output / size of the DFT.

  • If n < N, crop the input before DFT.

  • If n > N, zero-pad the input before DFT.

  • If n=None or equivalent n=N, no cropping/padding.

norm: str, optional

Normalization mode:

  • norm='ortho' (default): normalization with scale factor \(1/\sqrt{n}\). if square (\(N = n\)) the operator is unitary.

    Inverse FFT is the adjoint.

  • norm=None: no normalization. The operator is not unitary.

    Inverse FFT is the adjoint scaled by 1/n.

  • norm='1/n': the operator is scaled by \(1/n\). The operator is not unitary.

    Inverse FFT is the adjoint scaled by n.

See below the table for conversion to SciPy norm argument.

workers: int, optional

Number of workers (default is os.cpu_count()) to use for parallel computation.

See scipy.fft.fft for more details.

Returns:

LazyLinOp DFT

Raises:
ValueError

norm must be either ‘ortho’, ‘1/n’ or None.

Example:
>>> import lazylinop.wip.signal as lz
>>> import numpy as np
>>> import scipy as sp
>>> F = lz.fft(32)
>>> x = np.random.rand(32)
>>> np.allclose(F @ x, sp.fft.fft(x, norm='ortho'))
True
>>> # easy inverse
>>> F = lz.fft(32, norm=None)
>>> y = F @ x
>>> x_ = F.inv() @ y # inverse works for any norm
>>> np.allclose(x_, x)
True

Correspondence to scipy.fft.fft/ifft norm argument:

Forward operator

Lazylinop

SciPy

fft(N, norm=None) @ x

fft(x, norm='backward')

fft(N, norm='ortho') @ x

fft(x, norm='ortho')

fft(N, norm='1/n') @ x

fft(x, norm='forward')

Adjoint operator

Lazylinop

SciPy

fft(N, norm=None).H @ y

ifft(y, norm='forward')

fft(N, norm='ortho').H @ y

ifft(y, norm='ortho')

fft(N, norm='1/n').H  @ y

ifft(y, norm='backward')

lazylinop.wip.signal.rfft(N, n=None, norm='ortho', workers=None, fft_output=True)

Builds a Discrete Fourier Transform (DFT) LazyLinOp for real input.

Operator dimensions:
  • fft_output=True: \(n \times N\).

  • fft_output=False: \(L \times N\) with \(L = (n + 1) / 2\) if \(n\) is odd, \(L = (n / 2) + 1\) otherwise (take \(n = N\) if \(n\) is None).

SciPy real FFT is used as underlying implementation.

To compute the inverse real FFT, simply use rfft(...).inv() (see example below). It works for any norm.

Args:
N: int

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

n: int, optional

Crop/zero-pad the input to get a signal of size n to apply the DFT on. None (default) means n=N.

norm: str, optional

Normalization mode: 'ortho' (default), None or '1/n'. See fft() for more details.

workers: int, optional

Number of workers (default is os.cpu_count()) to use for parallel computation.

See scipy.fft.rfft for more details.

fft_output: bool, optional
  • True to get same output as fft (default).

  • False to get truncated output (faster but check() fails on forward - adjoint operators consistency).

Returns:

LazyLinOp real DFT

Raises:
ValueError

norm must be either ‘ortho’, ‘1/n’ or None.

Example:
>>> import lazylinop.wip.signal as lz
>>> import numpy as np
>>> import scipy as sp
>>> F = lz.rfft(32)
>>> x = np.random.rand(32)
>>> np.allclose(F @ x, sp.fft.fft(x, norm='ortho'))
True
>>> # easy inverse
>>> F = lz.rfft(32, norm=None)
>>> y = F @ x
>>> x_ = F.inv() @ y # inverse works for any norm
>>> np.allclose(x_, 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
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:

Signal processing appendices

Computing the inverse DST from the transpose operator

Computing inverse DST from DST of x depending on norm & orthogonalize

orthogonalize=True

DST

Inverse DST

L = dst(N, norm=norm_, type=1-4)

s * L.T

Scale factor s

norm_='ortho'

s=1

norm_=None and type=2-4

s = 1 / (2 * n)

norm_=None and type=1

s = 1 / (2 * (n+1))

norm_=1/(2n) and type=2-4

s = 2 * n

norm_=1/(2n) and type=1

s = 2 * (n+1)

orthogonalize=False

DST

Inverse DST

L = dst(N, norm=norm_, type=2, orthogonalize=False)

s * L.T @ I

I = eye(N);I[-1,-1]= 1 / 2

L = dst(N, norm=norm_, type=3, orthogonalize=False)

s * I @ L.T

I = eye(N);I[-1,-1] = 2

L = dst(N, norm=norm_, type=1|4, orthogonalize=False)

s * L.T

Computing the inverse DCT from the transpose operator

Computing inverse DCT from DCT of x depending on norm & orthogonalize

orthogonalize=True

DCT

Inverse DCT

L = dct(N, norm=norm_, type=1-4)

s * L.T

Scale factor s

norm_='ortho'

s=1

norm_=None and type=2-4

s = 1 / (2 * n)

norm_=None and type=1

s = 1 / (2 * (n-1))

norm_=1/(2n) and type=2-4

s = 2 * n

norm_=1/(2n) and type=1

s = 2 * (n-1)

orthogonalize=False

DCT

Inverse DCT

L = dct(N, norm=norm_, type=1, orthogonalize=False)

s * I2 @ L.T @ I1

I1 = eye(N)

I1[0,0]=.5;I1[-1,-1]=2

I2 = eye(N);I2[0,0]=2

L = dct(N, norm=norm_, type=2, orthogonalize=False)

s * L.T @ I

I = eye(N);I[0,0]= .5

L = dct(N, norm=norm_, type=3, orthogonalize=False)

s * I @ L.T

I = eye(N);I[0,0]= 2

L = dct(N, norm=norm_, type=4, orthogonalize=False)

s * L.T