# An ADMM algorithm to restore bivariate signals with joint time and covariance regularization

Yusuf Yigit Pilavci, Julien Flamant, Pierre-Antoine Thouvenin, Jérémie Boulanger, Pierre Chainais

## ► To cite this version:

Yusuf Yigit Pilavci, Julien Flamant, Pierre-Antoine Thouvenin, Jérémie Boulanger, Pierre Chainais. An ADMM algorithm to restore bivariate signals with joint time and covariance regularization. XXXe Colloque Francophone de Traitement du Signal et des Images (GRETSI 25), Aug 2025, Strasbourg, France. hal-05262781

HAL Id: hal-05262781

https://hal.science/hal-05262781v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An ADMM algorithm to restore bivariate signals with joint time and covariance regularization

Y. Yiğit PİLAVCI<fnref n="1" /> Julien FLAMANT<fnref n="2" /> Pierre-Antoine THOUVENIN<fnref n="1" /> Jérémie BOULANGER<fnref n="1" /> Pierre CHAINAIS<fnref n="1" />

**Résumé** – Les signaux bivariés occupent une place à part car ils permettent un traitement spécifique tenant compte de leurs propriétés géométriques. Ils correspondent aux ondes bidimensionnelles caractérisées par leur état de polarisation qui traduit la covariance entre les deux composantes du signal. Ce travail propose de faire bénéficier la résolution des problèmes inverses de restauration des signaux bivariés de termes de régularisation portant non seulement sur le signal lui-même, mais aussi sur la covariance entre composantes. Cet objectif amène à travailler avec le signal analytique et à considérer la présence d'un terme quartique dans la fonction de coût. Un problème d'optimisation sous contraintes linéaires est formulé. Une méthode ADMM est proposée impliquant des étapes de mise à jour efficaces. Les simulations numériques illustrent les bonnes performances et la rapidité de la méthode proposée.

**Abstract** – Bivariate signals occupy a special place, as they can be processed specifically to take account of their geometric properties. They correspond to two-dimensional waves characterized by their polarization state, which reflects the covariance between the two signal components. The aim of this work is to enable the resolution of inverse restoration problems for bivariate signals to benefit from regularization terms relating not only to the signal itself, but also to the covariance between components. This objective leads us to work with the analytical signal and to consider the presence of a quartic term in the cost function. An optimization problem under linear constraints is formulated. An ADMM method is proposed, involving efficient updating steps. Numerical simulations illustrate the good performance and speed of the proposed method.

## 1 Introduction

Bivariate signals are involved in many application across various fields, such as astronomy [2, 3], geophysics [14], oceanography [5], and neuroscience [13]. In most of these, the polarization, i.e., the shape of the trajectory drawn by the extremity of the signal vector (e.g., a segment corresponds to a linear polarization, a circle to circular polarization, and in general the instantaneous polarization is elliptical) is a key geometric property. For example, gravitational waves, a prominent instance of bivariate signals, exhibit different polarization properties depending on the properties of the physical sources, e.g., black-holes or large masses [7, 2].

The tools to analyze bivariate signals, and polarization properties in particular, have received growing interest in the last decades. A series of works [8, 12, 15] have theorized and adapted various signal processing tools such as instantaneous frequency, wavelets, and mode decompositions for bivariate signals. Later, the seminal work of [4] studied the use of instantaneous covariance matrices and Stokes parameters. These are alternative signal representations that can be used interchangeably. Each representation provides a different perspective while still conveying the same essential information on the polarization state of the signal. In optics or radar signals, these representations capture the geometrical properties of a wave. Therefore, when dealing with inverse problems, regularization in the polarization domain is necessary. The initial work of [3] formulates the problem for the bivariate signal with known polarization properties. While the solution to this problem is fast to compute, it requires the true polarization of the reconstructed signal to be known in advance, thereby restricting its

applicability. In a recent work [11], this problem is formulated using a regularization promoting a smooth polarization evolution. This approach however entails two main limitations; firstly, the cost function used in the optimization problem is defined in a heuristic manner which is not straightforward to adapt for generic inverse problems. Secondly, the algorithm used to solve the problem is based on gradient descent algorithm, which does not exploit the structure of the problem and whose convergence can be slow in practice.

In this work, we reformulate smoothness priors on the bivariate signals and their polarization. This formulation uses the instantaneous covariance matrices, which leads to a compact formulation of the regularization function. In turn, the overall problem is a quartic optimization problem. To address this issue, we propose an appropriate parameter splitting strategy [1] such that the overall problem can be solved efficiently with an ADMM algorithm. Primal update steps in the algorithm reduce to quadratic convex subproblems, which admit a closed-form solution fast to compute. The organization of the paper is as follows; in Section 2, necessary background and notation related to bivariate signal processing are summarized. Then, we detail the problem definition and the proposed regularization functions. Section 3 introduces the proposed parameter splitting strategy and the ADMM algorithm used to solve the optimization problem considered. Section 4 illustrates the performance of the method with simulations and discusses the results. Finally, concluding remarks are reported in Section 5.

<footnote n="1">¹Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000, Lille, France</footnote>

<footnote n="2">²Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France</footnote>
# 2 Restoration of bivariate signals

## 2.1 Formulation of the problem

Given a discrete, real, narrow-band, bivariate signal $\mathbf{X} = [\mathbf{X}[1]^\top, \dots, \mathbf{X}[N]^\top]^\top \in \mathbb{R}^{N \times 2}$ with $\mathbf{X}[n] = [u[n], v[n]]^\top$, this paper considers the inverse problem of the signal reconstruction by taking both signal and polarization properties into account. A generic forward model for signal reconstruction is:

$$ \mathbf{Y} = \Phi(\mathbf{X}) + \varepsilon \in \mathbb{R}^{N \times D} \quad (1) $$

where $\Phi$ denotes the forward model, $\mathbf{Y}$ is the noisy measurements and $\varepsilon = [\varepsilon_d]_{d \in 1, \dots, D}$ represents the noise with $\varepsilon_d \sim \mathcal{N}(0, \Theta_d)$. For each measurement device $d$, the noise $\varepsilon_d$ is assumed stationary, so that the covariances can be written $\Theta_d = \mathcal{F}^H \Delta_d \mathcal{F}$, with $\mathcal{F}$ the discrete Fourier transform and $\Delta_d$ a diagonal matrix containing the power spectral density (PSD) of the noise for device $d$. For known PSDs, the noise can be whitened in the frequency domain as

$$ \forall d \in \{1, \dots, D\}, \quad \tilde{\mathbf{y}}_d = \Delta_d^{-1/2} \mathcal{F}[\Phi(\mathbf{X})]_d + \tilde{\varepsilon}. \quad (2) $$

The corresponding fidelity term is $f_{\mathbf{Y},\Phi}(\mathbf{X}) = \sum_{d=1}^D \|\Delta_d^{-1/2} \mathcal{F}[\Phi(\mathbf{X})]_d - \mathbf{Y}_d\|_2^2$ Reconstructing $\mathbf{X}$ from $\mathbf{Y}$ can be formulated as the following optimization problem:

$$ \mathbf{X}^* \in \arg \min_{\mathbf{X}} f_{\mathbf{Y},\Phi}(\mathbf{X}) + \sum_{r=1}^{R} \lambda_r g_r(\mathbf{X}), \quad (3) $$

where the data fitting term $f$ includes information about the acquisition system and the distribution of the noise while the second term reflects prior assumptions on the signal $\mathbf{X}$ such as smoothness, sparsity or hard constraints. The following sections details the various proposed regularization penalties.

## 2.2 Regularizations

Each component of the true bivariate signal $\mathbf{X}$ is expected to vary smoothly with time. This prior leads to regularization functions which penalize high variations within each component. In this work, we choose to use $g_1(\mathbf{X}) = \|\mathbf{D}\mathbf{X}\|_F^2 = \sum_{n=2}^N \|\mathbf{X}[n] - \mathbf{X}[n-1]\|_2^2$ where $\mathbf{D}$ is the discrete finite difference operator [6]. Similarly, the polarization of the true signal, thus the corresponding polarizations, smoothly vary. We represent the polarization by the instantaneous covariance matrices which read at time $n$

$$ \Sigma[n] = \mathbf{X}_a[n] \mathbf{X}_a[n]^H \quad (4) $$

where $\mathbf{X}_a$ is the analytic version of $\mathbf{X}$, obtained via the linear transformation

$$ \mathbf{X}_a = \mathbf{A}\mathbf{X}, \quad (5) $$

where $\mathbf{A} = \mathcal{F}^{-1} \Delta_{\omega>0} \mathcal{F}$, $\mathcal{F} \in \mathbb{C}^{N \times N}$ is the discrete Fourier transform and $\Delta_{\omega>0}$ is a diagonal matrix in which each coefficient are 2 for positive frequencies $\omega > 0$, 1 for the DC component, and 0 otherwise [9]. While this matrix encodes the correlations between two components, it has a direct connection with the instantaneous polarization. In other words, its evolution through time encodes the evolution of the polarization properties of the signal [4]. To encode the smoothness

prior, we formulate the second regularization function as:

$$ \begin{aligned} g_2(\mathbf{X}_a) &= \sum_{n=2}^{N} \|\Sigma[n] - \Sigma[n-1]\|_F^2 \\ &= \sum_{n=2}^{N} \|\mathbf{X}_a[n]\mathbf{X}_a[n]^H - \mathbf{X}_a[n-1]\mathbf{X}_a[n-1]^H\|_F^2 \\ &= \sum_{n=2}^{N} \|\mathbf{X}_a\mathbf{J}_n\mathbf{X}_a^H\|_F^2 = \sum_{n=2}^{N} \mathrm{tr}(\mathbf{X}_a^H\mathbf{X}_a\mathbf{J}_n\mathbf{X}_a^H\mathbf{X}_a\mathbf{J}_n), \end{aligned} \quad (6) $$

with $\mathbf{e}_n$ the n-th canonical basis vector and $\mathbf{J}_n = \mathbf{e}_n \mathbf{e}_n^\top - \mathbf{e}_{n-1} \mathbf{e}_{n-1}^\top$. In the rest, we re-denote this regularization function by $\tilde{g}_2(\mathbf{X}_a \mathbf{X}_a^H)$ as it is a function of $\mathbf{X}_a \mathbf{X}_a^H$. This notation will help to explain the proposed method in Section 3.

The regularization $\tilde{g}_2(\mathbf{X}_a \mathbf{X}_a^H)$ is a quartic function of $\mathbf{X}$, which requires the evaluation of $\mathbf{X}_a$ via a dense linear transformation $\mathbf{A}$. This would require large-scale dense linear systems of equations to be solved in the inference algorithm. To avoid this issue, we propose to solve the problem at hand in the domain of analytic signals. This is made possible by adding a hard constraint in the overall optimization problem which appears as an indicator function

$$ g_3(\mathbf{X}) = \iota_A(\mathbf{X}) = \begin{cases} 0 & \text{if } \frac{1}{2}\mathbf{A}\mathbf{X} = \mathbf{X} \\ +\infty & \text{else.} \end{cases} \quad (7) $$

The problem considered is of the form (3), with

$$ \mathbf{X}^* \in \arg \min_{\mathbf{X} \in \mathbb{C}^{N \times 2}} f_{\mathbf{Y},\Phi}(\mathbf{X}) + \lambda_1 g_1(\mathbf{X}) + \lambda_2 \tilde{g}_2(\mathbf{X}\mathbf{X}^H) + \iota_A(\mathbf{X}). \quad (8) $$

This is a constrained non-convex quartic problem, which does not admit a closed-form solution. We assume that the data fidelity term $f_{\mathbf{Y},\Phi}(\mathbf{X})$ is equivalent to $f_{\mathbf{A}\mathbf{Y},\Phi}(\mathbf{A}\mathbf{X})$ when reformulated in the analytic signal domain. This is verified if $\mathbf{A}\Phi(\mathbf{X}) = \Phi(\mathbf{A}\mathbf{X})$ which is often the case in applications of bivariate signal processing. For example, in Section 4, we use the forward model of gravitational waves, $\Phi(\mathbf{X}) = \mathbf{X}\mathbf{F}$, which verifies this assumption.

# 3 Variable splitting and ADMM

The problem (8) is prone to an ADMM approach (See Alg. 1). The resulting variable splitting strategy is twofold. The first applies to the quartic regularization term to simplify it to a quadratic one. In particular, we replace $\tilde{g}(\mathbf{X}\mathbf{X}^H)$ with $\tilde{g}_2(\mathbf{Z}_2\mathbf{X}^H)$ which is quadratic w.r.t. $\mathbf{X}$ and the splitting variable $\mathbf{Z}_2$. The second divides the objective functions into certain blocks so that the corresponding subproblems are cheap to solve. Specifically, two more auxilary variables $\mathbf{Z}_1$ and $\mathbf{Z}_3$ respectively replaces $\mathbf{X}$ as $f(\mathbf{Z}_1)$, $g_1(\mathbf{Z}_1)$ and $\iota_A(\mathbf{Z}_3)$. This is heavily inspired by a classical optimization scheme called (regularized) consensus optimization [1] in which the global variable $\mathbf{X}$ is replaced by the auxilary variables $\mathbf{Z}_i$'s for each block of the objective function. After splitting, (8) becomes

$$ \mathbf{X}^* \in \arg \min_{\substack{\mathbf{X} \in \mathbb{C}^{N \times 2} \\ \mathbf{Z}_i = \mathbf{X}, \\ i \in \{1, 2, 3\}}} \frac{f_{\mathbf{Y},\Phi}(\mathbf{Z}_1) + \lambda_1 g_1(\mathbf{Z}_1) + \lambda_2 \tilde{g}_2(\mathbf{Z}_2 \mathbf{X}^H) + \iota_A(\mathbf{Z}_3)}{\mathcal{L}(\mathbf{X}, \mathbf{Z}_1, \mathbf{Z}_2, \mathbf{Z}_3)}. \quad (9) $$
<caption>**Algorithm 1: ADMM for Eq. (9)**</caption>

**Input:** Y
**Initialize randomly:** $X^{(0)}, \{Z_i^{(0)}, U_i^{(0)}\}_{1 \le i \le 3}$
**for** $k = 1, 2, \dots$ **until** convergence **do**
> **Primal parameter update**
> &nbsp;&nbsp;&nbsp;&nbsp; $X^{(k+1)} = \arg \min_{\mathbf{X}} \mathcal{L}(\mathbf{X}, Z_{1,2,3}^{(k)}) + \frac{\rho}{2} \sum_{i=1}^{3} \|Z_i^{(k)} - \mathbf{X} + U_i^{(k)}\|_F^2$
> **Dual parameter updates**
> &nbsp;&nbsp;&nbsp;&nbsp; $Z_1^{(k+1)} = \arg \min_{Z_1} \mathcal{L}(X^{(k+1)}, Z_1, Z_{2,3}^{(k)}) + \frac{\rho}{2} \|Z_1 - X^{(k+1)} + U_1^{(k)}\|_F^2$
> &nbsp;&nbsp;&nbsp;&nbsp; $Z_2^{(k+1)} = \arg \min_{Z_2} \mathcal{L}(X^{(k+1)}, Z_1^{(k+1)}, Z_2, Z_3^{(k)}) + \frac{\rho}{2} \|Z_2 - X^{(k+1)} + U_2^{(k)}\|_F^2$
> &nbsp;&nbsp;&nbsp;&nbsp; $Z_3^{(k+1)} = \arg \min_{Z_3} \mathcal{L}(X^{(k+1)}, Z_{1,2}^{(k+1)}, Z_3) + \frac{\rho}{2} \|Z_3 - X^{(k+1)} + U_3^{(k)}\|_F^2$
> **Multiplier updates**
> &nbsp;&nbsp;&nbsp;&nbsp; $U_i^{(k+1)} = U_i^{(k)} + Z_i^{(k+1)} - X^{(k+1)}$, for $i \in \{1,2,3\}$
**end**
**return** $X^{(k)}$

This divide-to-conquer strategy permits to split the overall problem into subproblems that admit closed-form solutions. For example, the problem composed by $f$ and $g_1$ often leads to signal filtering which can be implemented efficiently either in time or frequency domain. Similarly, separating quadratic input parameters in $\tilde{g}_2$ leads to solving quadratic problems w.r.t. $X$ or $Z_2$ instead of a quartic one. Finally, solving for $Z_3$ boils down a projection onto the space of analytic signals which can be efficiently implemented. The update step for $X$ at iteration $k$ requires to solve a sparse linear system:

$$ \underbrace{\left( \lambda_2 \sum_{i=2}^{N} \mathbf{J}_t \mathbf{Z}_2^{(k)} (\mathbf{Z}_2^{(k)})^{\mathrm{H}} \mathbf{J}_t^{\top} + 3\rho \mathbf{I} \right)}_{\mathbf{M}_x} \mathbf{X} = \underbrace{\sum_{i=1}^{3} \rho (\mathbf{Z}_i^{(k)} - \mathbf{U}_i^{l})}_{\mathbf{b}_x} \quad (10) $$

where $\mathbf{M}_x$ is a triangular matrix. Similarly, the $Z_2$-update is obtained by solving $\mathbf{M}_{z_2} \mathbf{Z}_2 = \mathbf{b}_{z_2}$ with $\mathbf{M}_{z_2} = (\lambda_2 \sum_{i=2}^{N} \mathbf{J}_t \mathbf{X}^{(k+1)} (\mathbf{X}^{(k+1)})^{\mathrm{H}} \mathbf{J}_t^{\top} + \rho \mathbf{I})$ and $\mathbf{b}_{z_2} = \rho(\mathbf{X}^{(k+1)} + \mathbf{U}_2^l)$. The solution to $Z_3$-update is independent from $\rho$ and a projection onto the analytic signal subspace, i.e., $Z_3 = \mathrm{proj}_A(\mathbf{X} + \mathbf{U}_3^{(k)}) = \frac{1}{2}\mathbf{A}(\mathbf{X} + \mathbf{U}_3^{(k)})$. Finally, the updates on $Z_1^{(k+1)}$ take the form of the Tikhonov regularization, which can be often interpreted as a spectral filtering. This will be exemplified later in Section 4 for simulation results.

## 4 Experimental results

For the illustrative results, we consider the following linear forward model of the gravitational waves [3]:

$$ \mathbf{Y} = \mathbf{X}\mathbf{F} + \varepsilon \quad (11) $$

where $\mathbf{F} = [\mathbf{f}_1 | \dots | \mathbf{f}_D] \in \mathbb{R}^{2 \times D}$ is the linear projection model that depicts the behaviours of different antennas and the noise in each channel/antenna $d$, denoted by $\varepsilon_d$, is assumed to be colored but known up to its power spectral density $\Delta_d$.

In the experiments, the bivariate signals are randomly generated by using the model in [4, (1.57)] using smooth polarization parameters. The generated signals are narrow-band, AM-FM-PM signals with the main carrying frequency linearly increasing in [25, 100] Hz. The forward model matrix $\mathbf{F}$ is

<figure><img src="image_2.png" /><figcaption>Figure 1 – The proposed method with different parameter configurations: MLE solution i.e., $\lambda_{1,2} = 0$ (top), covariance i.e., $\lambda_1 = 0$ (second row) and time smoothing i.e., $\lambda_2 = 0$ (third row), the proposed method i.e., $\lambda_{1,2} > 0$ (bottom)</figcaption></figure>

<table><caption>Table 1 – Runtime (sec) for varying sample size $N$.</caption><thead><tr><th>N</th><th>MLE</th><th>STS</th><th>Prop. method</th></tr></thead><tbody><tr><td>512</td><td>0.0003</td><td>3.48</td><td>1.78</td></tr><tr><td>1024</td><td>0.0004</td><td>3.90</td><td>2.38</td></tr><tr><td>4096</td><td>0.0013</td><td>6.20</td><td>6.91</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2 – Reconstruction SNR of the different methods.</figcaption></figure>
entry-wise sampled from i.i.d normal distributions. To gen-
erate noise $\varepsilon_d$, we randomly sample its amplitude $\alpha_d[f] \sim U(0, \sigma)$ and its phase $\varphi_d[f] \sim U(0, 2\pi)$ for each frequency $f$. The noise is then obtained as $\varepsilon_d = F^{-1}\{\alpha_d \exp(j\varphi_d)\}$. The corresponding spectral whitening operator reads $T_d = \text{diag}(\alpha_d)^{-1}$ where diag generates a diagonal matrix whose diagonal are given by the input vector. The range of pa- rameter $\sigma$ is set to $\{10^{-1}, \dots, 10^1\}$ and the simulations are repeated for several realizations of $N = 512, 1024, 4096$. The hyperparameters $\lambda_1 \in \{10^{-1}, 1, 10, 10^2, 10^3\}$ and $\lambda_2 \in \{10^2, 10^3, 10^4, 10^5, 10^6\}$ are set by a grid-search. The pro- posed method is compared with the denoising approach pro- posed in [11] refered as Stokes and signal smoothing (STS). The algorithm in [11] requires a noisy version of the bivariate signal as input thus we pass the maximum likelihood solution for the forward model in (11) ($\lambda_{1,2} = 0$). Finally, we com- pare these approaches in terms of runtime and reconstruction SNR i.e., r-SNR(·) = 10 log $\frac{\|\mathbf{X}\|_F^2}{\|\mathbf{X}_-\|_F^2}$.

Fig. 1 illustrates the restoration performance of the proposed method for different parameter configurations, namely covariance smoothing ($\lambda_1 = 0$), time smoothing ($\lambda_2 = 0$), MLE ($\lambda_{1,2} = 0$) and the proposed method with all regularizations activated ($\lambda_{1,2} > 0$). The best performances are achieved by regularizing both in signal and polarization domains (r-SNR=8.53 dB) whereas smoothing only in time domain (r-SNR=5.59 dB) or only in polarization (r-SNR=7.07 dB) is not sufficient. In addition, Fig. 2 shows that the covariance smoothing is more relevant in high SNR regimes whereas in low SNR regime, time smoothing is more relevant. Moreover, the proposed method usually performs better than the STS method, while taking comparable or less amount of time to reach to solution (see Table 1).

5 Conclusion

A new formulation is proposed to take the polarization prop-
erties into account for solving inverse problems involving
bivariate signals. This formulation admits an efficient ADMM
algorithm where each step has a closed-form solution and can
be calculated with fast algorithms. The existing method [11]
relies on an algorithm based on gradient descent, while the
proposed approach exploits the quartic structure of the prob-
lem. As a result, the proposed method yields superior signal
restoration results, while scaling comparably with existing
algorithms. In the future, we plan to adapt suitable hyperpa-
rameter tuning methods [10] for selecting $\lambda_1$ and $\lambda_2$ and apply
this approach for real-life problems such as the reconstruction
of gravitational waves [2].

References

[1] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers. *Foundations and Trends*® in *Machine learning*, 3(1):1–122, 2011.

[2] Cyril Cano. *Mathematical tools and signal processing algorithms for the study of gravitational waves polariza-

tion. PhD thesis, Université Grenoble Alpes [2020-....],
2022.

[3] Cyril Cano, Nicolas Le Bihan, and Eric Chassande-Mottin. Contrainte de polarisation pour la reconstruction d'ondes gravitationnelles. In *GRETSI 2023-XXIXème Colloque Francophone de Traitement du Signal et des Images*, 2023.

[4] Julien Flamant. Une approche générique pour l'analyse et le filtrage des signaux bivariés. PhD thesis, Ecole Centrale de Lille, 2018.

[5] Joseph Gonella. A rotary-component method for analysing meteorological and oceanographic vector time series. In *Deep sea research and oceanographic abstracts*, volume 19, pages 833–846. Elsevier, 1972.

[6] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction, 2017.

[7] Maximiliano Isi. Parametrizing gravitational-wave polarizations. *Classical and Quantum Gravity*, 40(20):203001, 2023.

[8] Jonathan M Lilly and Sofia C Olhede. Bivariate instantaneous frequency and bandwidth. *IEEE Transactions on Signal Processing*, 58(2):591–603, 2009.

[9] Lawrence Marple. Computing the discrete-time" analytic" signal via FFT. *IEEE Transactions on signal processing*, 47(9):2600–2603, 1999.

[10] Marcelo Pereyra, José M Bioucas-Dias, and Mário AT Figueiredo. Maximum-a-posteriori estimation with unknown regularisation parameters. In *2015 23rd European Signal Processing Conference (EUSIPCO)*, pages 230–234. IEEE, 2015.

[11] Yusuf Yiğit Pılavci, Jérémie Boulanger, Pierre-Antoine Thouvenin, and Pierre Chainais. Denoising bivariate signals via smoothing and polarization priors. In *2024 32nd European Signal Processing Conference (EUSIPCO)*, pages 2602–2606. IEEE, 2024.

[12] Gabriel Rilling, Patrick Flandrin, Paulo Gonçalves, and Jonathan M Lilly. Bivariate empirical mode decomposition. *IEEE signal processing letters*, 14(12):936–939, 2007.

[13] Vangelis Sakkalis. Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG. *Computers in biology and medicine*, 41(12):1110–1117, 2011.

[14] JC Samson. Pure states, polarized waves, and principal components in the spectra of multiple, geophysical time-series. *Geophysical Journal International*, 72(3):647–664, 1983.

[15] Adam M Sykulski, Sofia Charlotta Olhede, Jonathan M Lilly, and Jeffrey J Early. Frequency-domain stochastic modeling of stationary bivariate or complex-valued signals. *IEEE Transactions on Signal Processing*, 65(12):3136–3151, 2017.