text:
\begin{document}

\title{Colored Gaussian Multiple Descriptions: Spectral-Domain Characterization and Time-Domain Design}

\author{Jan {\O}stergaard,~\IEEEmembership{Member,~IEEE,}
 Yuval Kochman,~\IEEEmembership{Member,~IEEE,}
         and~Ram~Zamir,~\IEEEmembership{Fellow,~IEEE}
 \thanks{This work was presented in part at the IEEE Data Compression Conference, Snowbird, Utah, 2008.}
\thanks{J.~{\O}stergaard (janoe@ieee.org) is with the Department of Electronic Systems, Aalborg University, Aalborg, Denmark.
 The work of~J.~{\O}stergaard is supported by the Danish Research Council for Technology and Production Sciences, grant no.\ 274-07-0383.}
 \thanks{Y.~Kochman (yuvalko@mit.edu) is with the Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.}
 \thanks{R.\ Zamir (zamir@eng.tau.ac.il) is with the Department of Electrical Engineering-Systems, Tel Aviv University, Tel Aviv, Israel.}
 }

\markboth{Submitted to IEEE Transactions on Information Theory, May 2010}{}

\maketitle

\begin{abstract}

It is well known that Shannon's rate-distortion function (RDF) in the colored quadratic Gaussian (QG) case, can be parametrized via a single Lagrangian variable (the ``water level'' in the reverse water filling solution).
In this work, we show that the symmetric colored QG multiple-description (MD) RDF in the case of two descriptions, can be parametrized via two Lagrangian variables.
To establish this result, we use two key ideas.
First, we propose a new representation for the MD test channel, and show that the mutual information rate across this channel coincides with the QG MD RDF. Second, we use variational calculus to obtain a spectral domain representation of the test channel's optimal side and central distortion spectra given the source spectral density and the side and central distortion constraints.
The distortion spectra are specified via two Lagrangian parameters, which control the trade-off between the side distortion, the central distortion, and the coding rate.

We also show that the symmetric colored QG MD RDF can be achieved by
combining noise-shaped predictive coding, dithered quantization, and memoryless entropy coding.
In particular, we show that the proposed MD test channel can
be materialized by embedding two source prediction loops, one for
each description, within a common noise shaping loop whose parameters are explicitly found from the spectral-domain characterization.
The source prediction loops exploit the source memory, and thus reduce the
coding rate. The noise-shaping loop controls the trade-off between the side and the
central distortions, by shaping the quantization noise.

\end{abstract}

\begin{IEEEkeywords}
Multiple-description coding, rate-distortion theory, predictive coding, noise shaping, delta-sigma quantization,
optimization, calculus of variations.
\end{IEEEkeywords}

\IEEEpeerreviewmaketitle

\section{Introduction}
\IEEEPARstart{T}{he}
traditional multiple-description (MD) problem
\cite{elgamal:1982,ozarow:1980} considers a source, which is encoded
into two descriptions that are transmitted over separate
channels. Either one of the channels may break down and thereby cause
a description loss at the receiver.
The decoder knows which channels that are working but the encoder does
not.
The problem is then to design the two descriptions so that they individually represent the source to within some desired distortion level and yet are able to refine each other. Thus, combining the descriptions improves upon their individual performances.
In the \emph{symmetric} situation, the two descriptions are balanced, i.e., they are encoded at the same coding rate, and leads to the same distortion when used separately at the decoder.

The achievable rate-distortion region for the MD problem is only
completely known for very few cases.
El-Gamal and Cover~\cite{elgamal:1982} presented an achievable rate region for two
descriptions and memoryless sources and Ozarow~\cite{ozarow:1980}
showed that in the white quadratic Gaussian (QG) case, i.e., for white
Gaussian sources and MSE distortion, the El-Gamal and Cover
region is tight.
An achievable rate region for the case of general (time-correlated) stationary Gaussian sources and mean-squared error (MSE) distortion was recently characterized by
Chen et al.~\cite{chen:2009}.
In particular, it was shown in~\cite{chen:2009} that the achievable rate region forms a closed and convex set and that the minimal description rates can be found by extremizing over all distortion spectra satisfying the individual side and central distortion constraints. No explicit solution to the optimal distortion spectra was found. However, some intuition towards a spectral domain characterization was provided. Specifically,
it was shown that the optimal rates for stationary Gaussian sources
can be expressed as the sum of rates of parallel channels, each one
representing a frequency band. Each of the channels must be tuned to a
minimum Ozarow MD rate for some frequency dependent distortion
level. In some sense, this can be seen as a reverse ``water-filling''
approach where instead of having a flat water level as in the
conventional single-description (SD) case, the water level is frequency dependent. The authors also pointed out that obtaining an explicit spectral domain solution from their results is technically non-trivial.
Instead it was argued that the optimal rates can be found through numerical optimization by approximating the source spectral density by piece-wise constant functions. However, in general, for arbitrarily shaped sources, this becomes an infinite-dimensional optimization problem.

In this paper, we present a parametrization of the symmetric colored QG MD rate-distortion function (RDF). While Shannon's RDF in the SD case can be parametrized by a single Lagrangian variable~\cite{kolmogorov:1956} (usually referred to as a ``water level''), we show here that the symmetric colored
QG MD RDF can be parametrized via two Lagrangian variables.\footnote{In our case, however, the two parameters cannot generally be interpreted as ``water levels''.}
To establish this result, we use two key ideas.
We propose a new representation for the MD test channel (see e.g., Fig.~\ref{DPCM_DSM_2}), and show that the mutual information rate across this channel coincides with the QG MD RDF. Moreover, the mutual information rate is shown to be equal to the scalar mutual information over an additive white Gaussian noise (AWGN) channel, and the test channel can therefore be realized with white Gaussian quantization (e.g., high dimensional lattice quantization).

This still leaves open the optimization of the spectral parameters,
which minimize the mutual information rate.

Instead of taking the conventional approach of diagonalizing the colored Gaussian source  and thereby obtain an infinite number of independent sources (which might result in an infinite-dimensional optimization problem) we show that it is possible (and feasible) to directly optimize over the continuum of the test channels' side and central distortion spectral densities through the use of calculus of variations~\cite{weinstock:1974}.
The resulting distortion spectra are then specified via two Lagrangian parameters, which control the trade-off between the side distortion, the central distortion, and the coding rate.
Thus, we avoid extremizing over all \emph{functions} representing admissible distortion spectra subject to the two distortion constraints.
Instead, our results reveal that this otherwise intractable infinite-dimensional optimization problem, can be cast as a two-dimensional optimization problem over
two non-negative Lagrangian parameters subject to the same distortion constraints.

In~\cite{ostergaard:2009}, it was shown that Ozarow's white Gaussian
MD RDF can be achieved by noise-shaped coding based on dithered Delta-Sigma
quantization (DSQ) and then followed by memoryless entropy coding. Furthermore, by exploiting the fact
that Ozarow's test channel becomes asymptotically optimal for
stationary sources in the high-rate regime~\cite{zamir:2000}, it was
shown in \cite{ostergaard:2009} that, at high resolution, the
stationary MD RDF is achievable by DSQ and \emph{joint} entropy
coding.
In~\cite{chen:2009} it is demonstrated how one can achieve any
point on the boundary of the colored Gaussian achievable rate
region by a \emph{frequency-domain} scheme, where the source is divided
into sub-bands, and in each sub-band the ``quantization-splitting'' scheme for white Gaussian sources presented in~\cite{chen:2006}
is applied.

In this paper, we propose a \emph{time-domain} approach: We show that the symmetric colored QG MD RDF can be achieved by noise-shaped predictive coding and \emph{memoryless} dithered quantization (in the limit of high dimensional quantization)
at all resolutions and all side-to-central distortion ratios.
We establish this result by forming a nested prediction / noise-shaping structure
containing a dithered DSQ scheme similar to~\cite{ostergaard:2009}
in the outer loop and a predictive coder per each description in
the inner loop, see for example Fig.~\ref{high_res}.
Each of the
predictive coders has the structure of the differential pulse-code modulation (DPCM) scheme, shown to be
optimal in the SD setting in~\cite{zamir:2008}.
\footnote{The idea of exploiting prediction in MD coding has previously been
proposed by other authors, see for example the following related
works~\cite{ingle:1995,regunthan:2000,vaishampayan:1999,nathan:2001}.
All these works faced the basic problem: Since DPCM uses prediction
from the reconstruction rather than from the source itself, and this
prediction should be reproduced at the decoder, it is not clear
which of the possible reconstructions should be used for prediction.
The present work solves this problem.}
The role of the DSQ loop is to shape the quantization noise so that
a desired trade-off between the side distortions and the central
distortion is achieved. It was shown in \cite{ostergaard:2009} that
the central distortion is given by the power of the noise that falls
within the in-band spectrum  (i.e.\ the part of the frequency
spectrum which overlaps the source spectrum) whereas the side
distortion is given by the power of the complete noise spectrum,
i.e.\ the in-band and the out-of-band noise spectrum.
It was furthermore shown that any ratio of side-to-central distortion can
be obtained by proper shaping of the quantization noise.
We establish a similar result here.
In particular, the predictive coders exploit the source memory
and thereby minimize the coding rate and make sure that memoryless
entropy coding is optimal. Moreover, the DSQ loop performs the noise
shaping, which is required in order to achieve any desired pair of
side and central distortions. At general resolutions, the optimal
noise shaping depends upon the source spectrum. However, at high
resolutions, the optimal noise shaping becomes independent of the source
spectrum and converges to a piece-wise constant function with a single jump
discontinuity.

This paper is organized as follows.
In section~\ref{sec:prelim}, we provide the preliminaries.
Then, in Section~\ref{sec:qrdf}, we propose a test channel, which provides a
new interpretation of the QG MD RDF.
We derive the optimal distortion spectra in Section~\ref{sec:spectral}.
With the test channel in mind, we present, in Section~\ref{sec:mask}, an SD time-domain scheme which encodes a source subject to a distortion mask.
Then, in Section~\ref{sec:mdscheme}, we extend the SD time-domain scheme of Section~\ref{sec:mask} to the MD case.
Conclusions are in Section~\ref{sec:conclusion}.
Longer proofs are deferred to the appendices.

\section{Preliminaries}\label{sec:prelim}

\subsection{Stationary Gaussian Processes and Spectral Decomposition}

Let
$X=\{X[n]\}_{n=0}^{\infty}$ be a zero-mean discrete-time stationary autoregressive Gaussian process with power spectral density
$S_X=\{\s\}_{\omega=-\pi}^{\pi}$.
We assume that $S_X$ obeys the Paley-Wiener conditions \cite{van_trees:1968}, such that it has a positive entropy-power $0<P_e(S_X)<\infty$, where the entropy power of
a spectrum $S_X$ is defined as $P_e(S_X)\triangleq\frac{1}{2\pi e}e^{2h(X)}$. In the stationary Gaussian case, $h(X)=\frac{1}{2}\log(2\pi e) + \frac{1}{4\pi}\int_{-\pi}^{\pi}
\log(S_X(e^{j\omega})) d\omega$ and the entropy power is therefore given by
\begin{equation}\label{eq:entropy_power}
P_e(S_X) = \exp\bigg( \frac{1}{2\pi}\int_{-\pi}^{\pi}{\log\big(\s\big)} \, d\omega \bigg),
\end{equation}
where here and onwards all logarithms are taken to the natural base
unless explicitly stated otherwise.
The source may be represented in the time-domain by
\begin{equation}\label{eq:Xn}
X[n] = \sum_{k=1}^\infty a_k X[n-k] + I[n],
\end{equation}
or in the $z$-domain by
\begin{equation}
X(z) = \frac{1}{1-A(z)}I(z),
\end{equation}
where  $\{I[n]\}_{n=0}^{\infty}$ is i.i.d.\ zero-mean Gaussian and where
\begin{equation}\label{eq:az}
A(z) \triangleq \sum_{i=1}^{\infty}a_kz^{-k}
\end{equation}
is the \emph{optimal predictor}\footnote{We use the term \emph{optimal
    predictor} to denote the unique filter, which when used to predict
  the source from its infinite past, minimizes the variance of the
  prediction error. If the source is Gaussian, the prediction error,
  i.e., $\{I[n]\}$, is a white Gaussian process~\cite{makhoul:1975}.}
associated with $S_X$.
Using this notation, a spectrum $S_X$ has a spectral decomposition~\cite{makhoul:1975}:
\beq{spectral_decomposition}
\s = \left. \frac{P_e(S_X)}{
(1-A(z))
(1-{A}^*\left(\frac{1}{z^*}\right)) }\right|_{z=e^{j\omega}}. \eeq

\subsection{Symmetric Quadratic-Gaussian MD RDF}
Let $f_i^{(k)}$ be the encoder
for the $i$th description, i.e., $\{Y_i[n]\}_{n=0}^{k-1} =
f_i^{(k)}(\{X[n]\}_{n=0}^{k-1})$ and let
$g_i^{(k)}$ denote the $i$th decoder, i.e., $\{\hat{X}_i[n]\}_{n=0}^{k-1} = g_i^{(k)}(\{Y[n]\}_{n=0}^{k-1})$.
Similarly, let $g_c$ be the joint decoder that produces the central
description $\{\hat{X}_C[n]\}_{n=0}^{k-1}$.

We will use the time-averaged mean squared error (MSE) as fidelity
criterion.

A rate pair $({R}_1,{R}_2)$ is said to be achievable with respect to a
distortion triplet $({D}_1,{D}_2,{D}_C)$ if there exists encoding functions
$f_i^{(k)}, i=1, 2,$ where for all sufficiently
large $k$,
\begin{equation}
\frac{1}{k}\log|f_i^{(k)}|\leq {R}_i, \quad i=1,2,
\end{equation}
and there exists decoding functions $g_c^{(k)}, g_i^{(k)}, i=1,2,$ such that
the side and central distortions satisfy
\begin{equation}
\lim_{k\to \infty}\frac{1}{k}\sum_{n=0}^{k-1}{E}[(X[n] -
\hat{X}_i[n] )^2 ] \leq {D}_i, \quad i=1, 2,
\end{equation}
and
\begin{equation}
\lim_{k\to \infty}\frac{1}{k}\sum_{n=0}^{k-1}\mathbb{E}[(X[n] -
\hat{X}_C[n] )^2 ] \leq D_C.
\end{equation}

In this work, we are interested in the symmetric case,
where $R_1=R_2\triangleq R$ and $D_1=D_2\triangleq D_S$. The symmetric
MD RDF is defined as the minimum rate $R$ per description, which is achievable
with respect to the distortion pair $(D_S,D_C)$. When the source is
white, we replace the spectrum $S_X$ by the variance
$\sigma_X^2$. In the colored case, we refer to $R(\ej)$ as the rate
spectral density, i.e., the rate required for coding $\sx$ at some
frequency $\omega\in [-\pi;\pi]$.
Similarly, we use the notations $D_S(\ej)$ and $D_C(\ej)$ for the side and
central distortion spectral densities, respectively. With a slight
abuse of notation, $D_S(\ej)$ either refers to the complete spectrum
$\{D_S(\ej)\}_{\omega=-\pi}^\pi$ or to the scalar distortion value at a single
frequency $\omega$. Moreover, $D_S$ either refers to the side
distortion or to the side distortion constraint. In all cases, the meaning should be clear from
the context.

A high-resolution MD RDF is the limit where the coding rate diverges
to infinity, and where the side and central distortions both tend to
zero. Thus, under high resolution conditions, we have that $\sx \gg
D_C(\ej)$ and $\sx\gg D_S(\ej)$,
for all $\omega\in[-\pi;\pi]$ where $\sx>0$.

\subsection{Additional Notation}
For $x$ real or complex, $\sqrt[n]{x}$ has $n$ roots.
For $n=2$ and $0\leq x\in\mathbb{R}$ we define $\sqr{x} \triangleq |\sqrt{x}|$, i.e., it is always non-negative. For $0>x\in\mathbb{R}$ we define $\sqr{x} \triangleq i|\sqr{|x|}|$, i.e., we take the principal complex root. For $n=3$ and $x\in \mathbb{R}$ we let $\sqrr{x} \triangleq sign(x)|\sqrr{|x|}|$ denote the unique real cubic root of $x$, e.g., $\sqrr{-8} = -2$. If $x\in\mathbb{C}$ and $imag(x)\neq 0$, we let $\sqrr{x}$ denote the principal complex root, i.e., it has a positive imaginary part.
We use the notation $\xi_i^\Xi$ to indicate the $i$th root of the
function $\Xi$. If $\varphi$ is a function of $\zeta$, we use the
notation $\varphi|_{\zeta = \lambda}$ to indicate that the function
$\varphi$ is evaluated at the point $\zeta=\lambda$. If $\Phi$ is a
matrix, we use $|\Phi|$ to denote the determinant of $\Phi$.

\section{The Quadratic Gaussian Symmetric MD Rate Revisited}\label{sec:qrdf}
In this section we re-state known results about the QG MD achievable rate in the symmetric case, in order to gain
some insight and prepare the ground for what follows. In the high
resolution limit, these results also hold for general sources with finite differential entropy rate \cite{zamir:1999}. We consider four cases: we start with a white source, moving from high resolution to general resolution, and then we repeat the same for colored sources.

\subsection{White Source, High Resolution}
Let us first define the following rate expression
\beq{Ozarow_rate_HR} R_{white,HR}(\sigma_X^2,D_C,D_S) \triangleq
\frac{1}{2}\log\left(\frac{\sigma_X^2}{2\sqrt{D_C(D_S-D_C)}} \right),
\eeq
which is only a
function of the source variance and the side and central distortions
and which is valid when the distortions are non-degenerate, i.e., for distortions satisfying:\footnote{Of
  course, in order for the high-resolution approximation to hold, it
  must be that $D_S\ll \sigma_X^2$. Still, one gets a valid expression
  as long as $D_S\leq\sigma_X^2$.}
\beq{non-degenerate_HR}  2 D_C \leq D_S \leq \sigma_X^2 \eeq
and
\begin{equation}\label{eq:dcbound}
D_C\geq 2D_S-\sigma_X^2.
\end{equation}
In the symmetric rate/side-distortions case
the quadratic-Gaussian MD RDF of Ozarow
\cite{ozarow:1980} converges to \eqref{Ozarow_rate_HR} in the
high-resolution limit where $\sigma_X^2\gg
\max\{D_C,D_S\}$ and when \eqref{non-degenerate_HR} and \eqref{eq:dcbound} are satisfied.

In the high-resolution limit, if the central decoder was to linearly combine two side descriptions
of mutually independent errors of variances $D_S$, it would
achieve exactly the maximal central distortion $D_C=D_S/2$. This
motivates the model of \emph{negatively correlated} side
noises (see \cite{ozarow:1980}). Indeed, in this limit,
the relation between the side and central noises can be
explained by the side noises having a correlation matrix:
\beq{corr1} \Phi = D_S \left[\begin{array}{cc}
            1 & \rho  \\
            \rho & 1 \\
        \end{array} \right],  \eeq
where \[ \rho=-\frac{D_S-2D_C}{D_S}\leq 0. \] With this notation,
\eqref{Ozarow_rate_HR} becomes:
\begin{align}\label{Ozarow_rate_det}
R_{white,HR}(\sigma_X^2,D_C,D_S)
&= \frac{1}{2}\log\left(\frac{\sigma_X^2}{\sqrt{|\Phi|}}\right) \\
&=
\frac{1}{2}\log\left(\frac{\sigma_X^2}{D_S}\right) +
\frac{1}{2}\log\left(\frac{1}{\sqrt{1-\rho^2}} \right) \\
& \triangleq
\frac{1}{2}\log\left(\frac{\sigma_X^2}{D_S}\right) + \frac{1}{2} \delta_{HR},
\end{align}
 where $\delta_{HR}$ is the high-resolution excess
rate \cite{zamir:1999}. Still in the high-resolution case, we take
another step: Without loss of generality, we can represent the
correlated noises as the sum of two mutually independent noises, one
is added to both branches while the other is added to one branch and
subtracted from the other, as depicted in \figref{Ozarow_high}. In the
figure, $\hat{X}_1$ and $\hat{X}_2$ refer to the two side
descriptions, and $\hat{X}_C$ refers to the central reconstruction.

\begin{figure}[t]
\begin{center}

\psfrag{X}{$\scriptstyle X$}
\psfrag{X1}{$\scriptstyle \hat{X}_1$}
\psfrag{X2}{$\scriptstyle \hat{X}_2$}
\psfrag{Xc}{$\scriptstyle \hat{X}_C$}
\psfrag{Z+Z}{$\scriptstyle Z_+ + Z_-$}
\psfrag{Z-Z}{$\scriptstyle Z_+ - Z_-$}
\psfrag{1/2}{$\scriptstyle \ \frac{1}{2}$}
\includegraphics[width=6cm]{figs/ozarow.eps}
\caption{A differential form of Ozarow's double-branch test channel for high resolution coding.}
\label{Ozarow_high}
\end{center}
\end{figure}

Note that the averaging eliminates $Z_-$ from the central
reconstruction $\hat{X}_C$. If we denote the variances of the noises $Z_+$ and
$Z_-$ as $\Theta_+$ and $\Theta_-$, respectively, then we can
re-write \eqref{corr1} as: \beq{corr2} \Phi =
\left[\begin{array}{lr}
            \Theta_++\Theta_- \ & \ \Theta_+-\Theta_- \\
            \Theta_+-\Theta_- \ & \ \Theta_++\Theta_- \\
        \end{array} \right], \eeq with
        \beq{det_phi} |\Phi| = 4 \Theta_+\Theta_-. \eeq
Further, the distortions are related to these variances by:
\beqn{distortions_high}
D_S &=& \Theta_+ + \Theta_-,  \nonumber
\\ D_C &=& \Theta_+. \eeqn
Substituting in \eqref{eq:dcbound}, we find that the non-degenerate
distortion pairs correspond to the region:
\beq{theta_region} 0 \leq  \Theta_+ \leq \Theta_- \leq \frac{\sigma_X^2}{2} \eeq
(which, on account of \eqref{corr2}, confirms that only negative correlation is considered).
\begin{definition}
For variance $\sigma_X^2$ and $2\times 2$ matrix $\Phi$, the $R_0$
function is defined as:
\beq{rate_det} R_0(\sigma_X^2,\Phi) \triangleq
\frac{1}{2} \log \left(\frac {\sigma_X^2}{\sqrt{|\Phi|}}\right). \eeq
\end{definition}
The following proposition gives the high-resolution optimal rate, based on the observations above.
\begin{proposition} \label{prop_white_HR} \textbf{Parametric representation of $\mathbf{R_{white,HR}(\sigma_X^2,D_S,D_C)}$}. For a white Gaussian source of variance $\sigma_X^2$, the high-resolution white MD RDF for \emph{all} non-degenerate distortion pairs is parametrically given by varying $(\Theta_+,\Theta_-)$ over the range \eqref{theta_region} and setting:
\begin{enumerate} \item $R_{white,HR}(\sigma_X^2,D_S,D_C)=R_0(\sigma_X^2,\Phi)$, where $R_0$ is given by \eqref{rate_det} and $\Phi$ is given by \eqref{corr2}. \item $(D_S,D_C)$ according to \eqref{distortions_high}. \end{enumerate}
\end{proposition}

\begin{figure}[t]
\begin{center}
\psfrag{X}{$\scriptstyle X$}
\psfrag{X1}{$\scriptstyle \hat{X}_1$}
\psfrag{X2}{$\scriptstyle \hat{X}_2$}
\psfrag{Xc}{$\scriptstyle \hat{X}_C$}
\psfrag{Z+Z}{$\scriptstyle Z_+ + Z_-$}
\psfrag{Z-Z}{$\scriptstyle Z_+ - Z_-$}
\psfrag{1/2}{$\scriptstyle \ \frac{1}{2}$}
\psfrag{as}{$\scriptstyle \alpha_S$}
\psfrag{ac}{$\scriptstyle \alpha_C$}
\psfrag{V1}{$\scriptstyle V_1$}
\psfrag{V2}{$\scriptstyle V_2$}
\psfrag{Vc}{$\scriptstyle V_C$}
\psfrag{U}{$\scriptstyle U$}

\includegraphics{figs/ozarow_filters.eps}
\caption{Ozarow's test channel with pre and post factors.}
\label{Ozarow}
\end{center}
\end{figure}

\subsection{White Source, General Resolution}

Generalizing our view to all distortion levels, the minimum
achievable symmetric side-descriptions rate was given by Ozarow
\cite{ozarow:1980}:

\begin{equation}\label{Ozarow_rate}
R_{white}(\sigma_X^2,D_C,D_S) \triangleq \frac{1}{4}
\log\left(
\frac{\sigma_X^2(\sigma_X^2-D_C)^2}{4D_C(D_S-D_C)(\sigma_X^2-D_S)} \right)
\end{equation}
for all non-degenerate distortion pairs, which are the ones that satisfy:
\beqn{non-degenerate}
D_S &\leq& \sigma_X^2  \nonumber \\
D_C &\geq& 2D_S - \sigma_X^2 \nonumber \\
 \frac{1}{D_C} &\geq& \frac{2}{D_S} - \frac{1}{\sigma_X^2}
 . \eeqn

A similar correlated-noises
model to \eqref{corr1} can be obtained by expressing $\rho$ in a
rather complicated form. However, we can greatly simplify such an
expression by proper use of pre- and post-factors, as depicted in
\figref{Ozarow} by $\alpha_S$ and $\alpha_C$.
In a point-to-point scenario, it is convenient to make the pre factors
equal to the post factors \cite{gallager:1968},\cite{zamir:1996}. However, this is
generally not possible in MD coding because the optimal post-factors
(Wiener coefficients) are different for the side and central
reconstructions. We choose the pre-factor to be equal to the
\emph{side} post-factor. While this choice seems arbitrary, it will
prove useful in the sequel. Thus we have:
\beqn{alphas} \alpha_S & \triangleq & \sqrt {\frac
{\sigma_X^2 - \Theta_+ - \Theta_-}
{\sigma_X^2}}, \nonumber \\
\alpha_C & \triangleq & \frac {\alpha_S \sigma_X^2}{\alpha_S^2 \sigma_X^2
+ \Theta_+} = \sqrt{\frac{\sigma_X^2(\sigma_X^2 - \Theta_+ -
\Theta_-)}{(\sigma_X^2 - \Theta_-)^2}}. \eeqn Noting that between $U$ and $\{V_1, V_2, V_C\}$ we have exactly the
high-resolution test channel of Proposition~\ref{prop_white_HR}, it is straightforward to show that under this choice, the distortions are given by:
 \beqn{distortions} D_S &=& \Theta_+ +
\Theta_-, \nonumber
\\ D_C &=& \frac{\sigma_X^2 \Theta_+}{\sigma_X^2 - \Theta_-},\eeqn
where as before $\Theta_+$ and $\Theta_-$ are the variances of the
noises $Z_+$ and $Z_-$ in \figref{Ozarow}.
Note that at high resolution conditions $\sigma_X^2 \gg \Theta_-$,
so (\ref{distortions}) reduces to (\ref{distortions_high}). The choice of factors we have made maintains the high-resolution rate expression at any resolution; the following is identical to Proposition~\ref{prop_white_HR}, except for replacing the distortion equations \eqref{distortions_high} by \eqref{distortions}.

\begin{proposition} \label{prop_white} \textbf{Parametric presentation of $\mathbf{R_{white}(\sigma_X^2,D_S,D_C)}$}. For a white Gaussian source of variance $\sigma_X^2$, the white MD RDF for \emph{all} non-degenerate distortion pairs is parametrically given by varying $(\Theta_+,\Theta_-)$ over the range \eqref{theta_region} and setting:
\begin{enumerate} \item $R_{white}(\sigma_X^2,D_S,D_C)=R_0(\sigma_X^2,\Phi)$, where $R_0$ is given by \eqref{rate_det} and $\Phi$ is given by \eqref{corr2}. \item $(D_S,D_C)$ according to \eqref{distortions}. \end{enumerate}
\end{proposition}

\begin{IEEEproof}  To show that it is enough to consider the triangular region \eqref{theta_region}, substitute \eqref{distortions} in  \eqref{non-degenerate} to assert:
\begin{align*}
\Theta_+ + \Theta_- & \leq \sigma_X^2 \nonumber \\
\frac{\sigma_X^2 \Theta_+}{\sigma_X^2 - \Theta_-} & \geq 2(\Theta_+ + \Theta_-) - \sigma_X^2 \\
\frac{\sigma_X^2 - \Theta_-}{\sigma_X^2 \Theta_+} & \geq \frac{2}{\Theta_+ + \Theta_-} - \frac{1}{\sigma_X^2}. \nonumber \end{align*}
Now the second and third equations can be written as:
\begin{align*} \left(\frac{\sigma_X^2}{2}-\Theta_-\right)(\sigma_X^2 - \Theta_+ - \Theta_-) & \geq 0 \nonumber \\
(\Theta_- - \Theta_+) (\sigma_X^2 - \Theta_+ - \Theta_-) & \geq 0, \nonumber  \end{align*} respectively, which, on behalf of the first equation, proves that $\sigma_X^2 / 2 \geq \Theta_- \geq \Theta_+$. Finally, from the first two equations it follows that if $\Theta_+<0$ then also $\Theta_++\Theta_-<0$, but since in this case $D_S<0$, it is not possible. Deriving the rate expression is straightforward, by substituting the distortions \eqref{distortions} in
\eqref{Ozarow_rate}.
\end{IEEEproof}

\subsection{Stationary Source, High Resolution}

We first re-state  the optimal rates in propositions
\ref{prop_white_HR} and \ref{prop_white} for white sources in terms of spectra. This will yield a single expression which is also valid for colored sources. We therefore define the following equivalent distortion spectrum (which plays a part in the time-domain solution even for \emph{white} sources, see \cite{ostergaard:2009}):
\beq{spectrum_def} \tTheta(\ej) \triangleq \left\{
\begin{array}{ll}
            2\Theta_+, & \hbox{$|\omega|< \frac{\pi}{2}$}, \\
            2\Theta_-, & \hbox{$\frac{\pi}{2}\leq  |\omega| < \pi$}. \\
        \end{array} \right.

        \eeq
The
entropy-power~\eqref{eq:entropy_power} of this spectrum is given
by:
\[P_e(\tTheta)  = 2 \sqrt{\Theta_+ \Theta_-} = \sqrt{|\Phi|}.  \]

\begin{definition}
For source spectrum $\sx$ and noise spectrum $\tilde{\Theta}(e^{jw})$,
the $\tilde{R}_0$ function is defined as:
\beq{rate_Pe} \tilde{R}_0(S_X,\tilde\Theta) \triangleq
\frac{1}{2} \log \left(\frac {P_e(S_X)}{P_e(\tTheta)}\right). \eeq
\end{definition}

Since the source is white with spectrum $S_X(\ej)=\sigma_X^2,  \forall \omega$, we have that $P_e(S_X)=\sigma_X^2$.
Consequently, $R_0 = \tilde{R}_0$, where $R_0$ is the optimal rate in
Propositions \ref{prop_white_HR} and \ref{prop_white} and given by
\eqref{rate_det} and $\tilde{R}_0$ is given by \eqref{rate_Pe}.

We now turn to general (colored) stationary Gaussian sources.
In the high resolution limit, it was shown in \cite{zamir:2000} that the
minimum rate is given by Ozarow's rate \eqref{Ozarow_rate_HR} with the
source variance $\sigma_X^2$ replaced by its entropy-power $P_e(S_X)$
\eqref{eq:entropy_power}. Recalling Proposition~\ref{prop_white_HR} and \eqref{rate_Pe}, this leads directly to the following.

\begin{proposition} \label{prop_HR} \textbf{Parametric presentation of $\mathbf{R_{HR}(S_X,D_S,D_C)}$}. For a stationary Gaussian source of spectrum $S_X(\ej)$, the high-resolution MD RDF for \emph{all} non-degenerate distortion pairs is parametrically given by varying $(\Theta_+,\Theta_-)$ over the range \eqref{theta_region} and setting:
\begin{enumerate} \item $R_{HR}(S_X,D_S,D_C)=\tilde{R}_0(S_X,\tTheta)$, where $\tilde{R}_0$ is given by \eqref{rate_Pe} using the spectrum $\tTheta(\ej)$ given by \eqref{spectrum_def}. \item $(D_S,D_C)$ according to \eqref{distortions}. \end{enumerate}
\end{proposition}

\subsection{Stationary Source, General Resolution}

\begin{figure}[t]
\begin{center}
\psfrag{X}{$\scriptstyle X[n]$}
\psfrag{X1}{$\scriptstyle \hat{X}_1[n]$}
\psfrag{X2}{$\scriptstyle \hat{X}_2[n]$}
\psfrag{Xc}{$\scriptstyle \hat{X}_C[n]$}
\psfrag{Z+Z}{\raisebox{-3mm}{$\scriptstyle Z_+[n] + Z_-[n]$}}
\psfrag{Z-Z}{\raisebox{3.5mm}{$\scriptstyle Z_+[n] - Z_-[n]$}}
\psfrag{1/2}{$\scriptstyle \ \frac{1}{2}$}
\psfrag{as}{$\scriptstyle \alpha_S$}
\psfrag{ac}{$\scriptstyle \alpha_C$}
\psfrag{V1}{$\hspace{-2mm} \scriptstyle V_1[n]$}
\psfrag{V2}{$\hspace{-2mm} \scriptstyle V_2[n]$}
\psfrag{Vc}{$\hspace{-2mm} \scriptstyle V_C[n]$}
\psfrag{U}{$\scriptstyle U$}
\psfrag{X}{$\scriptstyle X[n]$}
\psfrag{U}{$\hspace{-1mm} \scriptstyle U[n]$}
\psfrag{B}{$\scriptstyle B[n]$}
\psfrag{D}{$\scriptstyle D[n]$}
\psfrag{Y}{$\scriptstyle Y[n]$}
\psfrag{V}{$\scriptstyle V[n]$}
\psfrag{Z}{$\scriptstyle Z[n]$}
\psfrag{E}{$\scriptstyle E[n]$}
\psfrag{Et}{$\scriptstyle \tilde{E}[n]$}
\psfrag{Xh}{$\scriptstyle \hat{X}[n]$}
\psfrag{F(z)}{$\scriptstyle F(\ej)$}
\psfrag{A(z)}{$\scriptstyle \ A(z)$}
\psfrag{F*(z)}{$\hspace{-2mm} \scriptstyle \ F^{*}(\ej)$}
\psfrag{G*(z)}{$\hspace{-2mm} \scriptstyle \ G^{*}(\ej)$}
\psfrag{C(z)-1}{$\scriptstyle \ \ C(z)$}

\includegraphics{figs/ozarow_filters_new.eps}
\caption{Pre-post filtered test channel for a colored Gaussian source.}
\label{Ozarow_colored}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\psfrag{sx2}{\raisebox{-1mm}{\hspace{-2.5mm}$\frac{\sx}{2}$}}
\psfrag{O+}{$\Tp$}
\psfrag{O-}{$\Tm$}
\psfrag{Ob}{Boundary region}
\psfrag{O}{Region of support}
\includegraphics{figs/boundary1.eps}
\caption{The triangular support region of $\Tm$ and $\Tp$ as given by~(\ref{theta_region_spectra}).}
\label{fig:boundary}
\end{center}
\end{figure}

For general resolution, the achievable colored Gaussian MD rate
region was characterized by Chen et al. \cite{chen:2009}.
We can re-write the result (in the symmetric case) in the spirit of the above exposition. To that end, we present in \figref{Ozarow_colored} a test channel, which may be seen as a colored version of the channel of \figref{Ozarow}. In that channel, the noise processes $Z_+$ and $Z_-$ are stationary Gaussian, with spectra $\{\Theta_+(e^{j2\omega})\}$ and $\{\Theta_{-}(e^{j2\omega})\}$, respectively. In terms of these, we generalize \eqref{spectrum_def} to the form:\footnote{Notice that the
lowpass and highpass spectra of $\tilde{\Theta}$ are
formed by $\{\Theta_+(e^{j2\omega})\}$ and $\{\Theta_{-}(e^{j2\omega})\}$,
which are compressed versions (by a factor of two) of the spectra
$\Theta_+ = \{\Tp\}_{\omega=-\pi}^{\pi}$ and $\Theta_-= \{\Tm\}_{\omega=-\pi}^{\pi}$, respectively.}
\beq{theta}
\tTheta(\ej) = \left\{
\begin{array}{ll}
            2\Theta_+\left(e^{j2\omega}\right), & \hbox{$|\omega| < \frac{\pi}{2}$}, \\
            2\Theta_-\left(e^{j2 \left(\omega-\frac{\pi}{2}\right)}\right), & \hbox{$\frac{\pi}{2} < \omega \leq \pi$}, \\
            2\Theta_-\left(e^{j2 \left(\omega+\frac{\pi}{2}\right)}\right), & \hbox{$-\pi \leq \omega < -\frac{\pi}{2}$}.
        \end{array} \right.
        \eeq
We will be interested in spectra that satisfy \eqref{theta_region} at each frequency, i.e.,
\beq{theta_region_spectra} 0 \leq \Tp \leq \Tm \leq \frac{S_X(\ej)}{2}, \forall \omega. \eeq The per-frequency resulting triangular support region is shown in
Fig.~\ref{fig:boundary}.

The factors $\alpha_S$ and $\alpha_C$ of of \figref{Ozarow} are replaced by filters that satisfy:\footnote{Without loss of generality, one may take all filters to have real frequency response. However, the more general complex form allows more flexibility, e.g. some of the filters may be made causal.} \beqn{filters}
|F(\ej)|^2 &=& \frac{\sx -\Tp - \Tm}{\sx} \nonumber \\
|G(\ej)|^2 &=& \frac{\sx(\sx - \Tp - \Tm)}{(\sx - \Tm)^2} \nonumber \\
\angle{F(\ej)} &=& \angle{G(\ej)}, \eeqn
where $\angle{F(\ej)}$ denotes the phase of $F(\ej)$.
It directly follows, that the
side and central distortion spectral densities are given by: \beqn{distortion_spectra} D_S(e^{j\omega}) &\triangleq& S_{\hat X_1-X}(e^{j\omega}) =  S_{\hat X_2-X}(e^{j\omega}) =  \Tp+\Tm \nonumber \\  \Dc &\triangleq& S_{\hat X_C-X}(\ej)   =  \frac{\sx\Tp}{\sx-\Tm}, \eeqn which reduces to \eqref{distortions} at any fixed frequency $\omega$.
The distortion spectra are connected to the total distortion constraints $(D_S,D_C)$ via:  \beqn{spectra_constraints} \inthalf{\Ds} =  \inthalf{[ \Tp+\Tm]} &\leq& D_S \nonumber \\ \inthalf{\Dc} =  \inthalf{ \frac{\sx\Tp}{\sx-\Tm}} &\leq& D_C . \eeqn In these terms, we have the following:

\begin{proposition}  \label{prop_Chen} \textbf{Parametric presentation of $\mathbf{R(S_X,D_S,D_C)}$}. For a stationary Gaussian source of spectrum $S_X(\ej)$, the RDF for \emph{all} non-degenerate distortion pairs is parametrically given by taking the lower envelope (in the rate axis) of the $(R,D_S,D_C)$ region formed by scanning all spectra $\Tp$ and $\Tm$ satisfying \eqref{theta_region_spectra} and setting:
\begin{enumerate} \item $R=\tilde{R}_0(S_X,\Theta)$, where $\tilde{R}_0$ is given by \eqref{rate_Pe} using the spectrum $\tTheta(\ej)$ given by \eqref{theta}. \item $(D_S,D_C)$ according to \eqref{spectra_constraints}. \end{enumerate}
\end{proposition}

\begin{IEEEproof} In the symmetric case, Theorem 4 of \cite{chen:2009} reduces to:
\[ R(S_X,D_C,D_S)= \min_{\{\Ds\},\{\Dc\}} \inthalf{R_{white}(S_X(\ej),\Dc,\Ds)} , \]
where the minimization is carried over all distortion spectra satisfying \eqref{spectra_constraints}, and where $R_{white}$ is the MD rate for the white problem formed at each frequency $\omega$. Now without loss of generality, one may restrict the minimization to spectra that are everywhere non-degenerate as in \eqref{non-degenerate}, i.e.,  \beqn{non-degenerate_spectra}
\Ds &\leq& S_X(\ej)  \nonumber \\
\Dc &\geq& 2\Ds - S_X(\ej) \nonumber \\
\frac{1}{\Dc} &\geq& \frac{2}{\Ds} - \frac{1}{S_X(\ej)}
\eeqn for all $\omega$ (since if at some frequency this does not hold, we may reduce the contribution of that frequency to the total distortion without increasing its contribution to the total rate). Scanning all spectra satisfying this, amounts to scanning all spectra $\{\Tp\}$ and $\{\Tm\}$ satisfying \eqref{theta_region_spectra}. Now since at any fixed frequency $\omega$, the equation system \eqref{distortion_spectra} is just \eqref{distortions}, we have by Proposition~\ref{prop_white} that \beq{explicit_minimization} R(S_X,D_C,D_S)= \min_{\{\Tp\},\{\Tm\}} \inthalf{\frac{1}{2} \log \frac{S_X(\ej)}{2\sqrt{\Tp\Tm}}} , \eeq where the minimization is taken over all spectra satisfying \eqref{spectra_constraints} and \eqref{theta_region_spectra}. On account of \eqref{theta}, this completes the proof.
\end{IEEEproof}

Unlike Propositions~\ref{prop_white_HR}-\ref{prop_HR}, where any choice of $\Theta$ in the range \eqref{theta_region} led to an optimal rate, here we have an optimization over the spectra.
This result can be seen as follows: Given some distortion spectra, the test-channel at each frequency is the scalar test-channel of Proposition~\ref{prop_white} (it can be shown that this is the \emph{optimal} test channel under distortion-mask constraints, but the proof for that is not included in this work). Then, an optimization needs to be carried out over all distortion spectra such that the total distortion constraints are satisfied. The rate is thus specified as the solution of an optimization problem;
in the next section, we provide an explicit solution for the optimal distortion spectra $\{\Ds\}$ and $\{\Dc\}$, by minimizing \eqref{explicit_minimization} with respect to the noise spectra $\Tp$ and $\Tm$.

\begin{remark} In the high resolution limit, the optimal spectra $\Tp$,
$\Tm$ become flat, thus $\tTheta(\ej)$ becomes a two-step
spectrum, as in Proposition~\ref{prop_HR} and in \cite{ostergaard:2009}. \end{remark}

\begin{remark}\label{rem:epsilon}
If $X$ does not satisfy the Paley-Wiener condition, i.e. $P_e(S_X)=0$, then the rate expression \eqref{rate_Pe} is not well defined.
In this case, we may use the following:
For any $\epsilon>0$, let
$S_{X_\epsilon}(e^{j\omega}) = \max (S_X(e^{j\omega}) , \epsilon), \forall \omega$,
and
$D_\epsilon = \frac{1}{2\pi}\int_{-\pi}^{\pi} \max(0,\epsilon - S_X(e^{j\omega}))  d\omega$.
Then there exists some $\epsilon>0$ such that Proposition~\ref{prop_Chen} holds with
$S_X$, $D_S$, and $D_C$ replaced by
$S_{X_\epsilon}, D_S + D_\epsilon$, and $D_C+D_\epsilon$, respectively.
\end{remark}

\section{Optimal Distortion Spectra}\label{sec:spectral}
In this section we derive a spectral domain characterization of the
optimal distortion spectra.
From Proposition~\ref{prop_Chen} and in particular
from~(\ref{explicit_minimization}), it may be noticed that
finding the RDF is equivalent to finding a pair of noise spectra $\Theta_+=\{\Theta_+(e^{j\omega})\}_{\omega=-\pi}^{\pi}$ and $\Theta_-=\{\Theta_-(e^{j\omega})\}_{\omega=-\pi}^{\pi}$, which minimizes the description rate $R$
subject to the two target distortion constraints $D_S$ and $D_C$, cf.~(\ref{spectra_constraints}).
This constrained minimization problem can also be formulated as a Lagrangian unconstrained problem, which then provides a two-parameter characterization of the RDF.
Specifically, we consider the problem of minimizing the functional
\begin{equation}\label{eq:J}
J = \frac{1}{2 \pi}\int_{-\pi}^{\pi} \L d\omega, \end{equation} where \begin{equation}\label{eq:lagrangian}
\begin{split}
\L&=
\frac{1}{2}\log\left(\frac{\sx}{2\sqr{\Theta_+(e^{j\omega})\Theta_-(e^{j\omega})}}\right)
+
\lambda_1(\Theta_+(e^{j\omega})+\Theta_-(e^{j\omega})) \\
&\quad +
\lambda_2\frac{\sx\Theta_+(e^{j\omega})}{\sx-\Theta_-(e^{j\omega})}.
\end{split}
\end{equation}
The scalar weights $\lambda_1$ and $\lambda_2$ are Lagrangian variables~\cite{weinstock:1974,giaquinta:2004},
which provide a trade-off between rate, side distortion, and central
distortion.

Intuitively, letting $\lambda_1\approx 0,\lambda_2 \approx 0$, results in a rate close to zero since this is equivalent to solving a minimization problem without any constraints (except that of a non-zero rate). On the other hand, letting $\lambda_1 \gg 1$ or $\lambda_2\gg 1$ penalizes one of the distortions and corresponds to a high rate  situation. In particular, if $\lambda_1 \gg \lambda_2$, then the side distortion is severely penalized and is therefore forced to be small. The central distortion is of less concern in this case.  If $\lambda_2\gg \lambda_1$, then the central distortion is minimized and the side distortion is of less concern. Finally, if both $\lambda_1,\lambda_2 \gg 1$, then both the side and central distortions are small.
From the Kuhn-Tucker theorem, it follows that the Lagrangian variables are non-negative, i.e.\ $\lambda_1,\lambda_2\geq 0$, see~\cite{luenberger:1969} for details.
If either $\lambda_1=0, \lambda_2>0$ or $\lambda_1>0,\lambda_2=0$, the
optimization problem reduces to a conventional SD
problem. If $\lambda_1=\lambda_2=0$, we have an unconstrained problem,
with a trivial zero rate solution. Thus, without loss of generality
we may assume that $\lambda_1>0,\lambda_2>0$.

Fortunately, the unconstrained minimization problem~(\ref{eq:J}) may be solved by optimizing the per-frequency function $\L$, as we formally state in the sequel. To that end, we obtain
the following partial derivatives of~(\ref{eq:lagrangian}):
\begin{equation}\label{eq:diftheta+}
\frac{\partial \L}{\partial \Theta_+(e^{j\omega})} =
-\frac{1}{4\Theta_+(e^{j\omega})} + \lambda_1 + \lambda_2\frac{\sx}{\sx-\Theta_-(e^{j\omega})}
\end{equation}
and
\begin{equation}\label{eq:diftheta-}
\frac{\partial \L}{\partial \Theta_-(e^{j\omega})} =
-\frac{1}{4\Theta_-(e^{j\omega})} + \lambda_1 + \lambda_2\frac{\sx\Theta_+(e^{j\omega})}{(\sx-\Theta_-(e^{j\omega}))^2}.
\end{equation}
Equating both~(\ref{eq:diftheta+}) and (\ref{eq:diftheta-}) to zero and jointly solving, yields
\begin{equation}\label{eq:theta2}
\Theta_-(e^{j\omega}) = \Psi_{\lambda_1,\lambda_2}^\dagger(e^{j\omega})
\end{equation}
and
\begin{equation}\label{eq:theta1}
\Theta_+(e^{j\omega}) =
\frac{\sx - \Psi^\dagger_{\lambda_1,\lambda_2}(e^{j\omega})}
{4\sx(\lambda_1+\lambda_2) - 4\lambda_1\Psi^\dagger_{\lambda_1,\lambda_2}(e^{j\omega})},
\end{equation}
where for a fixed pair $(\lambda_1, \lambda_2)\in \mathbb{R}^2$, $\Psi^\dagger_{\lambda_1,\lambda_2}(e^{j\omega})$ denotes a real (and positive) root of the third-order polynomial
\begin{equation}\label{eq:psipoly1}
\begin{split}
\Psi^3(e^{j\omega}) &- \frac{4\lambda_1\lambda_2\sx + 8\lambda_1^2\sx + \lambda_1}{4\lambda_1^2}\Psi^2(e^{j\omega}) \\
&+ \frac{2\lambda_1 \sx + 2\lambda_2 \sx + 4\lambda_1\lambda_2\sxx{2} + 4\lambda_1^2\sxx{2}}{4\lambda_1^2}\Psi(e^{j\omega}) -\frac{\sxx{2}(\lambda_2+\lambda_1)}{4\lambda_1^2}.
\end{split}
\end{equation}
Since~(\ref{eq:psipoly1}) is a real third-order polynomial in $\Psi(e^{j\omega})$, three solutions are possible (of which two might be complex conjugates).
Given a real polynomial $x^3(\ej) +
a_2(\ej)\,x^2(\ej)+a_1(\ej)\,x(\ej)+a_0(\ej)$, where the $a_i(\ej)$'s
follow from~(\ref{eq:psipoly1}), we denote its discriminant $\Xi(\ej)$ by
\begin{align}\label{eq:Xi}
\Xi(\ej) &= q^2(\ej) + p^3(\ej),
\end{align}
where
\begin{equation} \label{eq:p}
\begin{split}
p(\ej) &=\frac{a_1(\ej)}{3} - \frac{a_2^2(\ej)}{9} \\
  &= -\frac{1}{144\lambda_1^2}(-8\lambda_1\sx-16\lambda_2\sx+16\lambda_1\lambda_2\sxx{2} + 16\lambda_1^2\sxx{2} \\
&\quad +16\sxx{2}\lambda_2^2 + 1)
\end{split}
\end{equation}
and
\begin{align}\label{eq:q}
q(\ej) &= \frac{1}{6}(a_1(\ej)\,a_2(\ej)-3a_0(\ej)) -\frac{a_2^3(\ej)}{27} \\ \notag
  &= -\frac{1}{1728\lambda_1^3}\big(
96\lambda_1\lambda_2\sxx{2} - 48\lambda_1^2\sxx{2} - 64\lambda_2^3\sxx{3} - 96\lambda_2^2\sxx{3}\lambda_1 \\
&\quad+ 96\lambda_2^2\sxx{2}  +96\lambda_2\sxx{3}\lambda_1^2 + 24\lambda_2\sx+64\lambda_1^3\sxx{3}+12\lambda_1\sx-1
\big).
\end{align}
Moreover, let $s_1(\ej) = \sqrr{q(\ej) + \sqr{\Xi(\ej)}}$
and $s_2(\ej) = \sqrr{q(\ej) - \sqr{\Xi(\ej)}}$. Then, the
three solutions are given by (see e.g.,~\cite{abramovitz:1973})
\begin{align}\label{eq:x1}
x_1(\ej) &= (s_1(\ej)+s_2(\ej)) -\frac{a_2(\ej)}{3} \\ \label{eq:x2}
x_2(\ej) &= -\frac{1}{2}(s_1(\ej)+s_2(\ej)) -\frac{a_2(\ej)}{3} + \frac{i\sqr{3}}{2}(s_1(\ej)-s_2(\ej)) \\ \label{eq:x3}
x_3(\ej) &= -\frac{1}{2}(s_1(\ej)+s_2(\ej)) -\frac{a_2(\ej)}{3} - \frac{i\sqr{3}}{2}(s_1(\ej)-s_2(\ej)).
\end{align}

At each frequency $\omega$, the global minimum of $\L$ could be one of
these solutions, or lie on the boundary of the region of support given by
\eqref{theta_region_spectra} and
shown
in Fig.~\ref{fig:boundary}.

\begin{lemma}\label{lem:unique_sol}
For any $-\pi<\omega\leq \pi, \lambda_1>0$, and $\lambda>0$, there is at most one unique stationary solution within the
triangular support region \eqref{theta_region_spectra}. If such a
solution exists, it is given by $\Psi(\ej)$, where
\begin{equation}\label{eq:cases_x}
\Psi(\ej) =
\begin{cases}
x_2(\ej), \, \text{if}\ \Xi(\ej)<0,  \\
x_1(\ej), \, \text{if}\ \Xi(\ej)\geq 0,
\end{cases}
\end{equation}
and where $x_1(\ej)$ and $x_2(\ej)$ given by~\eqref{eq:x1} and
\eqref{eq:x2}, respectively, can be explicitly expressed in terms of
$\lambda_1$ and $\lambda_2$ by
\begin{equation}
x_1(\ej) = \sqrr{ q(e^{j\omega})\! +\! \sqr{\Xi(e^{j\omega})} } - \sqrr{ q(e^{j\omega})\! -\! \sqr{\Xi(e^{j\omega})} }
+ \frac{\sx}{3\lambda_1}(2\lambda_1 + \lambda_2)
+\frac{1}{12\lambda_1}
\end{equation}
and
\begin{equation}
x_2(\ej)=-\sqr{|p(e^{j\omega})|}(\cos(\phi(\ej)/3) + \sqr{3}\sin(\phi(\ej)/3)) + \frac{\sx}{3\lambda_1}(2\lambda_1 + \lambda_2)
+\frac{1}{12\lambda_1}.
 \end{equation}
\end{lemma}

The proof is given in Appendix~\ref{app:lem:unique_sol}. To resolve when the global minimum is given by~\eqref{eq:cases_x} and
when it is on the boundary, we define the set
$\Omega(\lambda_1,\lambda_2)$ of \emph{support} frequencies by:
\beq{eq:inband} \Omega(\lambda_1,\lambda_2) \triangleq \{-\pi <\omega\leq\pi : 2\lambda_1 \sx + 8 \lambda_2 \Psi(\ej)  >1\}. \eeq

\begin{lemma}\label{lem:boundary}
For any positive spectrum $S_X$, the triplet ($\omega,\lambda_1,\lambda_2)$, where $ -\pi<\omega\leq
\pi, \lambda_1 > 0$, and  $\lambda_2>0$, yields a stationary solution
$\Psi(\ej)$ when used  in \eqref{eq:cases_x}. If $\omega\in\Omega(\lambda_1,\lambda_2)$, then $\Psi(\ej)$  is
globally optimal. Otherwise, the globally optimal solution is $\Tp=\Tm=\sx/2$.
\end{lemma}
\begin{proof}
For any $\omega\in
[-\pi;\pi]$, the boundary of the triangular optimization region \eqref{theta_region_spectra} (see also Fig.~\ref{fig:boundary}) is given by the union of the following three faces:
\begin{enumerate}
\item $\Tp =0, 0\leq \Tm\leq \sx/2$.
\item $0<\Tp = \Tm \leq \sx/2$.
\item $\Tm = \sx/2, 0<\Tp\leq\Tm$.
\end{enumerate}
We now consider the behavior of $\L$ near these faces.
\begin{enumerate}
\item If $\Tp = 0$, then $R\to \infty$ for any non-negative
$\lambda_1,\lambda_2$ and $\Tm<\sx$. Thus, moving from the boundary
and inside the region can only reduce the function.
\item If $0<\Tp=\Tm<\sx/2$
then $D_S(e^{j\omega}) = 2\Theta_+(e^{j\omega})$ and
the spectral rate is $R(e^{j\omega})=\frac{1}{4\pi}\log\left(\frac{\sx}{D_S(e^{j\omega})}\right)$, which is optimal for the side descriptions, i.e., each side description is on the SD RDF.
The (unnormalized) directional derivative, i.e., the ``normal'' vector
to the surface $\Theta_+(e^{j\omega}) = \Theta_-(e^{j\omega})$, which
points into the support region, is given by the
following difference of partial derivatives:
\begin{equation} \label{eq:directional}
\begin{split}
\frac{\partial \L }{\partial \Tm} &- \frac{\partial \L }{\partial \Tp} \\
&= \frac{\lambda_2 \sx}{\sx-\Tp}\bigg( \frac{2\Tp-\sx}{\sx-\Tp}  \bigg),
\end{split}
\end{equation}
where we have used~(\ref{eq:diftheta+}) and~(\ref{eq:diftheta-}), substituting $\Tp= \Tm)$.
Since $\Tp <\sx/2$ it follows that~(\ref{eq:directional}) is always
negative. Thus, $\L$ can be reduced by going away from the boundary
and into the support region, i.e., by letting $0<\Theta_+(e^{j\omega})< \Theta_-(e^{j\omega})<\sx/2$.
When $\Tm= \Tp=\sx/2$, the spectral distortions satisfy $\Ds = \Dc=\sx$
and the spectral rate is $R(\ej) = 0$.
This is only optimal for $\lambda_1=\lambda_2=0$.
\item For $\Tm=\sx/2$,
the directional derivative is given by
\begin{equation}\label{eq:directional1}
\left.\frac{\partial \L}{\partial \Tm)}\right|_{\Tm=\sx/2} =
-\frac{1}{2\sx} + \lambda_1 +\frac{4\lambda_2\Tp}{\sx},
\end{equation}
which is positive along the whole boundary if and only if $\omega\in\Omega(\lambda_1,\lambda_2)$ \eqref{eq:inband}.
Thus for support frequencies the Lagrangian can
be reduced by going away from the boundary and into the support
region, while for other frequencies it is increased for some boundary points. \end{enumerate}

Now we complete the proof, using that inside the region, $\L$ is
differentiable, and that by Lemma~\ref{lem:unique_sol} it has at most
one extremum, which is given by $\Psi(\ej)$. For the first part of the lemma, $\omega\in\Omega(\lambda_1,\lambda_2)$. The function decreases by moving inwards from all boundary points, thus it must indeed have an extremum, which must be a local and global minimum. We now prove the second part of the lemma. For non-support frequencies the function increases by moving inwards from part of the boundary $\Tm = \sx/2, 0<\Tp\leq\Tm$. Thus, if it has a single extremum inside the region it cannot be a local minimum, and (regardless of the existence of an extremum) the global minimum must lie on that boundary. Now recall that
$\L$ is strictly convex in $\Tp$ for a fixed $0<\Tm<\sx$ and that $\Tp=0$ is a maximum. Moreover,
\begin{equation}\label{eq:di}
\left.\frac{\partial \L}{\partial \Tp)}\right|_{\Tm=\sx/2} =
-\frac{1}{4\Theta_+(e^{j\omega})} + \lambda_1 +2\lambda_2,
\end{equation}
is zero at $\Tp = 1/4/(\lambda_1+2\lambda_2)$ and corresponds to a global minimum of $\mathcal{L}(\Tp,\sx/2)$.
Let $\Tm=\sx/2, \Tp =
 1/4/(\lambda_1+2\lambda_2)$ and solve for $\lambda_2$ when
 equating~(\ref{eq:diftheta-}) to zero, i.e.,
 \begin{equation}
 \lambda_2 = \frac{1}{4\sx} - \frac{\lambda_1}{2},
 \end{equation}
 which implies that $\L$ is minimized for $\Tp = \sx/2$.
Thus, there is a single global minimum at $\Tp=\Tm=\sx/2$, which corresponds to a trivial zero-rate point.
\end{proof}

\begin{remark}
Since the rate spectral density $R(\ej)$, which is given by
\begin{equation}\label{eq:Rw}
R(e^{j\omega})=\frac{1}{2}\log\left(\frac{\sx}{2\sqr{\Theta_+(e^{j\omega})\Theta_-(e^{j\omega})}}\right),
\end{equation}
is zero for $\Tm=\Tp=\sx/2$, it follows that the set of support frequencies $\Omega(\lambda_1,\lambda_2)$~(\ref{eq:inband}) defines
a frequency dependent ``water level'' below which no rate is
allocated. See Figs.~\ref{fig:rdplots} and~\ref{fig:dist_ratios}  for
an example where zero rate is
allocated to a part of a source spectrum.
\end{remark}

We are now in a position to present the symmetric two-description RDF in a parametric form:

\begin{theorem}\label{theo:rdf}
 Let $X$ be stationary Gaussian with spectral density $S_X$, and
 having finite positive differential entropy rate, i.e.,
 $0<\bar{h}(X)<\infty$, and finite variance $\sigma_X^2$. Then
 the symmetric quadratic two-description RDF of $X$, $R(S_X,D_S,D_C)$, is given for all non-degenerate distortion pairs by taking the lower envelope (in the rate axis) of the $(R,D_S,D_C)$ region formed by scanning all $\lambda_1,\lambda_2 \geq 0$ and setting:
\begin{enumerate}
\item $R$ according to:
 \begin{align}\label{eq:rdf}
 R = \frac{1}{4\pi} \int_{-\pi}^\pi \log\left(\frac{\sx}{2\sqr{\Theta_+(e^{j\omega})\Tm}}\right)\, d\omega,
 \end{align}
\item $(D_S,D_C)$ according to \eqref{spectra_constraints}.
\end{enumerate}
In these expressions, for $\omega \in \Omega(\lambda_1,\lambda_2)$,
\begin{equation}
\Tm = \Psi(\ej)
\end{equation}
and
\begin{equation}\label{eq:theta1_opt}
 \Tp=
\frac{\sx - \Psi(\ej)}{4\sx(\lambda_1+\lambda_2) - 4\lambda_1
  \Psi(\ej)},
\end{equation} and otherwise \[ \Theta_-(\ej) = \Theta_+(\ej) = \frac{\sx}{2}. \]
where for any $\lambda_1,\lambda_2> 0$, $\Psi(\ej)$ and $\Omega(\lambda_1,\lambda_2)$  are given by~(\ref{eq:cases_x}) and \eqref{eq:inband}, respectively.
\end{theorem}

\begin{IEEEproof}
The constrained minimization problem of Proposition~\ref{prop_Chen} can be
cast as an extended isoperimetric problem and solved using standard techniques
of variational calculus~\cite{weinstock:1974,giaquinta:2004}. In
particular,
we may form the cost functional shown
in~\eqref{eq:J} and use the fact that the functions $\Tm$ and $\Tp$ which extremizes~(\ref{eq:J}) must also satisfy
\begin{equation}
\frac{\partial \L}{\partial \Tp} =0, \quad \frac{\partial \L}{\partial \Tm} =0,
\end{equation}
for all $-\pi<\omega\leq \pi$, where the Lagrangian $\L$ is given
by~\eqref{eq:lagrangian}. Moreover,
even though the side constraints are defined as integrals of functions, it
 is enough to consider constant Lagrangian multipliers
$(\lambda_1,\lambda_2)$~\cite[Thm.2, p.91]{giaquinta:2004}. Thus,
for a fixed pair of positive Lagrangian variables $(\lambda_1,
\lambda_2)$, the optimal spectra $\{\Tp\}_{\omega=-\pi}^{\pi}$ and
$\{\Tm\}_{\omega=-\pi}^{\pi}$ are determined by
Lemma~\ref{lem:boundary} by sweeping over $-\pi<\omega\leq\pi$. The
RDF at the distortion point $(D_S, D_C)$ is then obtained by choosing
that pair $(\lambda_1,\lambda_2)$, which results in the lowest
rate when used in~\eqref{eq:rdf} and where the corresponding
distortions~\eqref{spectra_constraints} satisfy the constraints $(D_S,D_C)$.
\end{IEEEproof}

To elucidate the behavior of the noise spectra $\Tm$ and $\Tp$ as a function of $\lambda_1$ and $\lambda_2$, we present the following results and examples.

\begin{proposition}[High-Rate Cases]\label{prop:hr}
For any $\omega\in [-\pi;\pi]$,
\begin{equation}\label{eq:tm}
\mathop{\lim_{\lambda_1,\lambda_2\to\infty}}_{\lambda_1/\sqrr{\lambda_2}\to 0}  \lambda_1 \Tm  = \frac{1}{4}
\end{equation}
and
\begin{equation}\label{eq:tp}
\mathop{\lim_{\lambda_1,\lambda_2\to\infty}}_{\lambda_1/\sqrr{\lambda_2}\to 0}  (\lambda_1+\lambda_2) \Tp  = \frac{1}{4}.
\end{equation}
\end{proposition}

\begin{IEEEproof}
See Appendix~\ref{app:proof:prop:hr}.
\end{IEEEproof}

\begin{remark}
The convergence requirement of $\lambda_1/\sqrr{\lambda_2}\to 0$ in Proposition~\ref{prop:hr} is a technicality needed in the proof.
As shown in Fig.~\ref{fig:highrate}, the limiting behavior of
(\ref{eq:tm}) and (\ref{eq:tp}) can also be observed for small
$\lambda_2$ and large $\lambda_1$. This shows that under
high-resolution conditions, the optimum noise spectra are flat,
independent of the source spectrum, and approximately given by $\Tm\approx \frac{1}{4\lambda_1}$ and $\Tp\approx \frac{1}{4(\lambda_1+\lambda_2)}$.
\end{remark}

\begin{figure}
\psfrag{lambda1}{\scriptsize$\lambda_1$}
\psfrag{Theta+}{$\scriptscriptstyle(\lambda_1+\lambda_2)\Tp$}
\psfrag{Theta-ThetaThetaThetaTheta}{$\scriptscriptstyle\lambda_1\Tm$}
\begin{center}
\includegraphics[width=8cm]{figs/highrate.eps}
\caption{High-rate convergence of $(\lambda_1+\lambda_2)\Tp$ and $\lambda_1\Tm$ as $\lambda_1\to \infty$ when $\sx = 1, \lambda_2=2$ and any $\omega$. }
\label{fig:highrate}
\end{center}
\end{figure}

\begin{example}
Let the source have a positive and monotonically decreasing spectrum given by
\begin{equation}\label{eq:sx}
\sx = \cos(\omega)+1,\quad 0\leq |\omega| < \pi,
\end{equation}
and shown in Fig.~\ref{fig:rdplots}.
Moreover, let the distortion constraints be $D_C = 0.08$ and $D_S= 0.4$.

\begin{figure}[th]

\psfrag{w}{\scriptsize$\omega$}
\psfrag{Sx}{\scriptsize$\sx$}
\psfrag{Theta-}{\scriptsize$\Theta_-(e^{j\omega})$}
\psfrag{Theta+Theta}{\scriptsize$\Theta_+(e^{j\omega})$}
\psfrag{Dc}{\scriptsize$D_C(e^{j\omega})$}
\psfrag{Ds}{\scriptsize$D_S(e^{j\omega})$}
\psfrag{Rw}{\scriptsize$R(e^{j\omega})$}
\psfrag{pi/3}{\scriptsize\hspace{0.5mm}$\frac{\pi}3$}
\psfrag{2/3 pi}{\scriptsize\hspace{0.5mm}$\frac{2\pi}3$}
\psfrag{pi}{\scriptsize$\pi$}
\begin{center}
\includegraphics[width=8cm]{figs/dist_spectra.eps}
\caption{Source $S_X$, noise $\Theta_+, \Theta_-$, and distortion $\{D_S(e^{j\omega})\}$, $\{D_S(e^{j\omega})\}$ spectra for $0\leq |\omega|<\pi$.}
\label{fig:rdplots}
\end{center}
\end{figure}

Using the closed-form expressions for the noise spectra provided by Theorem~\ref{theo:rdf},  we have numerically performed a simple grid search over $\lambda_1$ and $\lambda_2$. As $\lambda_1$ and $\lambda_2$ are varied, we compared the resulting side and central distortions given by~(\ref{spectra_constraints}) to the above mentioned distortion contraints $D_S$ and $D_C$.
The noise spectra that resulted in distortions closest to the constraints are shown in Fig~\ref{fig:rdplots}. The spectra were obtained with $\lambda_1=0.2380$, and $\lambda_2 = 2.700$, which resulted in
$D_C = 0.0801$ and $D_S= 0.4000$. Moreover, when using these spectra in~(\ref{eq:rdf}) the obtained per description rate is $R=0.7468$ bits/dim.
In Fig.~\ref{fig:rdplots}, we have also shown the resulting side and central distortion spectra.

To better illustrate the trade-off between central and side
distortions as a function of the source spectrum, we have shown the
ratio $\log(\Theta_-(e^{j\omega})/\Theta_+(e^{j\omega}))$ as well as
the ratio $\log(D_S(e^{j\omega})/D_C(e^{j\omega}))$ in
Fig.~\ref{fig:dist_ratios}. Also shown in Fig.~\ref{fig:dist_ratios},
is the sum-rate $2R(e^{j\omega})$ allocated to each frequency band,
where $R(e^{j\omega})$ denotes the per description rate spectral density~(\ref{eq:Rw}).
It may be noticed that zero rate is allocated for the part of the source spectrum, which lies below a certain threshold (as is also the case in conventional SD reverse water-filling).

\begin{figure}[th]

\psfrag{Theta-/Theta+Theta+ThetaTheta}{\scriptsize$\log(\Theta_-(e^{j\omega})/\Theta_+(e^{j\omega}))$}
\psfrag{Ds/Dc}{\scriptsize$\log(D_S(e^{j\omega})/D_C(e^{j\omega}))$}
\psfrag{Rw}{\scriptsize$2R(e^{j\omega})$}
\psfrag{pi/3}{\scriptsize\hspace{0.5mm}$\frac{\pi}3$}
\psfrag{2/3 pi}{\scriptsize\hspace{0.5mm}$\frac{2\pi}3$}
\psfrag{pi}{\scriptsize$\pi$}
\psfrag{w}{\scriptsize$\omega$}
\begin{center}
\includegraphics[width=8cm]{figs/dist_ratio.eps}
\caption{Noise spectra ratio $\log(\Theta_-(e^{j\omega})/\Theta_+(e^{j\omega}))$, and distortion spectra ratio $\log(D_S(e^{j\omega})/D_S(e^{j\omega}))$ as a function of the source spectrum given by~(\ref{eq:sx}).
Also shown is the sum-rate spectral density $2R(e^{j\omega})$ (in this latter case, the y-axis represents bits/dim.\ instead of distortion ratios).}
\label{fig:dist_ratios}
\end{center}
\end{figure}
\end{example}

\begin{figure}[th]
\psfrag{Ds-lambda2}{\scriptsize$D_S(e^{j\omega}), \lambda_1=3$}
\psfrag{Ds-lambda1}{\scriptsize$D_S(e^{j\omega}), \lambda_2=3$}
\psfrag{Dc-lambda2}{\scriptsize$D_C(e^{j\omega}), \lambda_1=3$}
\psfrag{Dc-lambda1lambda1}{\scriptsize$D_C(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda1}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda2-}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=3$}
\psfrag{Theta+(w)-lambda1}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta+(w)-lambda2++}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=3$}
\psfrag{lambda1,lambda2}{\hspace{5mm}\scriptsize$\lambda_1 \text{\ or\
  } \lambda_2$}
\psfrag{R(lambda1)-lambda}{\scriptsize$R(e^{j\omega}), \lambda_2=3$}
\psfrag{R(lambda2)}{\scriptsize$R(e^{j\omega}), \lambda_1=3$}
\psfrag{R(omega)}{\scriptsize$R(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{R}{\scriptsize$R(e^{j\omega})$}
\psfrag{Theta-(w)-lambda}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Theta+(w)-lambda}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Ds-lambda}{\scriptsize$D_S(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Dc-lambda}{\scriptsize$D_C(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{10log10() [dB]}{\scriptsize$10\log_{10}(\cdot)$ [dB]}
\begin{center}
\includegraphics[width=8cm]{figs/thetas.eps}
\caption{$\Theta_-(e^{j\omega})$ and $\Theta_+(e^{j\omega})$ as functions of $\lambda_1$ and as functions of $\lambda_2$. The circles indicate when $\Xi$ switches from being positive to become negative, which corresponds to the cases given by~(\ref{eq:cases_x}).
A single frequency bin where $\sx=2$ is considered.}
\label{fig:theta_spectra}
\end{center}
\end{figure}

\begin{figure}[th]
\psfrag{Ds-lambda2}{\scriptsize$D_S(e^{j\omega}), \lambda_1=3$}
\psfrag{Ds-lambda1}{\scriptsize$D_S(e^{j\omega}), \lambda_2=3$}
\psfrag{Dc-lambda2}{\scriptsize$D_C(e^{j\omega}), \lambda_1=3$}
\psfrag{Dc-lambda1lambda1+}{\scriptsize$D_C(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda1}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda2}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=3$}
\psfrag{Theta+(w)-lambda1}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta+(w)-lambda2}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=3$}
\psfrag{lambda1,lambda2}{\hspace{5mm}\scriptsize$\lambda_1  \text{\ or\
  }\lambda_2$}
\psfrag{R(lambda1)-lambda}{\scriptsize$R(e^{j\omega}), \lambda_2=3$}
\psfrag{R(lambda2)}{\scriptsize$R(e^{j\omega}), \lambda_1=3$}
\psfrag{R(omega)}{\scriptsize$R(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{R}{\scriptsize$R(e^{j\omega})$}
\psfrag{Theta-(w)-lambda}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Theta+(w)-lambda}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Ds-lambda}{\scriptsize$D_S(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Dc-lambda}{\scriptsize$D_C(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{10log10() [dB]}{\scriptsize$10\log_{10}(\cdot)$ [dB]}
\begin{center}
\includegraphics[width=8cm]{figs/dcds.eps}
\caption{$D_C(e^{j\omega})$ and $D_S(e^{j\omega})$ as functions of $\lambda_1$ and as functions of $\lambda_2$. A single frequency bin where $\sx=2$ is considered.}
\label{fig:DsDc_spectra}
\end{center}
\end{figure}

\begin{figure}

\psfrag{Ds-lambda2}{\scriptsize$D_S(e^{j\omega}), \lambda_1=3$}
\psfrag{Ds-lambda1}{\scriptsize$D_S(e^{j\omega}), \lambda_2=3$}
\psfrag{Dc-lambda2}{\scriptsize$D_C(e^{j\omega}), \lambda_1=3$}
\psfrag{Dc-lambda1lambda1}{\scriptsize$D_C(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda1}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda2}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=3$}

\psfrag{lambda1,lambda2}{\hspace{5mm}\scriptsize$\lambda_1 \text{\ or\
  } \lambda_2$}
\psfrag{R(lambda1)-lambda}{\scriptsize$R(e^{j\omega}), \lambda_2=3$}
\psfrag{R(lambda2)}{\scriptsize$R(e^{j\omega}), \lambda_1=3$}
\psfrag{R(omega)}{\scriptsize$R(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{R}{\scriptsize$R(e^{j\omega})$}
\psfrag{Theta-(w)-lambda}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Theta+(w)-lambda++++}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Ds-lambda}{\scriptsize$D_S(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Dc-lambda}{\scriptsize$D_C(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{10log10() [dB]}{\scriptsize$10\log_{10}(\cdot)$ [dB]}
\begin{center}
\includegraphics[width=8cm]{figs/thetadcds.eps}
\caption{$\Theta_-(e^{j\omega}), \Theta_+(e^{j\omega}), D_C(e^{j\omega})$, and $D_S(e^{j\omega})$ as
functions of $\lambda_1=\lambda_2$. A single frequency bin $\omega$ where $\sx=2$ is considered.}
\label{fig:ThetaDc}
\end{center}
\end{figure}

\begin{figure}

\psfrag{Ds-lambda2}{\scriptsize$D_S(e^{j\omega}), \lambda_1=3$}
\psfrag{Ds-lambda1}{\scriptsize$D_S(e^{j\omega}), \lambda_2=3$}
\psfrag{Dc-lambda2}{\scriptsize$D_C(e^{j\omega}), \lambda_1=3$}
\psfrag{Dc-lambda1lambda1}{\scriptsize$D_C(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda1}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta-(w)-lambda2}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=3$}
\psfrag{Theta+(w)-lambda1}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_2=3$}
\psfrag{Theta+(w)-lambda2}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=3$}
\psfrag{lambda1,lambda2}{\hspace{5mm}\scriptsize$\lambda_1 \text{\ or\
  }\lambda_2$}
\psfrag{R(lambda1)-lambda+}{\scriptsize$R(e^{j\omega}), \lambda_2=3$}
\psfrag{R(lambda2)}{\scriptsize$R(e^{j\omega}), \lambda_1=3$}
\psfrag{R(omega)}{\scriptsize$R(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{R}{\scriptsize$R(e^{j\omega})$}
\psfrag{Theta-(w)-lambda}{\scriptsize$\Theta_-(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Theta+(w)-lambda}{\scriptsize$\Theta_+(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Ds-lambda}{\scriptsize$D_S(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{Dc-lambda}{\scriptsize$D_C(e^{j\omega}), \lambda_1=\lambda_2$}
\psfrag{R [bits/dim.]}{\scriptsize$R$ [bits/dim.]}
\begin{center}
\includegraphics[width=8cm]{figs/rates.eps}
\caption{$R(e^{j\omega})$ as a function of $\lambda_1$, as a function of $\lambda_2$, and as a function of $\lambda_1=\lambda_2$.
A single frequency bin where $\sx=2$ is considered.}
\label{fig:Rw}
\end{center}
\end{figure}

\begin{example}
For a given $\omega \in [-\pi;\pi]$ let $\sx = 2$ and let $\lambda_1=3$. Then the noise spectral density components $\Theta_-(e^{j\omega})$ and $\Theta_+(e^{j\omega})$ are shown in Fig.~\ref{fig:theta_spectra} (expressed in dB) as a function of $\lambda_2$ in the range $[0;10]$.
For $\lambda_2 < 0.52$, the discriminant $\Xi$ is positive whereas for $\lambda_2>0.52, \Xi$ is negative.
We have indicated the switching point with circles in the figure.
Notice that already at e.g., $\lambda_2=4$,  the approximation $\frac{1}{4(\lambda_1+\lambda_2)}$ of~(\ref{eq:tp}), i.e., $-10\log_{10}(4(\lambda_1+\lambda_2))= -14.47$ dB provides a good approximation of $\Theta_+(e^{j\omega})$.
Also shown in Fig.~\ref{fig:theta_spectra}, are $\Theta_-(e^{j\omega})$ and $\Theta_+(e^{j\omega})$ as a function of $\lambda_1\in [0;10]$ for fixed $\lambda_2=3$ and $\sx=2$.
In this case, $\Xi<0$ for all $\lambda_1$. Notice that at $\lambda_1=4$, the approximation $\frac{1}{4\lambda_1}$ of~(\ref{eq:tm}), i.e., $-10\log_{10}(4\lambda_1) = -12.04$ and~(\ref{eq:tp}), i.e., $-10\log_{10}(4(\lambda_1+\lambda_2))= -14.47$, provide good approximations of $\Theta_-(e^{j\omega})$ and $\Theta_+(e^{j\omega})$, respectively. As $\lambda_1\to 0$, $\Theta_-(e^{j\omega})\to \sx/2$, and since $\lambda_2>0$, the rate is used entirely for reducing the central distortion. However, as $\lambda_1$ increases, more rate is spend on decreasing $\Theta_-(e^{j\omega})$ with less emphasis on $\Theta_+(e^{j\omega})$.

The side and central distortions $\Ds$ and $\Dc$ for the above $\Tm$ and $\Tp$ and given $\omega$ are shown in Fig.~\ref{fig:DsDc_spectra}.
In Fig.~\ref{fig:ThetaDc}, we have illustrated $\Theta_-(e^{j\omega}), \Theta_+(e^{j\omega}), D_S(e^{j\omega})$, and $D_C(e^{j\omega})$ as functions of $\lambda_1=\lambda_2$, and the corresponding description rate spectral densities $R(e^{j\omega})$, given by~(\ref{eq:Rw}), are shown in Fig.~\ref{fig:Rw}.

\end{example}

\section{Time-Domain SD Source Coding Subject to a Distortion Mask}\label{sec:mask}

We take a detour to a SD problem that is suggested by
Proposition \ref{prop_Chen}; coding of a source subject to a maximum
distortion \emph{mask} $\{D(\ej)\}$, rather than subject to a total
distortion constraint. That is, we wish to compress a stationary Gaussian source $X$ such that the spectrum of the reconstruction error satisfies:\footnote{In principle, there is no guarantee that the reconstruction error is stationary. The condition can be thought of this way: pass the reconstruction error through a bank of filters, and measure the MSE as in the regular definition in Section~\ref{sec:prelim}; now take the limit of narrow filters.} \beq{distortion_mask} S_{\hat X - X}(\ej) \leq D(\ej), \forall -\pi < \omega \leq \pi . \eeq

Without loss
of generality\footnote{Otherwise, there is just wasted allowed
distortion which does not serve to reduce the rate.}, we assume that
$D(\ej)\leq \s, \forall \omega$. It is easy to verify, that the
minimum rate for this problem is given by (recall
\eqref{rate_Pe}): \beq{RD_mask}
R\big(S_X,D\big) = \frac{1}{2}
\log\left(\frac{P_e(S_X)}{P_e\big(D\big)}\right). \eeq

\begin{figure*}[t]
\begin{center}

\psfrag{X}{$\scriptstyle X[n]$}
\psfrag{U}{$\scriptstyle U[n]$}
\psfrag{B}{$\scriptstyle B[n]$}
\psfrag{D}{$\scriptstyle D[n]$}
\psfrag{Y}{$\scriptstyle Y[n]$}
\psfrag{V}{$\scriptstyle V[n]$}
\psfrag{Z}{$\scriptstyle Z[n]$}
\psfrag{E}{$\scriptstyle E[n]$}
\psfrag{Et}{$\scriptstyle \tilde{E}[n]$}
\psfrag{Xh}{$\scriptstyle \hat{X}[n]$}
\psfrag{F(z)}{$\scriptstyle F(z)$}
\psfrag{A(z)}{$\scriptstyle \ A(z)$}
\psfrag{F*(z)}{$\scriptstyle \ F^{*}(z)$}
\psfrag{C(z)-1}{$\scriptstyle \ \ C(z)$}
\includegraphics{figs/nested_channel.eps}
\caption{A DSQ/DPCM equivalent channel for SD coding subject to a
distortion mask.} \label{DPCM_DSM}
\end{center}
\end{figure*}

\figref{DPCM_DSM} presents a \emph{time-domain} test channel which
achieves this rate. Motivated by the ratio of entropy powers
\eqref{RD_mask}, we strive to achieve the optimal rate by the
combination of \emph{source prediction} in order to produce a prediction error of power $P_e(S_X)$, and \emph{noise
shaping} in order to shape the white quantization noise of power
$P_e(D)$ into the spectrum $D$.\footnote{An alternative
time-domain approach, is to accommodate for the distortion mask by
changing the pre and post-filters. However, we choose the
noise-shaping approach for the sake of extending this test channel to the
MD setting.}
These two tasks, we perform by a DPCM loop
\cite{zamir:2008} and a noise-shaping loop
\cite{ostergaard:2009}, respectively. In this test channel, $Z[n]$ is AWGN
of variance $P_e(D)$ and $A(z)$, which is given by~(\ref{eq:az}), is the optimal predictor
of the source spectrum $S_X$.\footnote{We assume that the optimal predictor $A(z)$ for the source spectrum exists. If not, then we may use the procedure outlined in Remark~\ref{rem:epsilon} in order to construct a predictor, which satisfies the assumption.} Moreover, the noise-shaping filter is given by: \begin{equation}\label{eq:C(z)}
C(z) = \frac{Q(z)}{1-Q(z)}.
\end{equation}
where $Q(z)$
is the optimal predictor \eqref{eq:az} for the distortion mask $\{D(\ej)\}$, i.e., for $z=\ej$,
\beq{D-Q-C}
D(\ej) =    \frac{P_e(D)}{\left|1-Q(\ej)\right|^2}
= P_e(D) \left| 1+C(\ej) \right|^2.
\eeq

Note that
$E[n]$, the input to the noise-shaping filter, is equal to $Z[n]$.
The pre-filter $F(\ej)$ satisfies: \beq{pre_filter}|F(\ej)|^2 =
\frac{\s-D(\ej)}{\s}. \eeq

\begin{theorem} \label{thm_DPCM_DSM} For all source and distortion spectra satisfying $S_X(\ej)\geq D(\ej), \forall \omega$ and $\infty>P_e(S_X)\geq P_e(D)>0$, the channel of
\figref{DPCM_DSM} with the choices above, satisfies:
\beq{DPCM_DSM_dist} S_{\hat X - X}(\ej) = D(\ej), \quad -\pi \leq \omega \leq \pi, \eeq
with the \emph{scalar} mutual information
$I(D[n];Y[n])=R\big(S_X,D\big)$ of \eqref{RD_mask}.
\end{theorem}

\begin{IEEEproof} \[V[n]-B[n] = Y[n]-D[n] = Z[n], \] thus it is a white process with variance $P_e(D)$. Now \[ V[n]-U[n] = ( V[n]-B[n] ) * (1+c_n) \] thus by \eqref{D-Q-C} \[ S_{V-U}(\ej) = D(\ej). \] By the choice of $F(\ej)$ this is also the spectrum of $\hat X[n] - X[n]$, which completes the proof of the first claim. For the second claim, straightforward analysis of the scheme shows that \[Y[n] = (1-a_n)*(U[n]+(1+c_n)*Z[n]). \] Using the independence of $\{Z[n]\}$ from $\{U[n]\}$ and the choice of $F(\ej)$, $A(z)$ and $C(z)$, one finds that $Y[n]$ is white with variance $P_e(S_X)$. On account of the Gaussian channel capacity, the scalar mutual information $I(D[n];Y[n])$ follows (the \emph{scalar} mutual information is not changed by the channel input $D[n]$ not being white). \end{IEEEproof}

\begin{remark}
The above proof does not use information-theoretic arguments, and does
not offer intuition as to why the structure or the choice of filters
are good. The same result can also be derived using an alternate proof, which shows that the \emph{instantaneous} mutual information over the AWGN from $D[n]$ to $U[n]$ must equal the mutual information \emph{rate} between the processes $U[n]$ and $V[n]$.
The interested reader is referred to \cite{zamir:2008} and
\cite{ostergaard:2009}  for examples of such proof.

\end{remark}

\begin{remark}
In the special case of a white distortion mask
$D$, the constraint becomes (by the water-filling principle)
equivalent to a regular quadratic distortion constraint. Indeed, the
channel collapses in this case to the pre/post filtered DPCM channel
of \cite{zamir:2008}. Much of the analysis there remains valid for
this problem as well. In particular, we can construct an optimal
coding scheme using the test channel of \figref{DPCM_DSM}, substituting the AWGN for an ECDQ,
and the scalar mutual information $I(D[n];Y[n])$ is also equal to
the directed mutual information $I(D[n]\rightarrow Y[n])$. See the end of Section~\ref{sec:mdscheme} for more details.
\end{remark}

\section{Optimal Time-Domain Colored MD Coding}\label{sec:mdscheme}

We now turn back to the original colored MD problem. Proposition~\ref{prop_Chen} tells us that if we optimize for the spectrum $\tTheta(\ej)$, then the test-channel of \figref{Ozarow_colored}  describes the optimal $R(S_X,D_S,D_C)$. However, we would like a test channel that provides a recipe for a coding scheme, i.e., the additive noise elements should be replaced by quantizers. Unfortunately, the noises in \figref{Ozarow_colored} are colored, and furthermore they are correlated, making the task hard. In this section we present a test-channel with a single white additive noise, in the spirit of the SD channel of Section~\ref{sec:mask}. After proving that the channel gives the ``right'' rate-distortion tradeoff, we present a corresponding coding scheme.

\subsection{Test Channel}

\begin{figure*}[t]
\begin{center}
\psfrag{X}{$\scriptstyle X[n]$}
\psfrag{U}{$\scriptstyle U[m]$}
\psfrag{B}{$\scriptstyle B[m]$}
\psfrag{D}{\hspace{-0.5mm}$\scriptstyle D[m]$}
\psfrag{Y}{\hspace{-0.5mm}$\scriptstyle Y[m]$}
\psfrag{V}{$\scriptstyle V[m]$}
\psfrag{Z}{$\scriptstyle Z[m]$}
\psfrag{E}{\hspace{-1mm}$\scriptstyle E[m]$}
\psfrag{Et}{\hspace{-1mm}$\scriptstyle \tilde{E}[m]$}
\psfrag{Xh}{$\scriptstyle \hat{X}_S[n]$}
\psfrag{Xc}{$\scriptstyle \hat{X}_C[n]$}
\psfrag{F(z)}{$\scriptstyle F(z)$}
\psfrag{A(z)}{$\scriptstyle A(z^2)$}
\psfrag{F*(z)}{$\scriptstyle \ F^{*}(z)$}
\psfrag{G*(z)}{$\scriptstyle \ G^{*}(z)$}
\psfrag{C(z)-1}{$\scriptstyle \ \ C(z)$}
\psfrag{LPF}{\scriptsize \ LPF}

\includegraphics[width=\textwidth]{figs/nested_channel_filters.eps}
\caption{A DSQ/DPCM equivalent channel for MD coding of a colored source. We use the index $n$
for sequences which are ``running'' at the source rate, and the
index $m$ when referring to the upsampled rate.}
\label{DPCM_DSM_2}
\end{center}
\end{figure*}

\figref{DPCM_DSM_2} shows the adaptation of the distortion-mask
equivalent channel to the MD problem. Following
\cite{ostergaard:2009}, we combine interpolation by a factor of two (upsampling and perfect low-pass filtering)
with the noise-shaping loop, forming a DSQ loop. $A(z)$ is the optimal predictor \eqref{eq:az} of
the spectrum $S_X(\ej)$.
It is worth pointing out that the predictor is applied
on the noisy output of the quantizer. Thus, the DPCM loop uses pure
predictions only in the limit of high resolutions,
see~\cite{zamir:2008} for details.
Note also that we apply an upsampled version of the source predictor, namely
$A(z^2)$.
The DSQ
loop, on the other hand, works in the upsampled rate and the
noise-shaping filter $C(z)$ is given by~(\ref{eq:C(z)}) for the
spectrum $\tTheta(\ej)$ given by \eqref{theta}. Recall from
\eqref{theta} that $\Theta_+(\omega)$ and $\Theta_-(\omega)$ now play
the roles of ``in-band noise'' and ``out of band noise''.
In the high-resolution case, the DSQ is independent of the shape of
the source
spectrum and uses a two-step spectrum as in
the white case, see \cite{ostergaard:2009} for details.
The additive Gaussian noise $Z[m]$ is white with variance \[P_e(\tTheta) = 2 \sqrt{P_e(\Theta_+)P_e(\Theta_-)}. \]
Finally, the filters $F(\ej)$ and $G(\ej)$ are chosen according to \eqref{filters}. The channel is completely characterized, then, by the choice of $\Tp$ and $\Tm$, which can be performed e.g. according to the optimization in Section~\ref{sec:spectral}. We show that regardless of the optimization, for any choice of spectra that satisfy \eqref{theta_region_spectra}, we have a rate (recall \eqref{explicit_minimization}):
\beq{before_minimization} R(S_X,\Theta_+,\Theta_-)  = \frac{1}{2} \log \left(\frac{P_e(S_X)}{2\sqrt{P_e(\Theta_+)P_e(\Theta_-)}}\right) =  \inthalf{\frac{1}{2} \log \frac{S_X(\ej)}{2\sqrt{\Tp\Tm}}}. \eeq

\begin{remark} For a white
source, $A(z)=0$ and the channel reduces to the DSQ MD test channel of
\cite{ostergaard:2009}, while for optimal side distortion, $C(z)=0$,
and the channel reduces to an upsampled version of the DPCM
equivalent channel of \cite{zamir:2008}. \end{remark}

\begin{theorem} \label{thm_DPCM_DSM_2} The channel of
\figref{DPCM_DSM_2}, with the above choice of filters, satisfies:
\beqn{DPCM_DSM_2_dist} S_{\hat X_C - X}(\ej) & = & \Dc , \quad -\pi \leq \omega \leq \pi,
\nonumber \\ S_{\hat X_S - X}(\ej) & = & \Ds , \quad -\pi \leq \omega \leq \pi,
\eeqn where the distortion spectra were defined in
(\ref{distortion_spectra}), while the \emph{scalar} mutual
information $I(D[m];Y[m])$ equals the rate \eqref{before_minimization}.
\end{theorem}

\begin{IEEEproof}
The channel between $U[m]$ and $V[m]$ is identical to the one between $U[n]$ and $V[n]$ in \figref{DPCM_DSM}. We can thus use the proof of Theorem~\ref{thm_DPCM_DSM} to assert:
\beq{eq:m1} U[m] - V[m] =  (1+c_m) * Z[m] \eeq  \beq{eq:m2} Y[m] = (1-\tilde a_m)*(U[m]+(1+c_m)*Z[m]), \eeq where $\tilde a_m$ corresponds to $A(z^2)$, i.e., equals $a_{m/2}$ for even $m$ and zero otherwise.
For the distortion spectra, \eqref{eq:m1} means that $S_{V-U}(\ej) = \tTheta(\ej)$. From here, we use the structure of rate conversion and low-pass filtering, in conjunction with \eqref{theta}, to obtain:
\beqn{eq:n} S_{\hat X_C - X}(\ej) &=& \frac{1}{2}  \tTheta(e^{j2\omega}) = \Tp \nonumber \\ S_{\hat X_C - X}(\ej) &=& \frac{1}{2} \left[ \tTheta(e^{j2\omega}) + \tTheta(e^{j2\omega - j \pi})\right] = \Tp + \Tm. \eeqn By the choice of $F(\ej)$ and $G(\ej)$, the distortion spectra follow. For the mutual information, using the independence of $\{Z[n]\}$ from $\{U[n]\}$ and the choice of $F(\ej)$, $A(z)$ and $C(z)$, one finds that $Y[n]$ is white with variance $P_e(S_X)$. On behalf of the noise variance and the AWGN channel capacity, the proof is completed.
\end{IEEEproof}

\subsection{Coding Scheme}

We now present a coding scheme based on the optimal time-domain test channel. The underlying principle is that an additive Gaussian white noise element in the test-channel may be replaced by an entropy-coded dithered quantizer (ECDQ), such that the rate of the quantizer is the scalar mutual information over the AWGN, and the variance of the (additive) quantization noise is the AWGN variance. This has been shown in~\cite{zamir:1996}, and extended to test channels with feedback loops in~\cite{zamir:2008}. We give a short account of these results, before turning to the MD scheme.

At all resolutions, an ECDQ can approach, at the limit of high lattice dimension $K\rightarrow\infty$,
a rate
\[ R = \frac{1}{2} \log \left(\frac{\sigma_Y^2}{\sigma_Z^2}\right), \] where $\sigma_Y^2$ and $\sigma_Z^2$ are the output and quantization noise variances, respectively. Furthermore, the quantization noise is independent of the quantizer input.\footnote{In the high-dimensional limit, the
quantization noise becomes approximately Gaussian distributed (in a
divergence sense) \cite{zamir:1996a}. However, for achievability results based on linear operations, such as the ones that appear in this paper, we do not need to rely on Gaussianity.}
Thus, in conjunction with optimal factors (i.e. Wiener estimation), the white Gaussian RDF is achieved \cite{zamir:1996}. In the presence of feedback loops, it is convenient to assume the
existence of a large number of identical and mutually
independent sources, equal to the lattice dimension $K$. These sources are treated
independently by the encoding/decoding scheme, except for the actual ECDQ which processes them
jointly. Thus we will only present the scheme for one source, but
the quantization noise has the properties of a high-dimensional ECDQ.
In a realistic setting, the multiple-source setting may be emulated by a single source which is divided into
$K$ long blocks and jointly encoded as $K$ parallel sources,
see~\cite{zamir:2008} for details.

\begin{figure*}[t]
\begin{center}
\psfrag{X}{$\scriptstyle X[n]$}
\psfrag{U}{$\scriptstyle U[m]$}
\psfrag{B}{$\scriptstyle B[m]$}
\psfrag{B1}{$\scriptstyle B_1[n]$}
\psfrag{B2}{$\scriptstyle B_2[n]$}
\psfrag{D1}{\hspace{-0.5mm}$\scriptstyle D_1[n]$}
\psfrag{D2}{\hspace{-0.5mm}$\scriptstyle D_2[n]$}
\psfrag{V}{$\scriptstyle V[m]$}
\psfrag{V1}{$\scriptstyle V_1[n]$}
\psfrag{V2}{$\scriptstyle V_2[n]$}
\psfrag{E}{\hspace{-1mm}$\scriptstyle E[m]$}
\psfrag{Q}{\hspace{-2mm}$\scriptstyle \mathcal{Q}(\cdot)$}
\psfrag{Y1}{\hspace{-0.5mm}$\scriptstyle Y_1[n]$}
\psfrag{Y2}{\hspace{-0.5mm}$\scriptstyle Y_2[n]$}
\psfrag{Et}{\hspace{-1mm}$\scriptstyle \tilde{E}[m]$}
\psfrag{Xh}{$\scriptstyle \hat{X}[n]$}
\psfrag{F(z)}{$\scriptstyle F(z)$}
\psfrag{A(z)}{$\scriptstyle \ A(z)$}
\psfrag{C(z)-1}{$\scriptstyle \ \ \ C(z)$}
\psfrag{LPF}{\scriptsize LPF}

\includegraphics{figs/dsq_dpcm_enc.eps}
\caption{Nested DSQ/DPCM MD encoder.}
\label{high_res}
\end{center}
\end{figure*}

\begin{figure*}[t]
\psfrag{X1}{$\scriptstyle \hat{X}_1[n]$}
\psfrag{X2}{$\scriptstyle \hat{X}_2[n-1]$}
\psfrag{Xc}{$\scriptstyle \hat{X}_C[n]$}
\psfrag{Y1}{$\scriptstyle Y_1[n]$}
\psfrag{Y2}{$\scriptstyle Y_2[n]$}
\psfrag{A(z)}{$\scriptstyle \ A(z)$}
\psfrag{LPF}{\scriptsize LPF}
\psfrag{F*(z)}{$\scriptstyle \ F^{*}(z)$}
\psfrag{G*(z)}{$\scriptstyle \ G^{*}(z)$}
\psfrag{z}{$\scriptstyle z^{-1}$}
\begin{center}
\includegraphics{figs/dsq_dpcm_dec.eps}

\caption{DSQ/DPCM MD decoder.}
\label{high_res_decoder}
\end{center}
\end{figure*}

The encoder and decoder which materialize the optimal test channel
are presented in \figref{high_res} and \figref{high_res_decoder},
respectively. All of the switches in the encoder and the decoder are
synchronized.\footnote{
It is to be understood that the switches change their positions with the upsampled rate ($m$).
Thus, in the encoder shown in Fig.~\ref{high_res}, the even samples $B_1[n]$ of $B[m]$ will go on the upper branch and the odd samples $B_2[n]$ will go on the lower branch.}
The up sampling operation followed by lowpass filtering introduces a half-sample delay on the odd samples.
This delay is corrected at the decoder by the delay operator $z^{-1}$ combined with the pair of up and downsamplers, see Fig.~\ref{high_res_decoder}.
The outputs of the quantizer blocks $\mathcal{Q}(\cdot)$ are the reconstructed values $Y_1[k]$ and $Y_2[k]$. Moreover, at each time $k$, the codeword of quantizer 1 (quantizer 2) is entropy-coded (conditioned upon the dither signal) and put into packet 1 (packet 2). The packet encoding operation is reversed at the decoder in order to obtain $Y_1[k]$ or $Y_2[k]$.

Since the two side descriptions consist of the even and odd
instances of $V[m]$, applying the predictor
$A(z)$ to each description in the original source rate is equivalent to the predictor $A(z^2)$ in the test channel.
The two quantizer
dither sequences are taken to be mutually independent, thus the two ECDQs are
equivalent to the additive noise $Z[m]$ of the test channel.
If each quantizer block is taken to be a
high-dimensional ECDQ with the required rate, and
Consequently, the two descriptions $Y_1[n]$ and $Y_2[n]$ are
equivalent to the odd and even samples, respectively, of $Y[m]$ in
the equivalent channel, and finally the whole channel from the source
to the central and side reconstructions is equivalent to the channel
from $X[n]$ to $\hat X_C[n]$ and $\hat X_S[n]$, respectively.

Since we have seen that the mutual information in the test channel achieves (with optimized spectra $\Tp$ and $\Tm$) the colored MD RDF, the encoder/decoder pair of Figs.~\ref{high_res} and \ref{high_res_decoder} is able to achieve the complete RDF for stationary Gaussian sources at all resolutions and for any desired side-to-central distortion ratio.

\begin{remark}
In the scheme shown in Fig.~\ref{high_res}, the two prediction loops are embedded within a common noise shaping loop. Alternatively, one may alter the nesting order and let the common noise shaping loop be embedded within the two prediction loops. At high-resolution conditions, there is no loss of performance by switching the nesting order. However, at general resolution, the latter approach is suboptimal.

The reason is, that for white quantization noise, the DPCM loop also shows to the outside a total white noise
(by the basic DPCM equality~\cite{jayant:1984}), while the DSQ loop
shapes the noise. Since the DPCM loop assumes white noise for
optimality \cite{zamir:2008}, it cannot be built around the shaped DSQ noise.
\end{remark}

\section{Conclusions and Discussion}\label{sec:conclusion}
A parametric formulation of the two-description symmetric RDF for
stationary colored Gaussian sources and MSE was presented.
This result was established by providing a spectral domain characterization of the optimum side and central distortion spectra. For white Gaussian sources, the optimum distortion spectral density is a two step function. For colored sources, the optimum distortion spectral density is generally not piece-wise flat but depends upon the source spectral density and the desired resolution (i.e., the desired central and side distortion levels).
It was furthermore shown that the symmetric MD RDF could be achieved by a time-domain approach based on prediction and noise-shaping.
The time domain implementation demonstrated that, at high resolutions,
it was possible to separate the mechanism responsible
for exploiting the source memory (DPCM) from the
mechanism controlling the MD coding parameters
(noise shaping).

\appendices

\section{Proof of Lemma~\ref{lem:unique_sol}}\label{app:lem:unique_sol}

There are three cases to consider, depending upon the sign of $\Xi(\ej)$.
If $\Xi(\ej)>0$, then there is one real root and two complex roots.
If $\Xi(\ej)<0$, there are three real distinct roots. Finally, if $\Xi(\ej)=0$, there is a single real triple root (if $q(\ej)=0$) or one real root and one real double root (if $q(\ej)\neq 0$)~\cite{abramovitz:1973}.
Thus, for every choice of $(\lambda_1,\lambda_2)$, one may identify
the admissible solutions
of~(\ref{eq:x1}) -- (\ref{eq:x3}), i.e., the ones that are inside the triangular region~(\ref{theta_region_spectra}) (recall Figure~\ref{fig:boundary}).

It easy to show that the discriminant \eqref{eq:Xi} satisfies:
\begin{equation} \label{eq:Xi-xi}
\Xi(\ej)=-\frac{\sxx{4}(4\lambda_1^2\sxx{2}+1)}{432\lambda_1^6}(\lambda_2-\xi^\Xi_0(\ej))(\lambda_2-\xi^\Xi_1(\ej))(\lambda_2-\xi^\Xi_2(\ej))(\lambda_2-\xi^\Xi_3(\ej)),
\end{equation}
where $\{\xi^\Xi_i(\ej)\}$ are the four real roots of $\Xi(\ej)$ given by $\xi^\Xi_0(\ej)=0, \xi^\Xi_1(\ej)=-\lambda_1$,
\begin{align}\label{eq:xi2}
\xi^\Xi_2(\ej) &= -\frac{2\sx\lambda_1 + 8\sxx{3}\lambda_1^3 - 16\sxx{2}\lambda_1^2 -3 +2\sqr{2(2\sxx{2}\lambda_1^2 + 1)^3}}{4\sx(4\sxx{2}\lambda_1^2+1)}, \\ \label{eq:xi3}
\xi^\Xi_3(\ej) &= -\frac{2\sx\lambda_1 + 8\sxx{3}\lambda_1^3 - 16\sxx{2}\lambda_1^2 -3 -2\sqr{2(2\sxx{2}\lambda_1^2 + 1)^3}}{4\sx(4\sxx{2}\lambda_1^2+1)}.
\end{align}

Since we only have to consider non-negative multipliers, $\xi_1^\Xi(\ej)\leq 0$.
Clearly, $\xi^\Xi_2(\ej)<\xi^\Xi_3(\ej), \forall \omega$. The following lemma states the signs of $\xi_2^\Xi(\ej)$ and $\xi_3^\Xi(\ej)$.

\begin{lemma}\label{lem:sign-xi3}
For $\lambda_1>0, \xi^{\Xi}_3(\ej) >0$. Moreover, the sign of $\xi_2^\Xi(\ej)$~(\ref{eq:xi2}) is given by:
\begin{equation}\label{eq:signofD}
sign(\xi^\Xi_2(\ej))= sign\left( \frac{1}{4\sx} - \lambda_1 \right).
\end{equation}
\end{lemma}
\begin{IEEEproof}

We first show that $\xi_3^\Xi(\ej)$~(\ref{eq:xi3}) is non-negative. To do so,
we show that
\begin{equation}
2\sx\lambda_1 + 8\sxx{3}\lambda_1^3 - 16\sxx{2}\lambda_1^2 -3 -2\sqr{2(2\sxx{2}\lambda_1^2 + 1)^3} \leq 0,
\end{equation}
which means that~(\ref{eq:xi3}) is non-negative.
Let $\varphi_1(\ej)= 2\sx\lambda_1 + 8\sxx{3}\lambda_1^3$ and \\ $\varphi_2(\ej) =
2\sqr{2(2\sxx{2}\lambda_1^2 + 1)^3}$ and notice that it is enough to show that $\varphi_1(\ej)<\varphi_2(\ej), \forall \omega$.
Since $\varphi_1(\ej)$ and $\varphi_2(\ej)$ are both positive functions, we may work on their squares, i.e., $\varphi_1^2(\ej)= 4\sxx{2}\lambda_1^2 + 32\sxx{4}\lambda_1^4 + 64\sxx{6}\lambda_1^6$ and
$\varphi_2^2(\ej) = 64\sxx{6}\lambda_1^6 + 192\sxx{4}\lambda_1^4 + 192\sxx{2}\lambda_1^2 + 64$. Forming the inequality $\varphi_2^2(\ej)>\varphi_1^2(\ej)$ and collecting similar terms yields $64>-188\sxx{2}\lambda_1^2 -160\sxx{4}\lambda_1^4$ which is always satisfied for $\lambda_1\in\mathbb{R}$. This proves the first part of the lemma.

We now consider the sign of $\xi_2^\Xi(\ej)$~(\ref{eq:xi2}).
Let $\varphi_1=2x + 8x^3 - 16x^2 -3$ and $\varphi_2=2\sqr{2(2x^2 + 1)^3}$.
The discriminant of $\varphi_1$ is strictly positive so $\varphi_1$ has only a single real root, which is located at $\xi=1.96973$ where we note that $\xi>\frac{1}{4}$. Moreover, $x=0\Rightarrow \varphi_1=-3$ and it follows that $\varphi_1<0$ for $x<\xi$ and $\varphi_1>0$ for $x>\xi$. Notice also that $\varphi_2>0$ for $x>0$.

At this point we let $h=\varphi_1^2 - \varphi_2^2 = -256 (z-\frac{i}2)(z+\frac{i}2)(z-\frac{1}4)^3$,
which is a fifth-order polynomial having a pair of complex conjugate roots at $x=\pm i/2$ and a real (triple) root at $x=1/4$. Thus, $h$ crosses the real line only once. Since $h=1$ for $x=0$ it follows that $h>0$ for $x<1/4$ and $h<0$ for $x>1/4$. Furthermore, $h=0$ for $x=1/4$.

Since $h>0$ for $x<1/4$ it follows that $\varphi_1^2>\varphi_2^2$ which implies that $\varphi_1+\varphi_2<0$ since $\varphi_1<0$ for $x<1/4$. The first case of~(\ref{eq:signofD}) now follows by inserting $x=\lambda_1\sx$ in $\varphi_1$ and remembering the additional sign from~(\ref{eq:xi3}). Since $h=0$ implies that $\varphi_1+\varphi_2=0$, it immediately follows that $sign(\xi^\Xi_2(\ej))=0$ for $x=\lambda_1\sx=1/4$. Finally, for $h<0$ we have $\varphi_2^2>\varphi_1^2$ which implies that $\varphi_2>\varphi_1$ since $\varphi_2$ is positive and it follows that $\xi^\Xi_2(\ej)<0$ for $\lambda_1\sx>1/4$. This proves the remaining parts of the lemma.
\end{IEEEproof}

In Fig.~\ref{fig:roots_of_Xi} we illustrate the possible sign behavior of $\Xi(\ej)$ as a function of $\lambda_2$, using Lemma~\ref{lem:sign-xi3} and the fact that
by \eqref{eq:Xi-xi}, $\lim_{\lambda_2\rightarrow \pm\infty} \Xi(\ej) = -\infty, \forall \omega$. Building on this, we prove Lemma~\ref{lem:unique_sol} by considering the following three cases.

 \begin{figure}[t]
   \begin{center}
  \psfrag{D}{$\Xi(\ej)$}
  \psfrag{l2}{$\lambda_2$}
  \psfrag{x0}{$\xi^\Xi_0$}
  \psfrag{x1}{$\xi^\Xi_1$}
  \psfrag{x2}{$\xi^\Xi_2$}
  \psfrag{x3}{$\xi^\Xi_3$}
  \psfrag{x12}{$\xi^\Xi_{1,2}$}
  \psfrag{x21}{$\xi^\Xi_{2,1}$}

 \subfigure[$\lambda_1<\frac{1}{4\sx}$]{\includegraphics[width=6cm]{figs/roots_of_D_1.eps}}\par
 \subfigure[$\lambda_1>\frac{1}{4\sx}$]{\includegraphics[width=6cm]{figs/roots_of_D_2.eps}}
 \caption{The possible zero locations for the discriminant $\Xi(\ej)$ as a function of $\lambda_2$ for a given $\omega$. We note that the $\xi_i^\Xi$'s are functions of $\omega$.}
 \label{fig:roots_of_Xi}
 \end{center}
 \end{figure}

\subsection{Negative Discriminant}\label{sec:negative_disc}
In this case
$\Xi(\ej) = q^2(\ej) + p^3(\ej) <0$ and we have three real solutions. It is easy to see that we must have $p(\ej)<0$ and $|p(\ej)|^3>q^2(\ej)$. Let $z_1(\ej)=q(\ej) + \sqr{\Xi(\ej)} = q(\ej) + i\sqr{-\Xi(\ej)}$
and $z_2(\ej)=q(\ej) - \sqr{\Xi(\ej)} = q(\ej) - i\sqr{-\Xi(\ej)}$ and notice that $s_i(\ej) = \sqrr{z_i(\ej)}, i=1,2.$ Since $\Xi(\ej), p(\ej)<0$, it follows that $|z_1(\ej)| = |z_2(\ej)| = \sqr{-p(\ej)^3} = \sqr{|p(\ej)|^3}$.
Moreover, the phase is given by
\begin{equation}\label{eq:phi1}
\phi_1(\ej) =
\begin{cases}
\arctan(\sqr{-\Xi(\ej)}/q(\ej)), & q(\ej)>0, \\
\pi + \arctan(\sqr{-\Xi(\ej)}/q(\ej)), & q(\ej)< 0, \\
\pi/2, & q(\ej) =0,
\end{cases}
\end{equation}
and $\phi_2(\ej) = -\phi_1(\ej)$. Thus, $\phi_1(\ej) \in [0;\pi]$. The value of $\phi_1(\ej)$ depends upon $q(\ej)$, which is
a third-order real polynomial in $\lambda_2$ having a negative (or zero) discriminant.
Thus, it has three real roots $\{\xi_i^q(\ej)\},i=0,1,2,$ which after some algebra can be shown to be given by
\begin{align}\label{eq:xi_q0}
\xi_0^q(\ej) &= \frac{1}{2\sx} \sqr{6}\sqr{2\sxx{2}\lambda_1^2+1}\cos\left(\frac{\phi_q(\ej)}3\right)-\frac{\lambda_1}{2}+\frac{1}{2\sx}, \\ \notag
\xi_1^q(\ej) &= -\frac{1}{4\sx} \sqr{6}\sqr{2\sxx{2}\lambda_1^2+1}\bigg(\sqr{3}\sin\left(\frac{\phi_q(\ej)}3 \right)  + \cos\left(\frac{\phi_q(\ej)}3\right) \bigg)
-\frac{\lambda_1}{2} \\  \label{eq:xi_q1}
&\quad +\frac{1}{2\sx}, \\ \notag
\xi_2^q(\ej) &= \frac{1}{4\sx} \sqr{6}\sqr{2\sxx{2}\lambda_1^2+1}\bigg(-\sqr{3}\sin\left(\frac{\phi_q(\ej)}3 \right)  + \cos\left(\frac{\phi_q(\ej)}3\right) \bigg)
-\frac{\lambda_1}{2} \\ \label{eq:xi_q2}
&\quad +\frac{1}{2\sx},
\end{align}
where
\begin{equation}\label{eq:phiq}
\phi_q(\ej) = \arctan\left(\frac{1}{9}\sqr{768\sxx{6}\lambda_1^6 + 1152\sxx{4}\lambda_1^4 + 576\sxx{2}\lambda_1^2+15}\right).
\end{equation}
Moreover, $\lim_{\lambda_2\to \infty} q(\ej) = \infty$ and $\lim_{\lambda_2\to -\infty} q(\ej) = -\infty$. Thus, $q(\ej)>0$ if $\lambda_2>\xi_0^q(\ej)$ or if $\xi_1^q(\ej)<\lambda_2<\xi_2^q(\ej)$.

With this, it is easy to show that the solutions (roots), $\{x_i(\ej)\}_{i=1}^3,$ as given by~(\ref{eq:x1})--(\ref{eq:x3}), can be written as
\begin{align}\label{eq:x1t}
x_1(\ej) &= 2\sqr{|p(\ej)|}\cos(\phi_1(\ej)/3) -\frac{a_2(\ej)}{3}, \\ \label{eq:x2t}
x_2(\ej) &= -\sqr{|p(\ej)|}\big( \cos(\phi_1(\ej)/3) +\sqr{3}\sin(\phi_1(\ej)/3) \big) -\frac{a_2(\ej)}{3}, \\ \label{eq:x3t}
x_3(\ej) &= -\sqr{|p(\ej)|}\big( \cos(\phi_1(\ej)/3) -\sqr{3}\sin(\phi_1(\ej)/3) \big) -\frac{a_2(\ej)}{3}.
\end{align}

We note that
\begin{align}
\min_{\zeta \in [0;\pi]} 2\cos(\zeta/3) \geq  \max_{\zeta \in [0;\pi]}  -(\cos(\zeta/3)-\sqr{3}\sin(\zeta/3))
\end{align}
and that
\begin{align}
\min_{\zeta \in [0;\pi]}  -(\cos(\zeta/3)-\sqr{3}\sin(\zeta/3))  \geq  \max_{\zeta \in [0;\pi]}  -(\cos(\zeta/3)+\sqr{3}\sin(\zeta/3),
\end{align}
which implies that $x_1(\ej)\geq x_3(\ej)\geq x_2(\ej)$  for any pair $(\lambda_1,\lambda_2)$ and for all $\omega\in [-\pi;\pi]$.

Let us first consider the solution $x_3(\ej)$ given by~(\ref{eq:x3t}). From Lemma~\ref{lem:x3_large} below, it follows that $x_3(\ej)>\sx/2$, which violates the spectral constraint. Moreover, since $x_1(\ej)>x_3(\ej)$, we deduce that $x_2(\ej)$
given by~(\ref{eq:x2t}) is the only admissible candidate solution.

\begin{lemma}\label{lem:x3_large}
Let $\Xi(\ej)<0$. Then, for any positive $\lambda_1$ and $\lambda_2, x_3(\ej)>\sx/2, \forall \omega$, where $x_3(\ej)$ is given by~(\ref{eq:x3t}).
\end{lemma}

\begin{IEEEproof}
Since $\Xi(\ej)<0$, we only need to consider $p(\ej)<0$. Thus, $\sqr{|p(\ej)|} = \sqr{-p(\ej)}$.
Moreover, $\frac{\partial^2}{\partial \lambda_2^2} \sqr{|p(\ej)|} = \frac{c(\ej)\,(4\sxx{2}\lambda_1^2-1)}{\sqr{|p(\ej)|^3}}$, for some everywhere positive function $c(\ej)$.
It follows that for $\lambda_1>\frac{1}{2\sx}$, $\sqr{|p(\ej)|}$ has increasing slope in $\lambda_2$. Taking the limit $\lambda_2\to\infty$ shows that the maximum slope of $\sqr{|p(\ej)|}$ is $\frac{\sx}{3\lambda_1}$.
On the other hand, $-\frac{a_2(\ej)}{3}> 3\sx/2$ and is increasing in $\lambda_2$ with constant slope $\frac{\sx}{3\lambda_1}$. Moreover,
$\cos(\phi_1(\ej)/3) -\sqr{3}\sin(\phi_1(\ej)/3) \leq 1$. Thus, $x_3(\ej) \geq -\sqr{|p(\ej)|} -\frac{a_2(\ej)}{3} > 3\sx/2$ for $\lambda_1>1/(2\sx)$.

Let us now instead assume that $\lambda_1<1/(2\sx)$. We will lower-bound $\Theta_-(\ej)$, showing that this is a high-distortion case.
To that end,  we first equate~(\ref{eq:diftheta-}) to
zero and isolate $\lambda_2$, which when inserted
into~(\ref{eq:theta1}) leads to:
\begin{equation}\label{eq:theta_plus}
\Theta_+(e^{j\omega}) =
\sx - \Tm + \frac{1}{2\lambda_1}
- \frac{\sx}{4\lambda_1\Tm},
\end{equation} where we have also substituted
$\Psi_{\lambda_1,\lambda_2}^\dagger(e^{j\omega})=\Tm$.
Using the non-negativity of $\Tm$ and $\Tp$ and
inserting~\eqref{eq:theta_plus} into $\Tm+\Tp>0$ leads to
\[\Theta_- \geq \frac{\sx}{2(2\lambda_1\sx+1)}, \] which on behalf of the assumption $\lambda_1<1/(2\sx)$ implies that $\Theta_-(e^{j\omega})>\sx/4, \forall \lambda_2$. In this case $-\sqr{|p(\ej)|} -a_2(\ej)/3$ is concave in $\lambda_2$.
Clearly, $[-\sqr{|p(\ej)|} -a_2/3]_{\lambda_2=0} >\sx/2$ and we know from above that in the limit $\lambda_2\to\infty, -\sqr{|p(\ej)|} -a_2(\ej)/3 > \sx/2$. Thus, since we can lower bound a concave function by an affine function, it follows that $-\sqr{|p(\ej)|} -a_2(\ej)/3>\sx/2$. Thus,  $x_3(\ej)>\sx/2$ as was to be proven.
\end{IEEEproof}

\subsection{Positive Discriminant}
In this case $\Xi(\ej)>0$ and we have only a single real solution given by $x_1(\ej)$ (\ref{eq:x1}).

\subsection{Zero Discriminant}\label{sec:zero_disc}
In this case $\Xi(\ej)=0$, which is possible if $-p^3(\ej)= q^2(\ej)$. The positive zeros of $\Xi(\ej)$ are given by $\xi_2^\Xi(\ej),\xi_3^\Xi(\ej)$, where the former is only positive if $\lambda_1<\frac{1}{4\sx}$.
Since $q(\ej)\neq 0$ when $\Xi(\ej)=0$, there are two real solutions i.e., $x_1(\ej)$ and $x_2(\ej)=x_3(\ej)$. We now show that $x_1(\ej)$ is the desired solution.

Let $\lambda_2=\xi_3^\Xi(\ej)$ for some $\omega$. Then it is easy to show that $s_1(\ej)=s_2(\ej)<0$ and clearly $-a_2(\ej)/3 > \sx/2$. Thus, $x_2(\ej)=x_3(\ej)=-\frac{1}{2}(s_1(\ej)+s_2(\ej))-a_2(\ej)/3>\sx/2$, which is not an admissible solution.

Now let $\lambda_2=\xi_2^\Xi(\ej)$, which is positive if and only if $0<\lambda_1<1/(4\sx)$. Moreover, assume that $q(\ej)<0$ since otherwise $x_2(\ej)=x_3(\ej)$ is clearly greater than $\sx/2$. We first show that $p(\ej)<0$ for $\lambda_2=\xi_2^\Xi(\ej)$. After some algebra, we get
\begin{equation}
p(\ej)|_{\lambda_2=\xi_2^\Xi(\ej)} =
\frac{1}{72\lambda_1^2(4\sxx{2}\lambda_1^2+1)^2}
(\varphi_1(\ej)-\varphi_2(\ej)),
\end{equation}
where
\begin{equation}
\varphi_1(\ej) = 16\sqr{2}\sqr{(2\sxx{2}\lambda_1^2+1)^3}
\sxx{2}\lambda_1^2+2\sqr{2}\sqr{(2\sxx{2}\lambda_1^2+1)^3}
\end{equation}
and
\begin{equation}
\varphi_2(\ej) = 3+26\sxx{2}\lambda_1^2+104\lambda_1^4\sxx{4}+128\sxx{6}\lambda_1^6.
\end{equation}
Since $\varphi_1(\ej)$ and $\varphi_2(\ej)$ are both positive functions, we can work on their squares, and form the inequality $\varphi_2^2(\ej) - \varphi_1^2(\ej)$, that is
\begin{align}\notag
\varphi_2^2(\ej) - \varphi_1^2(\ej) &=
1-20\sxx{2}\lambda_1^2-76\lambda_1^4\sxx{4}+1504\sxx{6}\lambda_1^6+10304\sxx{8}\lambda_1^8 \\
&\quad +22528\lambda_1^{10} \sxx{10}+16384\sxx{12}\lambda_1^{12},
\end{align}
which is clearly positive for all $0<\lambda_1\leq 1/(4\sx)$.

Let us now consider the zeros of $\Xi(\ej)=q^2(\ej) + p^3(\ej)$ and $q^2(\ej)$ from a geometric point of view.
First, $q(\ej)$ is a third-order polynomial and $q^2(\ej)$ is a non-negative sixth-order polynomial that shares zeros with $q(\ej)$.
Moreover, $p(\ej)<0$ as we established above and therefore $p^3(\ej)<0$ for $\lambda_2=\xi_2^\Xi(\ej)$.
If we consider the middle zero of $q^2(\ej)$, i.e., $\xi_2^q(\ej)$, then adding the negative function $p^3(\ej)$ to $q^2(\ej)$ will result in two zeros around the point $\lambda_2=\xi_2^q(\ej)$ instead of a single zero at $\lambda_2=\xi_2^q(\ej)$.
In fact, the smaller of the zeros becomes $\xi_1^\Xi(\ej)$ and the larger zero becomes $\xi_2^\Xi(\ej)$.
But then clearly $\xi_0^q(\ej)>\xi_2^\Xi(\ej) \geq \xi_2^q(\ej)$ and it follows that $q(\ej)<0$ for $\lambda_2=\xi_2^\Xi(\ej)$.
This shows that $x_2(\ej)=x_3(\ej)>\sx/2$ also for $\lambda_2=\xi_2^\Xi(\ej)$.

Thus, the candidate solution for $\Xi(\ej)=0$, is therefore given by $x_1(\ej)$.

\section{Proof of Proposition~\ref{prop:hr}.}\label{app:proof:prop:hr}
\subsection{Case $\lambda_1>0$ and $\lambda_2\gg 1$}\label{sec:case1}
In this case, we note that
\begin{equation}\label{eq:Xi_order}
\Xi(\ej) = -\frac{\lambda_2^4\sxx{4}}{432\lambda_1^6}(1 + 4\sxx{2}\lambda_1^2) + \mathcal{O}\big( \frac{\lambda_2^3}{\lambda_1^{3}}\big).
\end{equation}
It follows that $\Xi(\ej)<0$ for large $\lambda_2$. Furthermore,
\begin{equation}
q(\ej) = \frac{\lambda_2^3 \sxx{3}}{27\lambda_1^3} + \mathcal{O}\big( \frac{\lambda_2^{2}}{\lambda_1^2} \big)
\end{equation}
and $\phi_1(\ej) = \arctan\big(\frac{\sqr{-\Xi(\ej)}}{q(\ej)}\big).$

We use the solution $x_2(\ej)$ given by~(\ref{eq:x2}) and need to carefully address its limiting behavior in $\lambda_2$, since the dominating terms cancel.
The first-order Taylor approximation of $\arctan(x)$ is $\arctan(x)=x + \mathcal{O}(x^2), \forall |x|\leq 1$. Thus,
\begin{equation}\label{eq:arctanapprox}
\phi_1(\ej) = \arctan\bigg(\frac{\sqr{-\Xi(\ej)}}{q(\ej)}\bigg) = \frac{3\sqr{3}}{4\sx\lambda_2}\sqr{1+4\sxx{2}\lambda_1^2}  + \mathcal{O}\left(\frac{\lambda_1^3}{\sqr{\lambda_2^{3}}}\right),
\end{equation}
where the approximation becomes an equality in the limit as $\lambda_2 \to \infty$ since this implies that $\phi_1(\ej)\to 0$.
Similarly, for all $x$,
\begin{equation}\label{eq:cosapprox}
\cos(x/3)+\sqr{3}\sin(x/3) = 1 + \frac{\sqr{3}}{3}x + \mathcal{O}(x^2),\quad
(\cos(x/3)+\sqr{3}\sin(x/3))^2 = 1 + \frac{2\sqr{3}}{3}x + \mathcal{O}(x^2).
\end{equation}
Let $\alpha(\ej)=\cos(\phi_1(\ej)/3) + \sqr{3}\sin(\phi_1(\ej)/3)$. Then, we can write
\begin{align}
x_2(\ej) & = -\sqr{|p(\ej)|}\alpha(\ej) -\frac{a_2(\ej)}{3} \\ \label{eq:x2sol1}
& = \frac{|p(\ej)|\alpha(\ej)^2 - \frac{a_2(\ej)^2}{9}}{-\sqr{|p(\ej)|}\alpha(\ej) + \frac{a_2(\ej)}{3}}.
\end{align}
From~(\ref{eq:cosapprox}) and using~(\ref{eq:arctanapprox}), it follows that
\begin{equation}
\alpha(\ej)^2
= 1 + \frac{3}{2\sx\lambda_2}\sqr{1+4\sxx{2}\lambda_1^2} + \mathcal{O}\left( \frac{\lambda_1^3}{\sqr{\lambda_2^{3}}}\right).
\end{equation}
With this, we can write the numerator of~(\ref{eq:x2sol1}) as
\begin{align} \label{eq:x2nom}
|p(\ej)|\alpha(\ej)^2 - \frac{a_2(\ej)^2}{9} & =
-\frac{\sx\lambda_2}{6\lambda_1^2}\big(
2\sx\lambda_1
+1 -\sqr{1+4\sxx{2}\lambda_1^2} \big) + \mathcal{O}\big(\frac{\lambda_1}{\sqr{\lambda_2}}\big).
\end{align}
On the other hand, since $\lim_{\lambda_2\to\infty}\alpha(\ej) = 1, \forall \omega$, the denominator of~(\ref{eq:x2sol1}) can be written as (for large $\lambda_2$)
\begin{equation}\label{eq:x2den}
-\sqr{|p(\ej)|}\alpha(\ej) + \frac{a_2(\ej)}{3} \approx - \frac{2\sx\lambda_2}{3\lambda_1}.
\end{equation}
Substituting~(\ref{eq:x2nom}) and~(\ref{eq:x2den}) into~(\ref{eq:x2sol1}) yields
\begin{equation}\label{eq:dsconstant1}
\Theta_-(e^{j\omega})  = \frac{1}{4\lambda_1}\big(
2\sx\lambda_1 + 1 -\sqr{1+4\sxx{2}\lambda_1^2} \big)  + \mathcal{O}\left(\frac{\lambda_1^2}{\sqr{\lambda_2^3}}\right),
\end{equation}
so that
\begin{equation}
\lim_{\lambda_2\to\infty} \Theta_-(e^{j\omega}) = \lim_{\lambda_2\to\infty}x_2(\ej) = \frac{1}{4\lambda_1}\big(
2\sx\lambda_1 + 1 -\sqr{1+4\sxx{2}\lambda_1^2} \big).
\end{equation}

\subsection{Case $\lambda_1,\lambda_2\gg 1$}
We note that when assuming $\lambda_1/\sqrr{\lambda_2}\to 0$, then the results for finite $\lambda_1$ in Section~\ref{sec:case1} remain valid.
We rewrite~(\ref{eq:dsconstant1}) as
\begin{align}\notag
\Theta_-(e^{j\omega}) &=
\frac{1}{4\lambda_1}\big( 2\sx \lambda_1 + 1 -\sqr{1+4\sxx{2}\lambda_1^2} \big) \frac{ -(2\sx \lambda_1 + 1) -\sqr{1+4\sxx{2}\lambda_1^2}}{ -(2\sx \lambda_1 + 1) -\sqr{1+4\sxx{2}\lambda_1^2}}\\
&\quad  + \mathcal{O}\left(\frac{\lambda_1^2}{\sqr{\lambda_2^3}}\right) \\
&= \frac{1}{4\lambda_1}\,\frac{4\sx\lambda_1}{\sqr{1+4\sxx{2}\lambda_1^2} +2\sx\lambda_1+1} +  \mathcal{O}\left(\frac{\lambda_1^2}{\sqr{\lambda_2^3}}\right) \\ \label{eq:theta-order}
&=\frac{1}{4\lambda_1}c_{\lambda_1}(\ej) + \mathcal{O}\left(\frac{\lambda_1^2}{\sqr{\lambda_2^3}}\right),
\end{align}
where $c_{\lambda_1}(\ej)=\mathcal{O}(1)$ and $\lim_{\lambda_1\to\infty} c_{\lambda_1}(\ej) = 1, \forall \omega$.
Inserting this into~(\ref{eq:theta1_opt}) yields
\begin{equation}\label{eq:theta-order1}
\Theta_+(e^{j\omega}) = \frac{1}{4(\lambda_1+\lambda_2) + \mathcal{O}\left(\frac{\lambda_1^2}{\sqr{\lambda_2^3}}\right)}
+ \mathcal{O}\left(\frac{\lambda_1^2}{\lambda_1\sqrt{\lambda_2^3} + \lambda_2\sqrt{\lambda_2^3}} +\frac{1}{\lambda_1^2 + \lambda_1\lambda_2}\right).
\end{equation}
Finally, it follows from~(\ref{eq:theta-order}) and~(\ref{eq:theta-order1}) that
\begin{equation}\label{eq:hr1}
\mathop{\lim_{\lambda_1,\lambda_2\to\infty}}_{\lambda_1/\sqrr{\lambda_2}\to 0}  \lambda_1 \Theta_-  = \frac{1}{4},
\end{equation}
and
\begin{equation}\label{eq:hr2}
\mathop{\lim_{\lambda_1,\lambda_2\to\infty}}_{\lambda_1/\sqrr{\lambda_2}\to 0}  \lambda_1 \Theta_+(\lambda_1+\lambda_2)  = \frac{1}{4}.
\end{equation}

\hfill \IEEEQED

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}

\input{md_dsq_revised.bbl}

config: arxiv