ENH: Determine arbitrary dual windows for short-time Fourier Transform

Provide a new function `closest_STFT_dual_window()` to calculate the closest dual window in a least-squares sense for a given desired window as well as a new method `ShortTimeFFT.from_win_equals_dual()` to determine windows which are equal to their dual window.

doc/source/tutorial/signal.rst

@@ -1122,6 +1122,7 @@ inverse (as implemented using efficient FFT calculations in the :mod:`scipy.fft`
 is given by
 
 .. math::
+    :label: eq_SpectA_FFT
 
     X_l := \sum_{k=0}^{n-1} x_k \e^{-2\jj\pi k l / n}\ ,\qquad
     x_k = \frac{1}{n} \sum_{l=0}^{n-1} X_l \e^{2\jj\pi k l / n}\ .
@@ -1645,7 +1646,7 @@ Eq. :math:numref:`eq_STFT_DFT` can be expressed as
 .. math::
 
     \vb{s}_p = \vb{F}\,\vb{x}_p \quad\text{with}\quad
-    F[q,m] =\exp\!\big\{-2\jj\pi (q + \phi_m)\, m / M\big\}\ ,
+    F[q,m] = \frac{1}{\sqrt{M}}\exp\!\big\{-2\jj\pi (q + \phi_m)\, m / M\big\}\ ,
 
 which allows the STFT of the :math:`p`-th slice to be written as
 
@@ -1655,8 +1656,10 @@ which allows the STFT of the :math:`p`-th slice to be written as
     \vb{s}_p = \vb{F}\vb{W}_{\!p}\,\vb{x} =: \vb{G}_p\,\vb{x}
     \quad\text{with}\quad s_p[q] = S[p,q]\ .
 
-Note that :math:`\vb{F}` is unitary, i.e., the inverse equals its conjugate
-transpose meaning :math:`\conjT{\vb{F}}\vb{F} = \vb{I}`.
+Due to the scaling factor of :math:`M^{-1/2}`, :math:`\vb{F}` is unitary, i.e., the
+inverse equals its conjugate transpose meaning :math:`\conjT{\vb{F}}\vb{F} = \vb{I}`.
+Other scalings, e.g., like in Eq. :math:numref:`eq_SpectA_FFT`, are allowed as well,
+but would in this section make the notation slightly more complicated.
 
 To obtain a single vector-matrix equation for the STFT, the slices are stacked
 into one vector, i.e.,
@@ -1683,6 +1686,7 @@ equation the Moore-Penrose inverse :math:`\vb{G}^\dagger` can be utilized
 which exists if
 
 .. math::
+   :label: eq_STFT_MoorePenrose_DD
 
     \vb{D} := \conjT{\vb{G}}\vb{G} =
         \begin{bmatrix}
@@ -1758,26 +1762,131 @@ zeros, enlarging :math:`\vb{U}` so all slices which touch :math:`x[k]` contain
 the identical dual window
 
 .. math::
+   :label: eq_STFT_CanonDualWin
 
-    w_d[m] = w[m] \inv{\sum_{\eta\in\IZ} \big|w[m + \eta\, h]\big|^2}\ .
+    w_d[m] = w[m] \inv{\sum_{\eta\in\IZ} \big|w[m + \eta\, h]\big|^2}
+           =\ w[m] \inv{\sum_{l=0}^{M-1} \big|w[l]\big|^2 \delta_{l+m,\eta h}}\ ,
+             \quad \eta\in \IZ\ .
 
 Since :math:`w[m] = 0` holds for :math:`m \not\in\{0, \ldots, M-1\}`, it is
 only required to sum over the indexes :math:`\eta` fulfilling
-:math:`|\eta| < M/h`. The name dual window can be justified by inserting Eq.
-:math:numref:`eq_STFT_Slice_p` into Eq. :math:numref:`eq_STFT_istftM`, i.e.,
+:math:`|\eta| < M/h`. The second expression is an alternate form by summing from index
+:math:`l=0` to :math:`M` in index increments of :math:`h`.
+:math:`\delta_{l+m,\eta h}` is Kronecker delta notation for ``(l+m) % h == 0``. The
+name dual window can be justified by inserting Eq. :math:numref:`eq_STFT_Slice_p` into
+Eq. :math:numref:`eq_STFT_istftM`, i.e.,
 
 .. math::
+   :label: eq_STFT_WindDualCond0
 
     \vb{x} = \sum_{p=0}^{P-1} \conjT{\vb{U}_p}\,\conjT{\vb{F}}\,
                                                  \vb{F}\,\vb{W}_{\!p}\,\vb{x}
         = \left(\sum_{p=0}^{P-1} \conjT{\vb{U}_p}\,\vb{W}_{\!p}\right)\vb{x}\ ,
 
-showing that :math:`\vb{U}_p` and :math:`\vb{W}_{\!p}` are interchangeable.
-Hence, :math:`w_d[m]` is also a valid window with dual window :math:`w[m]`.
-Note that :math:`w_d[m]` is not a unique dual window, due to :math:`\vb{s}`
-typically having more entries than :math:`\vb{x}`. It can be shown, that
-:math:`w_d[m]` has the minimal energy (or :math:`L_2` norm) [#Groechenig2001]_,
-which is the reason for being named the  "canonical dual window".
+showing that :math:`\vb{U}_p` and :math:`\vb{W}_{\!p}` are interchangeable. Due
+:math:`\vb{U}_p` and :math:`\vb{W}_{\!p}` having the same structure given in Eq.
+:math:numref:`eq_STFT_WinMatrix`, Eq. :math:numref:`eq_STFT_WindDualCond0` contains a
+general condition for all possible dual windows, i.e.,
+
+.. math::
+   :label: eq_STFT_AllDualWinsCond0
+
+   \sum_{p=0}^{P-1} \conjT{\vb{W}_{\!p}}\,\vb{U}_p = \vb{I}
+   \quad\Leftrightarrow\quad
+   \sum_{p=0}^{P-1} \sum_{m=0}^{M-1}
+   \conj{w[m]} u[m]\, \delta_{m+ph,i}\, \delta_{m+ph,j} = \delta_{i,j}
+
+which can be reformulated into
+
+.. math::
+   :label: eq_STFT_AllDualWinsCond
+
+   \conjT{\vb{V}}\,\vb{u} = \vb{1}\ , \qquad
+   V[i,j] :=  w[i]\, \delta_{i, j+ph}\ , \qquad
+   \vb{V}\in \IC^{M\times h}\ ,
+
+where :math:`\vb{1}\in\IR^h` is a vector of ones. The reason
+that :math:`\vb{V}` has only :math:`h` columns is that the :math:`i`-th and
+:math:`(i+m)`-th row, :math:`i\in\IN`, in Eq. :math:numref:`eq_STFT_WindDualCond0` are
+identical. Hence there are only :math:`h` distinct equations and one non-zero entry per
+row in :math:`\vb{V}`.
+
+Of practical interest is finding the valid dual window :math:`\vb{u}_d` closest to a
+given vector :math:`\vb{d}\in\IC^M`. By utilizing an :math:`h`-dimensional vector
+:math:`\vb{\lambda}` of Lagrange multipliers, we obtain the convex quadratic
+programming problem
+
+.. math::
+
+    \vb{u}_d = \min_{\vb{u}}{ \frac{1}{2}\lVert\vb{d} - \vb{u}\rVert^2 }
+    \quad\text{w.r.t.}\quad  \conjT{\vb{V}}\vb{u} = \vb{1}
+    \qquad\Leftrightarrow\qquad
+     \begin{bmatrix} \vb{I} & \vb{V}\\
+                     \conjT{\vb{V}} & \vb{0}        \end{bmatrix}
+    \begin{bmatrix} \vb{u} \\ \vb{\lambda} \end{bmatrix} =
+    \begin{bmatrix} \vb{d}\\ \vb{1} \end{bmatrix}\ .
+
+A closed form solution can be calculated by inverting the :math:`2\times2` block matrix
+symbolically, which gives
+
+.. math::
+   :label: eq_STFT_ClosestDual
+
+   \vb{u}_d &= \underbrace{\vb{V}\inv{\conjT{\vb{V}}\vb{V}}}_{%
+                                                   =:\vb{Q}\in\IC^{M\times h}} \vb{1} +
+                \big(\vb{I} - \vb{Q}\conjT{\vb{V}}\big)\vb{d}
+            = \vb{w}_d + \vb{d} - \vb{q}_d\ ,\qquad
+   \vb{w}_d = \vb{Q}\vb{1}\ ,\qquad
+   \vb{q}_d = \vb{Q}\conjT{\vb{V}}\vb{d}\ ,\\
+   Q[i,j] &= \frac{ w[i] \delta_{i-j, \eta h} }{%
+                             \sum_{k=0}^M |w[k]|^2\,\delta_{i-k, \xi h} } ,\qquad
+   w_d[i] = \frac{w[i] }{ \sum_{l=0}^M \conj{w[l]} w[l] \delta_{i-l, \xi h} }\ ,\\
+   q_d[i] &= w[i] \frac{\sum_{j=1}^h \delta_{i-j, \eta h}
+                        \sum_{l=1}^M \conj{w[l]} d[l] \delta_{j-l, \xi h} }{%
+                        \sum_{l=0}^M \conj{w[l]} w[l] \delta_{i-l, \zeta h} }
+           = w_d[i] \sum_{l=1}^M \conj{w[l]} d[l] \delta_{i-l, \eta h}\ ,
+
+with :math:`\eta,\xi,\zeta\in\IZ`. Note that the first term :math:`\vb{w}_d` is equal
+to the solution given in Eq. :math:numref:`eq_STFT_CanonDualWin` and that the inverse
+of :math:`\conjT{\vb{V}}\vb{V}` must exist or else the STFT is not invertible. When
+:math:`\vb{d}=\vb{0}`, the solution :math:`\vb{w}_d` is obtained. Hence,
+:math:`\vb{w}_d` minimizes the :math:`L_2`-norm :math:`\lVert\vb{u}\rVert`, which is
+the justification for its name "canonical dual window". Sometimes it is more desirable
+to find the closest vector in regard to direction and ignoring the vector length. This
+can be achieved by introducing a scaling factor :math:`\alpha\in\IC` to minimize
+:math:`\lVert\alpha\vb{d} - \vb{u}\rVert^2`. Since Eq.
+:math:numref:`eq_STFT_ClosestDual` already provides a general solution, we can write
+
+.. math::
+
+    \alpha_{\min} = \min_{\alpha}{
+         \frac{1}{2}\big\lVert\alpha\vb{d} - \vb{u}(\alpha)\big\rVert^2 }\ ,\qquad
+   \vb{u}(\alpha) = \vb{w}_d + \alpha\vb{d} - \alpha\vb{q}_d
+   \qquad\Leftrightarrow\qquad
+   \alpha_{\min} = \conjT{\vb{q}_d}\vb{w}_d \big/\,
+                    \conjT{\vb{q}_d}\vb{q}_d \ .
+
+The case where the window :math:`w[m]` and the dual window :math:`u[m]` are equal can
+be easily derived from Eq. :math:numref:`eq_STFT_AllDualWinsCond0` resulting in
+:math:`h` conditions of the form
+
+.. math::
+   :label: eq_STFT_EqualWindDualCond
+
+
+   \sum_{p=0}^{\lfloor M / h \rfloor} \big|w[m+ph]\big|^2 = 1\ ,
+                                                     \qquad m \in \{0, \ldots, h-1\}\ .
+
+Note that each window sample :math:`w[m]` appears only once in the :math:`h` equations.
+To find a closest window :math:`\vb{w}` for given window :math:`\vb{d}` is
+straightforward: Partition :math:`\vb{d}` according to Eq.
+:math:numref:`eq_STFT_EqualWindDualCond` and normalize the length of each partition to
+unity. Note that if Eq. :math:numref:`eq_STFT_EqualWindDualCond` holds, the matrix
+:math:`\vb{D}` of Eq. :math:numref:`eq_STFT_MoorePenrose_DD` is the identity matrix
+making the STFT :math:`\vb{G}` a unitary mapping, i.e.,
+:math:`\conjT{\vb{G}}\vb{G}=\vb{I}`. In this case :math:`w[m]` is also a canonical dual
+window, which can be seen by investigating Eq. :math:numref:`eq_STFT_CanonDualWin`.
+
 
 
 .. _tutorial_stft_legacy_stft:

scipy/signal/__init__.py

@@ -280,9 +280,11 @@
    ShortTimeFFT   -- Interface for calculating the \
                      :ref:`Short Time Fourier Transform <tutorial_stft>` and \
                      its inverse.
+   closest_STFT_dual_window -- Calculate the STFT dual window of a given window \
+                               closest to a desired dual window.
    stft           -- Compute the Short Time Fourier Transform (legacy).
    istft          -- Compute the Inverse Short Time Fourier Transform (legacy).
-   check_COLA     -- Check the COLA constraint for iSTFT reconstruction.
+   check_COLA     -- Check the COLA constraint for iSTFT reconstruction (legacy).
    check_NOLA     -- Check the NOLA constraint for iSTFT reconstruction.
 
 Chirp Z-transform and Zoom FFT