\documentstyle[12pt,a4,aps]{revtex}

\textwidth=16cm
\textheight=24cm

%\font\newrm=xcmr12 scaled 1200
%\baselineskip 0.75cm plus 0.5mm minus 0.5mm

\newcommand{\R}{{\sf R\hspace*{-0.9ex}\rule{0.15ex}%
      {1.5ex}\hspace*{0.9ex}}}
\newcommand{\C}{{\sf C\hspace*{-0.9ex}\rule{0.15ex}%
       {1.3ex}\hspace*{0.9ex}}}

\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}

\begin{document}
\draft
%\newrm

\title{MAXIMUM ENTROPY ALGORITHM FOR SPECTRAL ESTIMATION PROBLEM
WITH A PRIORI INFORMATION}
\author{Roman Ugrinovsky}
\address{
Institute of Applied Physics, Russian Academy of Sciences\\
                46 Uljanov Str., 603600 Nizhny Novgorod, Russia \\
   Phone: (831-2) 38-42-81 \quad
   Fax: (831-2) 36-97-17  \\
   E-mail: roman@hydro.nnov.su }

\date{\today}

\maketitle

\begin{abstract}
The power density estimation problem with a priori information about spectrum
is considered. The problem is reduced to power density estimation problem
with spectral gaps. The solution of the problem is studied in the context of
the trigonometric moment problem with gaps.  We establish the necessary and
sufficient conditions for the existence of solution and describe all
solutions in indeterminate case.  A new definition of the entropy
functional is given.  The effective algorithm to construct maximum entropy
solution is proposed. A number of explicit illustrations are presented.
\end{abstract}

\pacs{02.30.Sa, 02.90.+p}

\section{Introduction}
The moment problem \cite{Achieser} is relevant to a variety of classical
situations of the signal processing. There are several standard methods
of dealing with this classical problem. A relative newcomer among these
involves the introduction of the concept of entropy into the function space,
selecting as preferred solution among those functions, which otherwise meet
constrains imposed by the data which yields the maximum entropy (ME) value
\cite{Agmon-Alhassid-Levin,Mead-Papanicolau}.

Following the original contributions, there has been a long debate concerning
the conceptual foundations of the ME for a series of the mathematical physics
problems \cite{Levine-Tribus} and for signal processing problems
\cite{Jaynes:1982}.
 The debate is currently  more
meaningful than ever in view of the augmenting list  of successful
practical applications \cite{Levine,Burg}, which have become possible
because of
the steadily increasing computing power available today. While our
aim  is not to engage  in further conceptual  ramifications  of the
rationale  maximum  entropy, we  shall   attempt  to  sharpen  its
mathematical foundations to extend its  applicability  to  various
concrete problems  encountered in signal processing.

   The concept of  entropy  has been known for a long time in information
theory \cite{Shannon,Jaynes}, but its introduction into signal processing
 was due to Burg \cite{Burg}. In 1967 he gave the new conception of spectral
estimation and proposed the method by which one directly estimates the
relative parameters $\gamma_i~,i=0,~1, \dots,~N$  from measurement data
$c_n~,n=0,~1, \dots,~N$,
and gets the power spectral density estimation (PSDE) by means of the equation
\begin{equation}
P(\omega)={c_0\over{| \sum\nolimits_{i=0}^N \gamma_i e^{-j \omega i}|^2}} .
\label{Somega}
\end{equation}

The theoretical advantage of the  ME  spectral analysis and
the high-resolution property of Burg's algorithm attracted great attention.
Many efforts have been made to explain the properties of ME estimation
\cite{Jaynes:1982}
and also to employ it in geophysics, speech and image reconstruction, radar,
sonar and other areas \cite{Burg,Kay-Marple}.
Furthermore, various modifications of Burg's algorithm and some better spectral
estimation methods have continuously been proposed
\cite{Wernecke,Notic-Knafel-Turchin-Ugrinovsky-Ugrinovsky}.

In this paper we consider the PSDE problem in the most general statement and
propose the novel method of its solution based on moment method technique.
The statement of problem assumes using additional a priori information about
spectrum. For example, the separate regions of spectrum or the original
spectrum support may be known.
In view of such information the higher-resolution method of spectral
estimation is proposed.

We establish necessary and sufficient conditions of the
solvability of the problem and for existence of the ME solution, and give the
constructive way to compute this ME solution.  To achieve the aim we describe
all solutions of the trigonometric moment problem with gaps corresponding to
the considered PSDE problem and introduce the new definition of the entropy
functional.  Next, we show that this definition coincides with the generally
used one in classical cases. The new definition of entropy functional allows
us to find the construction of the ME solution for nontrivial spectrum
supports in terms of polynomials orthogonal with respect to given moment
sequence. To compute the orthogonal polynomials corresponding to the known
data (i.e.,  moments) we suggest using the modified Chebyshev algorithm
\cite{Korzh-Ovcharenko-Ugrinovsky} or Levinson algorithm \cite{Levinson}.
A number of explicit illustrations are presented.


\section{Power density spectral estimation problem with
a~priori information}

Consider the problem of estimating a power spectral density function given
only that it is positive on the spectral support, zero outside and compatible
with the known values of autocorrelation sequences.
In addition a priori information about spectrum is assumed to be given,
i.e., we assume the spectrum is unknown in full interval of
frequencies, but one or several regions of spectrum are known, and consider
the PSDE problem with such a priori information.  The aim of spectral
estimation in this case is to find the function which agrees with the set of
the known values of autocorrelation function $c_k$, $k=1,~2,\ldots,~N$ and
given a priori information.  To be rigorous we assume the original spectrum
$S(\theta)$ is known on intervals $[r_1,t_1]$, $[r_2,t_2]$, \ldots,
$[r_m,t_m]$, where $-\pi\le r_1<t_1<r_2<\ldots<t_m\le\pi$, i.e., for all
$\theta\in \cup_{i=1}^m [r_i, t_i]=I_m$.
In view of the given a priori information, the original PSDE problem
reduces to PSDE problem with a priori information about spectrum,
that is to determine the nonnegative function
$P(\theta)$ such that
\begin{equation}
P(\theta)\geq 0,\quad\hbox{\rm for}\quad \theta\in [-\pi,\pi],
\label{Pge0}
\end{equation}
which is compatible with the known values $c_k$ of the
autocorrelation function:
\begin{equation}
{1 \over {2 \pi}} \int\limits_{- \pi}^\pi
 P(\theta) e^{-ik\theta}\,d\theta = c_k ,
\quad k=0,~1,\ldots,~N .
\label{c_k}
\end{equation}
and with the given a priori information about spectrum
\begin{equation} P(\theta)=S(\theta),\quad\hbox{\rm for}\quad\theta\in I_m .
\label{P=S}
\end{equation}

Let us note that the statement of problem (\ref{Pge0})--(\ref{P=S})
is the most general one of the PSDE problem. Really, this statement
allows us to consider PSDE problem with a priori information about
spectrum in its different modification as a uniform problem.
Moreover, this statement contains both a traditional PSDE problem
that corresponds to case $I_m=\emptyset$, and the PSDE problem with gaps,
which is to determine the nonnegative function $P(\theta)$ with the gaps
$I_m$, where $I_m$ is the given system of intervals, i.e., the function
$P(\theta)$, such that
$$P(\theta)\ge0, \quad P(\theta)=0\quad\hbox{\rm
for}\quad\theta\in I_m,$$
and which is compatible with the known values $c_k$
of the authocorrelation function:
$$ {1 \over {2 \pi}} \int\limits_{-
\pi}^\pi P(\theta) e^{-ik\theta}\,d\theta = c_k , \quad k=0,~1,\ldots,~N .
$$

The solution of PSDE problem (\ref{Pge0})--(\ref{P=S})
% This
means that we need spectrum only for frequencies outside
$I_m$ to be estimated.

One possible approach to solve the problem (\ref{Pge0})--(\ref{P=S})
is to refuse %reject
the given
a priori information and estimate spectrum using one of known methods
of spectral estimation, for example, maximum entropy method,
maximum likelihood method or other. This approach possesses a fundamental
disadvantage. First of all, the spectral estimation determined by means of
a known method may not be in agreement with the given a priori
information.  Besides, when the spectral estimation with good resolution is
necessary it seems unreasonable to refuse
the additional information about spectrum, because the given a priori
information can be used to improve the resolution.

An other possible approach is to take into account the given a priori
information about spectrum and consider the PSDE problem with a priori
information.
In this case the problem arises of constructing a high-resolution method
of spectral estimation which will use the given additional a priori
information.
Bellow we will consider this approach.

First of all let us show that PSDE problem
with the given a priori information can be reduced to PSDE problem with
gaps in spectrum.

Let $P(\theta)$ be any solution of the PSDE problem with a priori information
(\ref{Pge0})--(\ref{P=S}). For convenience we define the function
corresponding to a priori information
\begin{equation}
P_a(\theta)=\left\{ \begin{array}{ll}
                            S(\theta) & \mbox{ if $\theta\in I_m$} \\
                            0         & \mbox{otherwise}
                     \end{array}
             \right.
\label{P_a}
\end{equation}
and function with gaps $I_m$
\begin{equation}
P_g(\theta)=\left\{ \begin{array}{ll}
                            0         & \mbox{ if $\theta\in I_m$} \\
                            P(\theta) & \mbox{otherwise}
                     \end{array}
             \right.
\label{P_g}
\end{equation}
Then any solution $P(\theta)$ of the PSDE problem (\ref{Pge0})--(\ref{P=S})
can be represented by means of the sum of these functions, i.e.,
the representation
$$
P(\theta)=P_a(\theta)+P_g(\theta)
$$
holds for all $\theta\in[-\pi,\pi]$.

Therefore, for all $k=0,~1,~\ldots,~N$ we have
\begin{equation}
{1\over{2\pi}}\int\limits_{-\pi}^\pi P(\theta)e^{-ik\theta}d\theta=
 {1\over{2\pi}}\int\limits_{-\pi}^\pi P_a(\theta)e^{-ik\theta}d\theta+
 {1\over{2\pi}}\int\limits_{-\pi}^\pi P_g(\theta)e^{-ik\theta}d\theta
\label{V}
\end{equation}
Due to relation (\ref{c_k})
the integral on the left-side of the relation (\ref{V}) is equal to
$c_k$, $k=0,~1,~\ldots,~N$. Next, the first integral on the right-side
of the relation (\ref{V}) can be calculated for all integer $k$ because of
$$
{1\over{2\pi}}\int\limits_{-\pi}^\pi P_a(\theta)e^{-ik\theta}d\theta=
{1\over{2\pi}}\int\limits_{I_m} S(\theta)e^{-ik\theta}d\theta
$$
and $S(\theta)$ is known on the  set $I_m$.

Let $c_k^a$, $k=0,~1,~\ldots,~N$ be the correlation coefficients
that correspond to function $P_a(\theta)$, i.e., let
$$c_k^a=\int\limits_{I_m} S(\theta) e^{-ik\theta}d\theta$$
Then from (\ref{V}) we have
\begin{equation} {1\over{2\pi}}\int\limits_{-\pi}^\pi
P_g(\theta)e^{-ik\theta}d\theta= c_k^g, \quad k=0,~1,~\ldots,~N,
\label{rectangle}
\end{equation}
where
$$c_k^g=c_k-c_k^a$$
is known values.

Let us note that $P_g(\theta)$ is a nonnegative function with the gaps $I_m$.
Hence due to relations (\ref{rectangle}) this function can be considered as
a power spectral estimation corresponding to the correlation
coefficients $c_k^g$, $k=0,~1,~\ldots,~N$.
Hence, the most general statement of the PSDE problem is reduced to
PSDE problem with gaps.

So, to solve the PSDE problem with a priori information
(\ref{Pge0})--(\ref{P=S}) it is sufficient to find the nonnegative function
$P_g(\theta)$ with the gaps $I_m$, which  satisfies the conditions
(\ref{rectangle}), i.e., it is sufficient to solve the PSDE problem with gaps
in spectrum.

In next sections  we will consider the PSDE problem with some spectral
gaps in context of trigonometric moment problem with gaps.


\section{The trigonometric moment problem with gaps}
The trigonometric moment problem with gaps  consists of a collection of complex
numbers
\begin{equation}
c_k={1\over{2\pi}}\int\limits_{-\pi}^\pi P(\theta)e^{ik\theta}\,d\theta ,
\qquad k=0,~1,\ldots,~N,
\label{1.c_k}
\end{equation}
where $P(\theta)$ is an unknown positive function of $\theta$ with
the gaps $I_m$. Assuming that the moments $c_k$, $k=0,~1,\ldots,~N$ of
$P(\theta)$ and the gaps $I_m$ are exactly known, the problem is to determine
the function $P(\theta)$. Not any set of numbers $c_k$, $k=0,~1,\ldots,~N$
can be the moments of nonnegative function with the gaps $I_m$, i.e., the
trigonometric moment problem with the gaps $I_m$ for arbitrary sequence of
the moments $c_k$, $k=0,~1,\ldots,~N$ can be unsolvable. Below we will give
the effective criterion of the solvability of the considered moment problem
with the given gaps $I_m$.

Let $\{ {\cal P}_C\}$ be a space of polynomials  with complex coefficients
on the unit circle~$\C=\{z\colon |z|=1\}$,
i.e., a space of trigonometric polynomials, and let~$< , >$ denote a
scalar product on the space~$\{ {\cal P}_C\}$ with the property
\begin{equation}
  <P,Q>=<P{\overline Q},1>
\qquad \hbox{\rm for the polynomials}\quad P,Q\in\{ {\cal P}_C\}
\end{equation}

We shall assume that the Toeplitz determinants are nonzero:
   $$\Delta_n=\det\left| \, \matrix{
   {<1,1>}&{<z,1>}&\ldots&{<z^n,1>}\cr
   \overline {<z,1>}&{<1,1>}&\ldots&{<z^{n-1},1>}\cr
   \vdots&\vdots&\ddots&\vdots\cr
   \overline {<z^n,1,>}&\overline {<z^{n-1},1>}&\ldots&{<1,1>}\cr
   }\right|\neq0 ,\qquad n=0,~1,\ldots .$$

It is well known (see for example \cite{Achieser}) that there exists a
$< , >$-orthogonal system of polynomials~$\{P_n(z)\}$ such that

1) $\hbox{\rm deg} P_n(z)=n$

2) the highest-degree coefficient is positive

3) $<P_n(z),P_m(z)>=0,$ if~$ n\neq m.$

We will consider the monic system of polynomials $P_k(z)$, which are
characterized by the highest degree coefficient being equal to $1$.
The monic polynomials~$\{P_n(z)\}$ satisfy the recurrence relation
for the orthogonal sets of polynomials:
\begin{equation}
P_{n+1}(z)=zP_n(z)+b_{n+1}P_n^*(z),
\label{rec}
\end{equation}
where~$R^{*}(z) \buildrel \rm def \over = z^k \overline R(1/z)$ for any
polynomial~$R(z)$ of degree~$k$; $b_k$ are
some constants for which formal analytic relations in determinant
form can be obtained. Note that $b_0$ does not take part in the
recursion~(\ref{rec}), however, below it will be convenient to define
  $$|b_0|^2=1-<1,1>.$$
>From recurrence relations~(\ref{rec}) the following three-term recurrence
relation immediately follows:
\begin{equation}
  b_{n+1}P_{n+2}(z)=(zb_{n+1}+b_{n+2})P_{n+1}(z)-
                      b_{n+2}(1-|b_{n+1}|^2)zP_n(z).
\label{ttrec}
\end{equation}

The following criterion of solvability of the moment problem (\ref{1.c_k})
with the gaps $I_m$ holds \cite{Ugrinovsky}:
\begin{theorem}
The trigonometric moment problem (\ref{1.c_k}) with the  gaps
$I_m=\{z:\ z=e^{i\theta},\ \theta\in\cup_{i=1}^m [r_i, t_i],
\ -\pi\le r_1<t_1<r_2<\ldots<t_m\le\pi\}$ is solvable if and only if
the inequalities
$$|a_k|<1, \quad k=0,~1,\ldots,~N,
\qquad |b_k^{(l)}|<1, \quad k=0,~1,\ldots,~N-1,\quad l=1,\ldots,~m$$
hold, where $\{ a_k\}_{k=0}^N$ and $\{ b_k^{(l)}\}_{k=0}^{N-1}$ are
the  coefficients of the recurrence relation of polynomials orthogonal
with respect to the moment sequences $\{ c_k\}_{k=0}^N$ and
$\{ \mu_k^{(l)} \}_{k=0}^{N-1}$ $l=1,\ldots,~m$, respectively,
$$\mu_k^{(l)}=2\cos{{t_l-r_l}\over 2}-e^{i{{r_l+t_l}\over 2}}c_{k+1}
- e^{-i{{r_l+t_l}\over 2}}c_{k-1},\qquad l=1,~2,\ldots,~m\ .$$
\end{theorem}

Let us now assume that the necessary conditions of solvability of the
trigonometric moment problem with the gaps $I_m$ are true. Then all
functions $P(\theta)$, for which the relation (\ref{1.c_k}) is fulfilled
for all $k=0,~1,\ldots, N$, can be described by means of the linear-fractional
transformation.
\begin{theorem}
\label{XXX}
Let $P_N(z)$ and $Q_N(z)$ be the monic polynomials of first and second type
orthogonal with respect to given moment sequence $\{ c_k\}_{k=0}^N$, and
let $I_m$ be the given gaps. Let also the moment problem (\ref{1.c_k})
with the gaps $I_m$ be indeterminate. Then there exist the linear-fractional
functions
$$\Phi_k(z,g)={
               {x_k(z) g+y_k(z)}
              \over
               {u_k(z) g+v_k(z)}
              } , \quad k=0,~1,\ldots,~m-1 ,$$
where $x_k(z)$, $y_k(z)$, $u_k(z)$ and  $v_k(z)$ are the polynomials,
such that all solutions of the moment problem (\ref{1.c_k}) with gaps
can be described by means of the linear-fractional transformation
\begin{equation}
P_\varepsilon(z)={{c_0}\over{h_N}}\cdot
 \Re e           {
                  {zQ_N(z)\varepsilon(z)+Q^*_N(z)}
                 \over
                  {zP_N(z)\varepsilon(z)-P^*_N(z)}
                 }, \quad z=e^{i\theta}
\label{lft}
\end{equation}
of the function $\varepsilon(z)$, such that
$$\Phi_k(z,\varepsilon(z))\in B, \quad k=0,~1,\ldots,~m-1 .$$
$B$ is a class of functions holomorphic in $\C$ and bounded there
in modulus by 1; $P^*_n(z)=z^n \overline P_n(1/z)$,
$Q^*_n(z)=z^n \overline Q_n(1/z)$,
$h_n={{\Delta_n}\over{\Delta_{n-1}}}$.
\end{theorem}

\section{Maximum entropy solution}
One of the more interesting problem in the trigonometric moment problem
with gaps is to find the extreme (i.e., ME) solution. Obviously,
the consideration of the extreme problem is possible only for indeterminate
moment problem.

Let $\Sigma_{I_m, P(\theta)}$ be the set of solutions of the indeterminate
moment problem (\ref{1.c_k}) with the gaps $I_m$.
In correspondence with Theorem~\ref{XXX} the set
$\Sigma_{I_m, P(\theta)}$ is described by means of
a set $B_{I_m}$ of the holomorphic function $\varepsilon(z)$.

Let us consider the Burg's functional \cite{Arov-Krein,Dym-Gohberg}
\begin{equation}
h[P,\zeta]={1\over{4\pi}} \int\limits_{-\pi}^\pi
{{(1-|\zeta |^2)}\over{|e^{i\theta}-\zeta |^{2}}}
\ln P(\theta ) d\theta ,  \qquad (|\zeta|<1) .
\label{Burgf}
\end{equation}
It is obvious that this functional is undetermined for all solutions
of the moment problem (\ref{1.c_k}) with the gaps $I_m\neq \{\emptyset\}$,
i.e.,
$$h[P,\zeta]=+\infty\qquad P\in\Sigma_{I_m, P(\theta)} ,$$
and so the extreme problem corresponding to Burg's functional to find ME
solution of the moment problem with gaps can not be considered.

But, due to Theorem~\ref{XXX} any solution $P(\theta)$ of the moment problem
with gaps $I_m$ can be expressed by means of linear-fractional
transformation~(\ref{lft}) with the holomorphic function $\varepsilon(z)$
corresponding to this solution.
Using properties of orthogonal polynomials \cite{Achieser} it is easy to show
that
$$
P_\varepsilon(z)=
{
{1-|\varepsilon(z)|^2}
\over
{|z {P_(z)}\varepsilon(z)-{P_n^*(z)}|^2 }} , \quad z=e^{i\theta} .
$$

Then the Burg's functional $h[P,\zeta]$ can be written in the form
$$h[P,\zeta]={1\over{4\pi}} \int\limits_{-\pi}^\pi
{{(1-|\zeta |^2)}\over{|e^{i\theta}-\zeta |^{2}}}
\ln (1-|\varepsilon(e^{i\theta})|^2) d\theta +
{1\over{2\pi}} \int\limits_{-\pi}^\pi
{{(1-|\zeta |^2)}\over{|e^{i\theta}-\zeta |^{2}}}
\ln{|z {P_(z)}\varepsilon(z)-{P_n^*(z)}|} d\theta
$$

Let us now introduce into consideration the functional
\begin{equation}
\widehat h[\varepsilon,\zeta]={1\over{4\pi}} \int\limits_{-\pi}^\pi
{{(1-|\zeta |^2)}\over{|e^{i\theta}-\zeta |^{2}}}
\ln (1-|\varepsilon(e^{i\theta})|^2) d\theta ,
\qquad (|\zeta|<1) ,
\label{ruf}
\end{equation}
on the set $B_{I_m}$.
Obviously,
this functional %(\ref{ruf})
is defined for all function
$\varepsilon(z)\in B_{I_m}$.

The following proposition holds:
\begin{theorem}
\label{YYY}
Let (\ref{1.c_k}) be indeterminate moment problem with the gaps
$I_m=\{ z:\ z=e^{i\theta},\ \theta\in\cup_{i=1}^m [r_i, t_i],
\ -\pi\le r_1<t_1<r_2<\ldots<t_m\le\pi \} $; $\Sigma_{I_m, P(\theta)}$
be the set of solutions of the moment problem (\ref{1.c_k}); $B_{I_m}$
be a set of the holomorphic functions corresponding to the set
$\Sigma_{I_m, P(\theta)}$ due to relation (\ref{lft}) and let
\begin{equation}
\varepsilon_0 (z)=\arg\inf_{\varepsilon(z)\in B_{I_m}} |\varepsilon(z)|\ .
\end{equation}
Then the inequality
\begin{equation}
\widehat h[\varepsilon,\zeta] \leq  \widehat h[\varepsilon_0,\zeta]\ .
\end{equation}
holds for all functions $\varepsilon(z)\in B_{I_m}$.
\end{theorem}

Let us note, that for the moment problem without gaps
(i.e., when $I_m=\{\emptyset\}$) we have
$$\varepsilon_0(z)\equiv 0$$
and, hence,
$$P_{\varepsilon_0}(z)=-{{c_0}\over{h_N}}\Re e{{Q_N^*(z)}\over{P_N^*(z)}}
={1\over{|P_N(z)|^2}} , $$
i.e., we get the well-known ME solution of the classical trigonometric moment
problem without gaps. Moreover, in this case the following proposition holds:
\begin{theorem}
Let $I_m=\{\emptyset\}$. Then for all $\varepsilon(z)\in B$ the relation
$$h[P_\varepsilon,\zeta]=\widehat h[\varepsilon,\zeta]+C$$
is true, where the function $P_\varepsilon(z)$ is determined by (\ref{lft})
and $C$ is a constant.
\end{theorem}

Thus, in spite of the fact that Burg's functional is undetermined for the solutions of
the
trigonometric moment problem with gaps the functional (\ref{ruf}) is determined
and in case of classical moment problem (i.e., $I_m\equiv \{\emptyset\}$)
these functions have the same extreme solution taking into account of course
the linear-fractional transformation (\ref{lft}).


\section{Numerical procedure}

To use the ME algorithm described in the previous section the effective
numerical
procedure is necessary. The main problem is to construct the polynomials
$P_N(z)$ and $Q_N(z)$ and the holomorphic function $\varepsilon_0(z)$.
But as follows from Theorem~\ref{YYY} the function $\varepsilon_0(z)$
can be expressed through the polynomials $\{ P_n(z) \}$ and $\{ Q_n(z) \}$
and the gaps $I_m$. Thus, the main problem is to propose the effective technique to get
the set of polynomials $\{ P_n(z) \}$, orthogonal with respect to given moment
sequence, and polynomials $\{ Q_n(z) \}$.
Below the numerical procedure to construct
the set of monic polynomial $\{ P_n(z) \}$ is given.

Let $\{P_k(z)\}$ be the set of monic polynomials orthogonal with respect
to given moment sequence $\{ c_k\}_{k=0}^N$ with the recurrence relations
(\ref{ttrec}) and let $\{\pi_k(z)\}$ be another set of
monic orthogonal polynomials with the following recurrence relations
\begin{equation}
\beta_{n+1}\pi_{n+2}(z)=(z\beta_{n+1}+\beta_{n+2})\pi_{n+1}(z)-
                  \beta_{n+2}(1-|\beta_{n+1}|^2)z\pi_n(z) .
\label{ttrecpi}
\end{equation}
The coefficients~$\{\beta_n\}$ are  assumed to be known and
the coefficients~$\{ b_n\}$ in recurrence relation (\ref{ttrec}) are assumed
to be unknown.

\begin{definition}
The matrix~$T_{P,\pi}$ of modified moments,
where $$\|T_{P,\pi}\|_{k,l}=\sigma_{k,l}{\buildrel \rm def \over =}
<P_l,\pi_k>$$
will be called a {\it Chebyshev matrix}.
\end{definition}

Note that $\sigma_{k,l}=0$ for all~$k>l$, i.e. the Chebyshev matrix is
triangular.

The entries of the Chebyshev matrix~$T_{P,\pi}$ and the coefficients of
the three-term recurrence formulas (\ref{ttrec}) and~(\ref{ttrecpi})
are connected by
relations that can be considered as a recursion process to determine
entries of the Chebyshev matrix and unknown elements~$b_k$. Writing
it down in a recursive manner, we have:

{{\it the initial conditions:}~$\sigma_{-1,l}=0,\quad\sigma_{0,l}=<1,\pi_l>,
\quad l=0,~1,\ldots,N ,$ and $b_0=1$

{\it and for all} $k=1,~2,\ldots,~N$
  $$b_k={1\over{b_{k-1}}}\bigl[-{\sigma_{k-1,k}\over\sigma_{k-1,k-1}}
           +{\overline \beta_k\over\overline \beta_{k-1}}
           +(1-|b_{k-1}|^2){\sigma_{k-2,k-1}\over\sigma_{k-1,k-1}}
           -(1-|b_{k-1}|^2){\overline \beta_k\over\overline \beta_{k-1}}
            (1-|\beta_{k-1}|^2){\sigma_{k-2,k-2}\over\sigma_{k-1,k-1}}
             \big]$$
\rightline {{~}~$(m)-T_C$}
\[\sigma_{k,l}=\left\{
\matrix{
        0, & \qquad l=0,~1,\ldots,~{k-1};  \cr
       {\overline \beta_l\over\overline \beta_{l-1}}\sigma_{k,l-1}+
        {{b_k}\over {b_{k-1}}}\sigma_{k-1,l}-
        {{b_k}\over {b_{k-1}}}{\overline \beta_l\over\overline \beta_{l-1}}
                                 \sigma_{k,l-1}+\sigma_{k-1,l-1}& \cr
      -{\overline \beta_l\over\overline \beta_{l-1}}(1-|\beta_{l-1}|^2)\sigma_{k-1,l-2}-
        {{b_k}\over{b_{k-1}}}(1-|b_{k-1}|^2)\sigma_{k-2,l-1} \cr
      +{{b_k}\over{b_{k-1}}}(1-|b_{k-1}|^2)
        {\overline \beta_l\over\overline \beta_{l-1}}(1-|\beta_{l-1}|^2)\sigma_{k-2,l-2},
        & \qquad l=k,~k+1,\ldots,~N \cr
}
\right. \]
}

In the case~$\pi_k(z)=z^k$, i.e. $\beta_k=0$ for~$k=0,~1,\ldots,$ the
$(m)-T_C$-recurrence relations are transformed into the well-known Levinson
algorithm \cite{Levinson}:

{\it step~1:} $P_0(z)=1,\qquad\sigma_{0,0}=<1,1>,$

{\it step~2:} for~$k=1,~2,\ldots,~N$
\begin{eqnarray*}
  b_k     &=&-{1\over\sigma_{k-1,k-1}}<z P_{k-1}(z),1>,\cr
  \sigma_{k,k}&=&(1-|b_k|^2)\sigma_{k-1,k-1},\cr
  P_k(z)    &=&z P_{k-1}(z)+b_k P_{k-1}^*(z).
\end{eqnarray*}

Let us note now that if the orthogonal polynomial $P_n(z)$ of the first type
is known the orthogonal polynomial $Q_n(z)$ of the second type can be found
by means of the relation
\begin{equation}
\left[ \, \matrix{
          q_{n,0}  \cr
          q_{n,1}  \cr
          q_{n,2}  \cr
          \vdots   \cr
          q_{n,n}  \cr
}\right] =-{1\over{c_0}} \left[ \, \matrix{
            0    &   c_1  &   c_2  &   c_3  & \ldots &   c_n    \cr
            0    &   c_0  &  2c_1  &  2c_2  & \ldots & 2c_{n-1} \cr
            0    &    0   &   c_0  &  2c_1  & \ldots & 2c_{n-2} \cr
          \vdots & \vdots & \vdots & \vdots & \ddots & \vdots   \cr
            0    &    0   &    0   &    0   & \ldots &   c_0    \cr
}\right]
\left[ \, \matrix{
          p_{n,0}  \cr
          p_{n,1}  \cr
          p_{n,2}  \cr
          \vdots   \cr
          p_{n,n}  \cr
}\right] \ ,
\label{P_k-to-Q_k}
\end{equation}
where $[q_{n,0},~q_{n,1},\ldots,~q_{n,n}]^T$ is the vector of coefficients
of the polynomial  $Q_n(z)=\sum\limits_{m=0}^n q_{n,m} z^m $,
$[p_{n,0},~p_{n,1},\ldots,~p_{n,n}]^T$ is the vector of coefficients
of the polynomial  $P_n(z)=\sum\limits_{m=0}^n p_{n,m} z^m $ and
$c_0,~c_1,\ldots,~c_N$ are the given moments.

\section{Example}
In this section we give an illustrative example and consider  the PSDE
problem with one gap $[r,t]$, i.e., further we will suppose $m=1$. In this
case the above described numerical algorithm is transformed as follows:
\begin{description}
\item[Step 1] We find the system of monic polynomials $\{P_k(z)\}_{k=0}^N$
orthogonal with respect to given moment sequence by means of the
$(m)-T_C$-recursion or Levinson algorithm.
\item[Step 2] From relation (\ref{P_k-to-Q_k}) we find the polynomials
$Q_{N-1}(z)$ and $Q_N(z)$.
\item[Step 3]  We find the functions
\begin{eqnarray*}
{\cal A}(z,\zeta)&=&Q_N(z)P_N(z,\zeta)-P_N(z)Q_N(z,\zeta)     , \\
{\cal B}(z,\zeta)&=&P_N^*(z)Q_N(z,\zeta)+Q_N^*(z)P_N(z,\zeta) ,
\end{eqnarray*}
where
\begin{eqnarray*}
P_N(z,\zeta)&=&P_N(z)P_{N-1}(\zeta)-P_N(\zeta)P_{N-1}(z) ,     \\
Q_N(z,\zeta)&=&Q_N(z)P_{N-1}(\zeta)-Q_{N-1}(z)P_N(\zeta)
\end{eqnarray*}
for $\zeta=\zeta_r$ and $\zeta=\zeta_t$ ($\zeta_r=e^{-i r}$, $\zeta_t=e^{-i t}$)
and the functions
\begin{eqnarray*}
\alpha(z)&=&|z|^2 [|{z-\zeta_t}|^2 |{\cal A}(z,\zeta_r)|^2-
                   |{z-\zeta_r}|^2|{\cal A}(z,\zeta_t)|^2] ,   \\
\beta(z)&=&\overline z [|{z-\zeta_t}|^2
\overline{{\cal A}(z,\zeta_r)}{{\cal B}(z,\zeta_r)}
                        -|z-\zeta_r|^2
 \overline{{\cal A}(z,\zeta_t)} {{\cal B}(z,\zeta_t)}] ,   \\
\gamma(z)&=&|{z-\zeta_t}|^2 |{\cal B}(z,\zeta_r)|^2
          -|z-\zeta_r|^2 |{\cal B}(z,\zeta_t)|^2 .
\end{eqnarray*}
\item[Step 4] We calculate the function
$$\varepsilon_0(z)=\left({{|\beta(z)|}\over{\alpha(z)}}-R(z)\right)
e^{i {{\Im m\beta(z)}\over{\Re e\beta(z)}}}  , $$
where
$$O(z)={{\beta(z)}\over{\alpha(z)}}, \qquad
R(z)={{\sqrt{|\beta(z)|^2-\alpha(z)\gamma(z)}}\over{\alpha(z)}} ,$$
for all $z=e^{i\theta}$, $\theta\in [-\pi,\pi]$.
\item[Step 5] We find the PSDE by
$$
P_{\varepsilon_0}(z)={{c_0}\over{h_N}}\cdot
 \Re e           {
                  {zQ_N(z)\varepsilon_0(z)+Q^*_N(z)}
                 \over
                  {zP_N(z)\varepsilon_0(z)-P^*_N(z)}
                 }, \quad z=e^{i\theta}, \quad \theta\in [-\pi,\pi] .
$$
\end{description}

Let us note that the functions $P_N(z,\zeta)$, $Q_N(z,\zeta)$,
${\cal A}(z,\zeta)$, ${\cal B}(z,\zeta)$, $\alpha(z)$, $\beta(z)$ and
$\gamma(z)$ from Steps~1--3 are  polynomials and, therefore, it is sufficient
to perform all calculations only with the coefficients of these polynomials,
but not with their values. The computation of values of the functions
$\alpha(z)$, $\beta(z)$ and $\gamma(z)$ for $z=e^{i\theta}$  is necessary
only in Step~4 to calculate the function $\varepsilon_0(z)$ for $z=e^{i\theta}$
and then the PSDE $P_{\varepsilon_0}(z)$ corresponding to the function
$\varepsilon_0(z)$.

Now we will consider the elementary numerical example for correlation sequence
$$c_k=\sum_{l=1}^L q_l {{\alpha_l}\over{4\alpha_l^2-k^2}}
            {{\sin{{\pi k}\over{2\alpha_l}}}\over{{\pi k}\over{2\alpha_l}}}
e^{-ika_l}+\epsilon\delta_k ,$$
where $q_l, a_l$ and $\alpha_l$, $l= 1,~2,~\ldots,~L$ and $L$ are given values,
$\epsilon$ is a parameter of regularization. Let us note, that the original
spectral power density has the gaps
$[a_l+{\pi\over{2\alpha_l}},a_{l+1}-{\pi\over{2\alpha_l}} ]$, $l=0,\ldots,~L-1$
with the widths
$\Delta\theta=a_{l+1}-a_l-{\pi\over 2}
               \left( {1\over{\alpha_l}}+{1\over{\alpha_{l+1}}}\right)$.
The original power spectral density for $L=2$, $a_1=-\pi /7, a_2=\pi /9,
q_1=10, q_2=100, \alpha_1=16.1, \alpha_2=2.37$ is shown in Fig.~1%\ref{fig1}.

We  assume that finite number of moments is known, for example $N=21$.
The PSDE obtained by means of proposed technique for the pointed out parameters
and different values of the regularization parameter $\epsilon$ is presented
in Fig.~2.%\ref{fig2}.
We used the algorithm for $m=1$ and $r=a_1+{\pi\over{2\alpha_1}}+0.15$,
$t=a_2-{\pi\over{2\alpha_2}}-0.15$.
As one can see from the figures the proposed algorithm
helps to get the PSDE with good resolution and to keep gaps in the PSDE.

\section{Conclusion}
We considered the spectral power density estimation problem with gaps.
The maximum entropy (ME) approach developed in the paper to solve the problem
is universal and can be applied to the solution of other problems.
The specificity of PSDE problem prompts us to consider trigonometric
moment problem and construct the ME algorithm to find the extreme solution
of trigonometric moment problem with gaps. In complete analogy the classic
moment problem can be considered and the ME algorithm to find the extreme
solution of moment problem with gaps can be obtained.

\section*{Acknowledgements}
I am grateful to I.E.~Ov\^carenko for his direct as well as indirect
influence and very useful discussions.

This work was supported by the Russian Foundation of Fundamental Researches,
Grant No. 93-02-16043.


%\begin{figure}[s]
%\vspace{14cm}
%\caption{\label{fig1}}
%\end{figure}

%\begin{figure}[s]
%\vspace{14cm}
%\caption{\label{fig2}}
%\end{figure}

\begin{thebibliography}{20}
\bibitem{Achieser} {N.I.~Achieser,
{\it The Classical Moment Problem and Some related Questions in Analysis}
(Oliver \& Boyd, Edinburgh, 1965)}
\bibitem {Agmon-Alhassid-Levin} N.~Agmon, Y.~Alhassid, R.D.~Levin,
% An algorithm for finding the distribution of Maximum Entropy
J. of Comput. Phys. {\bf 30}, (2) 250 (1979)
\bibitem {Mead-Papanicolau} L.R.~Mead, N.~Papanicolaou,
% Maximum entropy in the problem of moments
J. Math. Phis. {\bf 25}, (8) 2404 (1984)
\bibitem {Levine-Tribus} {\it The Maximum Entropy Formalism},
edited by R.D.~Levine \& M.~Tribus, (MIT, Cambridge, MA, 1979)
\bibitem {Jaynes:1982} E.T.~Jaynes, Proc. IEEE {\bf 70}, 939 (1982)
\bibitem {Levine} R.D.~Levine, J. Phys.{\bf A}: Math. Gen. {\bf 13}, 91 (1980)
\bibitem {Burg} J.P.~Burg, in Proc. 37th Annu. Intern. Sos. Explor. Geophys.
(Oklahoma City, OK, 1967)
\bibitem {Shannon} C.~Shannon, Bell Syst. Tech. J. {\bf 27}, 379, 623 (1948)
\bibitem {Jaynes} E.T.~Jaynes, Phys. Rev.  {\bf 106}, 620 (1957)
\bibitem {Kay-Marple}S.M.~Kay, S.L.~Marple Jr., Proc. IEEE {\bf 69},
1380 (1981)
\bibitem {Wernecke} S.J.~Wernecke, Radio Sci. {\bf 12}, (5) 831 (1977)
\bibitem {Notic-Knafel-Turchin-Ugrinovsky-Ugrinovsky}
 A.I.~Notic, A.I.~Knafel, V.I.~Turchin, V.A.~Ugrinovsky, R.A.~Ugrinovsky,
 Radiotecnica i Electronica. {\bf 35}, (9) 1904 (1990) (in Russian)
\bibitem {Korzh-Ovcharenko-Ugrinovsky} S.A.~Korzh, I.E.~Ov\^carenko,
R.A.~Ugrinovsky,  {\it Chebyshev recursion --- some analytical, computation
and applied aspects.} Ukrainian Mathematical Journal
{\bf 45}, (5) (1993)
\bibitem {Levinson} N.~Levinson,
%The Wiener rms (root-mean-square) Error Criterion in Filter Design and
%Prediction, 
J. Math. Physics. {\bf 25}, 261 (1947)
\bibitem {Ugrinovsky} R.A.~Ugrinovsky, Ukrainian Mathematical Journal,
{\bf 44}, (5) (1992)
\bibitem {Arov-Krein} D.Z.~Arov, M.G.~Krein, Functional Anal. and Appl.
{\bf 15}, (2) 123 (1981)
\bibitem {Dym-Gohberg} H.~Dym, I.~Gohberg, Integral Equation and Operator
Theory {\bf 3}, (2) 143 (1980)
\end{thebibliography}

\end{document}