Signal module
-------------

.. automodule:: lazylinop.signal

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

..
  This includes the Discrete Fourier Transform (:py:func:`dft`), Discrete Wavelet Transform (:py:func:`dwt`), Discrete Cosine/Sine Transforms (:py:func:`dct`/:py:func:`dst`) of types I to IV, the Modified DCT (:py:func:`mdct`) and its inverse MDCT (:py:func:`imdct`), the Walsh-Hadamard Transform (:py:func:`wht`) and the Zak Transform (:py:func:`dzt`).
..
  We also provide a lazy linear operator for **convolution with a given filter** (:py:func:`convolve`), a Short-Time-Fourier-Transform (:py:func:`stft`) with a given window and its inverse STFT (:py:func:`istft`).

List of transforms
~~~~~~~~~~~~~~~~~~

1. :meth:`lazylinop.signal.convolve`
   The lazy convolution operator is obtained as ``L = convolve(N, filter, ...)`` with ``N`` the dimension of the input signal and ``filter`` the filter (given as a vector), so that ``y = L @ x`` yields the prescribed convolution. Whenever needed ``z = L.H @ y`` computes the adjoint of the convolution operator, appropriately taking into account all boundary conditions that you may have specified using the "mode" option of :py:func:`convolve`.
2. :meth:`lazylinop.signal.dct` Discrete Cosine Transform of types I to IV
3. :meth:`lazylinop.signal.dst` Discrete Sine Transform of types I to IV
4. :meth:`lazylinop.signal.mdct` Modified DCT
5. :meth:`lazylinop.signal.imdct` Inverse of the Modified DCT
6. :meth:`lazylinop.signal.dft` Discrete Fourier Transform
7. :meth:`lazylinop.signal.wht` Walsh-Hadamard Transform
8. :meth:`lazylinop.signal.dwt` Discrete Wavelet Transform
9. :meth:`lazylinop.signal.idwt` Inverse of the Discrete Wavelet Transform
10. :meth:`lazylinop.signal.dzt` Zak Transform
11. :meth:`lazylinop.signal.stft` Short-Time-Fourier-Transform with a given window
12. :meth:`lazylinop.signal.istft` Inverse of Short-Time-Fourier-Transform
13. :meth:`lazylinop.signal.nufft` Nonuniform Fast-Fourier-Transform
14. :meth:`lazylinop.signal.fnt` Fast-Noiselet-Transform

Padding or cropping
~~~~~~~~~~~~~~~~~~~

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

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

Inverse transforms
~~~~~~~~~~~~~~~~~~

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

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

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

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

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

Other normalizations
~~~~~~~~~~~~~~~~~~~~

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

.. code-block:: python

   >>> 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`` (``norm`` is then either ``'forward'`` or ``'backward'``) requires pre- and/or post-composing by diagonal operators that depend on the type.
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'``):

.. code-block:: python

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


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


Transforms
~~~~~~~~~~

.. autofunction:: lazylinop.signal.convolve
.. autofunction:: lazylinop.signal.dct
.. autofunction:: lazylinop.signal.dst
.. autofunction:: lazylinop.signal.mdct
.. autofunction:: lazylinop.signal.imdct
.. autofunction:: lazylinop.signal.dft
.. autofunction:: lazylinop.signal.wht
.. autofunction:: lazylinop.signal.dwt
.. autofunction:: lazylinop.signal.idwt
.. autofunction:: lazylinop.signal.dzt
.. autofunction:: lazylinop.signal.stft
.. autofunction:: lazylinop.signal.istft
.. autofunction:: lazylinop.signal.nufft
.. autofunction:: lazylinop.signal.fnt


Utility functions
~~~~~~~~~~~~~~~~~

.. autofunction:: lazylinop.signal.chunk
.. autofunction:: lazylinop.signal.decimate
.. autofunction:: lazylinop.signal.downsample
.. autofunction:: lazylinop.signal.overlap_add
.. autofunction:: lazylinop.signal.dwt_coeffs_sizes
.. autofunction:: lazylinop.signal.dwt_to_pywt_coeffs
