paper.tex

国外免费地震资料处理软件包

\lefthead{Sava \& Fomel}
\righthead{Spectral factorization}
\footer{SEP--100}

\title{Spectral factorization revisited}

\email{[e-mail], [e-mail]}
\input{macros}

\author{Paul Sava and Sergey Fomel}

\maketitle

\begin{abstract}
In this paper, we review some of the iterative methods for the square
root, showing that all these methods belong to the same
family, for which we find a general formula. We then explain how those 
iterative methods for real numbers can be extended to spectral
factorization of auto-correlations. The iteration based on
the Newton-Raphson method is optimal from the convergence stand point, though
it is not optimal as far as stability is concerned. Finally, we show
that other members of the iteration family are more stable, though
slightly more expensive and slower to converge.
\end{abstract}

\inputdir{Matlab}

\section{Introduction}
Spectral factorization has been recently revived by the advent of the
helical coordinate system. Several methods are reported in the
literature, ranging from Fourier domain methods, such as Kolmogoroff's
\cite[]{Claerbout.blackwell.92, kolmog}, to iterative methods, such as the
Wilson-Burg method \cite[]{gee,mywilson,Sava.sep.97.paul1}.
\par
In this paper, after reviewing the general theory of root estimation by
iterative methods, we derive a general square root relationship
applicable to both real numbers and to auto-correlation functions. We
introduce a new spectral factorization relation and show 
its relation to the Wilson-Burg method.

\section{The square root of real numbers}
This section briefly reviews some well known square root iterative
algorithms, and derives the Newton-Raphson and Secant methods. It also
shows that Muir's iteration for the square root \cite[]{Claerbout.bei.95}
belongs to the same family of iterative methods, if we make an
appropriate choice of the generating function.

\subsection{Root-finding recursions}
Given a function $f(x)$ and an approximation for one of its roots
$x_n$, we can find a better approximation for the root by linearizing
the function around $x_n$
\[f(x) \approx f(x_n)+(x_{n+1}-x_n) f'(x_n) \]
and by setting $f(x)$ to be zero for $x=x_{n+1}$. We find that
\beq \label{eqn:newton-raphson}
\fbox {$ \displaystyle 
x_{n+1}=x_{n}-\frac{f(x_n)}{f'(x_n)} 
$} \eeq
\begin {enumerate}
%
\item {\bf Newton-Raphson's method for the square root} 
\par
A common choice of the function $f$ is $f(x)=x^2-s$. This
function has the advantage that it is easily differentiable, with
$f'(x)=2x$. The recursion relation thus becomes
\[ x_{n+1}=x_{n}-\frac{x_n^2-s}{2x_n}
=\frac{x_n}{2}+\frac{s}{2x_n} \]
or
\[ x_{n+1}=\frac{1}{2}\left(x_n+\frac{s}{x_n}\right) \]
or, after rearrangement,
\beq \label{eqn:newton}
\fbox {$ \displaystyle
x_{n+1}=\frac{s+x_n^2}{2 x_n}  
$} \eeq
\par
The recursion (\ref{eqn:newton}) converges to $\pm \sqrt{s}$ depending
on the sign of the starting guess $x_0\ne 0$.
%
\item {\bf Secant method for the square root} 
\par
A variation of the Newton-Raphson method is to use a finite
approximation of the derivative instead of the differential form. In
this case, the approximate value of the derivative at step $n$ is
\[ f'(x_n)=\frac{f(x_n)-f(x_{n-1})}{x_{n}-x_{n-1}} \]
\par
For the same choice of the function $f$, $f(x)=x^2-s$, we
obtain
\[ x_{n+1}=x_{n}-\frac{x_n^2-s}{x_n+x_{n-1}}\]
and
\beq \label{eqn:secant}
\fbox {$ \displaystyle
x_{n+1}=\frac{s+x_{n}x_{n-1}}{x_n+x_{n-1}} 
$} \eeq
\par
In this case, recursion (\ref{eqn:secant}) also converges to $\pm
\sqrt{s}$ depending on the sign of the starting guesses $x_0$ and
$x_1$.
%
\item {\bf Muir's method for the square root} 

\par
Another possible iterative relation for the square root is Francis
Muir's, described by \cite{Claerbout.bei.95}:
\beq \label{eqn:muir}
\fbox {$ \displaystyle 
x_{n+1}=\frac{s+x_{n}}{x_n+1} 
$} \eeq
\par
This relation belongs to the same family of iterative schemes as
Newton and Secant, if we make the following special choice of the
function $f(x)$ in (\ref{eqn:newton-raphson}):
\beq \label{eqn:muf} 
f(x)= 
|x+\sqrt{s}|^{\frac{\sqrt{s}-1}{2\sqrt{s}}}
|x-\sqrt{s}|^{\frac{\sqrt{s}+1}{2\sqrt{s}}} \eeq
\par
Figure~\ref{fig:muf} is a graphical representation of the
function f(x).
%%
\sideplot{muf}{width=3.0in}{The graph of the function defined in
Equation~(\ref{eqn:muf}) used to generate Muir's iteration for the
square root (solid line). The dashed lines are the plot of the two
factors in the equation~\label{eqn:muf}.}
%%
\item {\bf A general formula for the square root}
\par
From the analysis of equations (\ref{eqn:newton}), (\ref{eqn:secant}),
and (\ref{eqn:muir}), we can derive the following general form for 
the square root iteration:
\beq \label{eqn:general}
\fbox {$ \displaystyle 
x_{n+1}=\frac{s+x_{n}\gamma}{x_n+\gamma} 
$} \eeq
where $\gamma$ can be either a fixed parameter, or the value of the
iteration at the preceding step, as shown in Table~\ref{recursion}. 
%%
\begin{table} \center
\caption{Recursions for the square root}
\label{recursion}
\begin{tabular}{|c|c|c|}\hline
\R        & $\gamma$   & Recursion \\\hline
\R Muir   & $1$        & $x_{n+1}=\frac{s+x_n        }{x_n+1       }$ \\\hline
\R Secant & $x_{n-1}$  & $x_{n+1}=\frac{s+x_n x_{n-1}}{x_n+x_{n-1} }$ \\\hline
\R Newton & $x_n$      & $x_{n+1}=\frac{s+x_n^2      }{2  x_n      }$ \\\hline
\R Ideal  & $\sqrt{s}$ & $x_{n+1}=\frac{s+x_n\sqrt{s}}{x_n+\sqrt{s}}$ \\\hline
\end{tabular}
\end{table}
%%\beq \label{eqn:recursion}
%%\begin{array}{|c|c|c|}          \hline
%%\R        & \gamma   & Recursion \\\hline
%%\R Muir   & 1        & x_{n+1}=\frac{s+x_n        }{x_n+1       } \\\hline
%%\R Secant & x_{n-1}  & x_{n+1}=\frac{s+x_n x_{n-1}}{x_n+x_{n-1} } \\\hline
%%\R Newton & x_n            & x_{n+1}=\frac{s+x_n^2      }{2  x_n      } \\\hline
%%\R Ideal  & \sqrt{s} & x_{n+1}=\frac{s+x_n\sqrt{s}}{x_n+\sqrt{s}} \\\hline 
%%\end{array} \eeq
%%
The parameter $\gamma$ is the estimate of the square root at the given
step (Newton), the estimate of the square root at the preceding step
(Secant), or a constant value (Muir). Ideally, this value should be as
close as possible to $\sqrt{s}$.
\end{enumerate}

\subsection{The convergence rate}
We can now analyze which of the particular choices of $\gamma$ is more
appropriate as far as the convergence rate is concerned.
\par
If we consider the general form of the square root iteration
\[x_{n+1}=\frac{s+x_{n}\gamma}{x_n+\gamma}\]
we can estimate the convergence rate by the difference between the
actual estimation at step $(n+1)$ and the analytical value
$\sqrt{s}$. For the general case, we obtain
\[ x_{n+1}-\sqrt{s} =
\frac{s+\gamma x_n -x_n \sqrt{s}-\gamma \sqrt{s}}{x_n+\gamma}\]
or
\beq \label{convergence}
\fbox {$ \displaystyle
x_{n+1}-\sqrt{s} =
\frac{(x_n-\sqrt{s})(\gamma-\sqrt{s})}{x_n+\gamma} 
$} \eeq
%%
\sideplot{sqroot}{width=3.0in}{Convergence plots for different
recursive algorithms, shown in Table~\ref{recursion}.}
%%
The possible selections for $\gamma$ from Table~\ref{recursion} 
clearly show that the recursions
described in the preceding subsection generally have a linear
convergence rate (that is, the error at step $n+1$ is proportional to
the error at step $n$), but can converge quadratically for an
appropriate selection of the parameter $\gamma$, as shown in 
Table~\ref{convergence}. Furthermore, the convergence is faster
when $\gamma$ is closer to $\sqrt{s}$.

\begin{table} \center
\caption{Convergence rate}
\label{convergence}
\begin{tabular}{|c|c|c|}\hline
\R        & $\gamma$   & Convergence     \\\hline
\R Muir   & $1$        & linear          \\\hline
\R Secant & $x_{n-1}$  & quasi-quadratic \\\hline
\R Newton & $x_n$      & quadratic       \\\hline
\end{tabular}
\end{table}
%%\beq \label{eqn:convergence}
%%\begin{array}{|c|c|c|}          \hline
%%\R        & \gamma   &        Convergence     \\\hline
%%\R Muir   & 1        &     linear       \\\hline
%%\R Secant & x_{n-1}  &  quasi-quadratic \\\hline
%%\R Newton & x_n            &   quadratic      \\\hline
%%\end{array} \eeq
\par
We therefore conclude that Newton's iteration has the potential to
achieve the fastest convergence rate. Ideally, however, we could use
a fixed $\gamma$ which is a good approximation to the square root. The
convergence would then be slightly faster than for the Newton-Raphson 
method, as shown in Figure~\ref{fig:sqroot}.

\section{Spectral factorization}
We can now extend the equations derived for real numbers to
polynomials of Z, with $Z=e^{i\omega t}$, and obtain spectral
factorization algorithms similar to the Wilson-Burg method
\cite[]{Sava.sep.97.paul1}, as follows:
\beq \label{eqn:genspecfac}
\fbox {$ \displaystyle 
X_{n+1}=\frac{S+X_{n} \bar{G}}{ \bar{X_n} + \bar{G}} 
$} \eeq
\par
If $L$ represents the limit of the series in (\ref{eqn:genspecfac}),  
\[ L \bar{L} + L \bar{G} = S + L \bar{G} \] 
and so
\[ L \bar{L} = S\]
\par
Therefore, $L$ represents the causal or anticausal part of the given
spectrum $S=X\bar{X}$.
\par
Table~\ref{specfac} summarizes the spectral factorization
relationships equivalent to those established for real numbers in
Table~\ref{recursion}.
%%
\begin{table} \center
\caption{Spectral factorization}
\label{specfac}
\begin{tabular}{|c|c|}\hline
\R General& $X_{n+1} =\frac{S+X_n \bar G      }{ \bar X_n+\bar G      }$\\\hline
\R Muir   & $X_{n+1} =\frac{S+X_n             }{ \bar X_n+1           }$\\\hline
\R Secant & $X_{n+1} =\frac{S+X_n \bar X_{n-1}}{ \bar X_n+\bar X_{n-1}}$\\\hline
\R Newton & $X_{n+1} =\frac{S+X_n \bar X_n    }{2\bar X_n             }$\\\hline
\R Ideal  & $X_{n+1} =\frac{S+X_n\sqrt{S}     }{ \bar X_n+\sqrt{S}    }$\\\hline
\end{tabular}
\end{table}
%%\beq \label{eqn:specfac}
%%\begin{array}{|c|c|}    \hline
%%\R General& X_{n+1} =\frac{S+X_n \bar G      }{ \bar X_n+\bar G      }\\\hline
%%\R Muir   & X_{n+1} =\frac{S+X_n             }{ \bar X_n+1           }\\\hline
%%\R Secant & X_{n+1} =\frac{S+X_n \bar X_{n-1}}{ \bar X_n+\bar X_{n-1}}\\\hline
%%\R Newton & X_{n+1} =\frac{S+X_n \bar X_n    }{2\bar X_n             }\\\hline
%%\R Ideal  & X_{n+1} =\frac{S+X_n\sqrt{S}     }{ \bar X_n+\sqrt{S}    }\\\hline
%%\end{array} \eeq
\par
The convergence properties are similar to those derived for
real numbers. As shown above, the Newton-Raphson method should have
the fastest convergence. 
%%The only better choice is a $G(Z)$, which is a good
%%approximation for the causal/anticausal part of the spectrum $S(Z)$.

\section{A comparison with the Wilson-Burg method}
For reasons of symmetry, we can take Newton's relation from
Table~\ref{specfac}
\[ X_{n+1} =\frac{S+X_n \bar X_n}{2\bar X_n}\]
and convert it to
\[ \frac{X_{n+1}}{2 X_n} =
\frac{S+X_n \bar X_n}{(2 X_n) (2\bar X_n)}. \]
We can then consider a symmetrical relation where on the left side we
insert the anticausal part of the spectrum, and obtain
\[ \frac{\bar X_{n+1}}{2 \bar X_n} =
\frac{S+X_n \bar X_n}{(2 X_n) (2\bar X_n)}. \]
Finally, we can sum the preceding two equations and get
\beq \fbox {$ \displaystyle
\frac{X_{n+1}}{2 X_n} + \frac{\bar X_{n+1}}{2 \bar X_n} = 
\frac{2S+ X_n \bar X_n + \bar X_n X_n}{(2 X_n) (2\bar X_n)}
$} \eeq
which can easily be shown to be equivalent to the Wilson-Burg
relation
\beq \label{eqn:wilsonburgSpF} 
\frac{X_{n+1}}{X_n} + \frac{\bar X_{n+1}}{\bar X_n} =
1 + \frac{S}{ X_n \bar X_n}
\eeq
\par
In an analogous way, we can take the general relation from
Table~\ref{specfac}
\[ X_{n+1} =\frac{S+X_n \bar G}{\bar X_n+\bar G} \]
and convert it to
\[ \frac{X_{n+1}}{X_n + G} =
\frac{S+X_n \bar G}{(X_n + G) (\bar X_n + \bar G)}\]
We can then consider a symmetrical relation where on the left side we
insert the anticausal part of the spectrum, and obtain
\[ \frac{\bar X_{n+1}}{\bar X_n + \bar G} =
\frac{S+\bar X_n G}{(X_n + G) (\bar X_n + \bar G)}\]
Finally, we can sum the preceding two equations and get
\beq \label{eqn:generalSpF} \fbox {$ \displaystyle
\frac{X_{n+1}}{X_n + G} + \frac{\bar X_{n+1}}{\bar X_n + \bar G} = 
\frac{2S+X_n \bar G + \bar X_n G}{(X_n + G) (\bar X_n + \bar G)}
$} \eeq
\par
Equation~(\ref{eqn:generalSpF}) represents our general formula for
spectral factorization. If we consider the particular case when $G$ is 
$X_n$, we obtain equation~(\ref{eqn:wilsonburgSpF}), which we have shown
to be equivalent to the Wilson-Burg formula.
\par
From the computational standpoint, our equation is more expensive than
the Wilson-Burg because it requires two more convolutions on the
numerator of the right-hand side. However, our equation offers more
flexibility in the convergence rate. If we try to achieve a quick
convergence, we can take $G$ to be $X_n$ and get the Wilson-Burg
equation. On the other hand, if we worry about the stability, 
especially when some of the
roots of the auto-correlation function are close to the unit circle,
and we fear losing the minimum-phase property of the factors, 
we can take $G$ to be some damping function, more tolerant of
numerical errors.
\par
Moreover, by using the Equation~(\ref{eqn:generalSpF}), we can achieve
fast convergence in cases when the auto-correlations we are
factorizing have a very similar form, for example, in nonstationary
filtering. In such cases, the solution at the preceding step can be
used as the $G$ function in the new factorization. Since $G$ is already
very close to the solution, the convergence is likely to occur quite
fast.

\section{Conclusions}
The general iterative formula for the square root that we derived can
be extended to the factorization of the auto-correlation
functions. The Wilson-Burg algorithm is a special
case of our more general formula. Using such a general formula
provides flexibility in choosing between fast convergence and
stability. We can achieve fast convergence
when factorizing auto-spectra that have a very similar form.
This improvement in convergence rate can have a useful
application, for instance, in nonstationary preconditioning.

\section{Acknowledgments}
We thank Jon Claerbout, who brought Muir's iterative scheme to our
attention, and suggested its application to spectral factorization.

\bibliographystyle{seg}
\bibliography{SEG,SEP2,paper}