\chapter{Random Processes}


\section{Processes in One Variable}

A wide-sense-stationary random process $v(n)$ may be specified in terms
of its power spectral density $S_v(f)$.  The power spectral density, in
turn, is expressed as the superposition of frequency-shifted copies of a
basis function $\Phi(Mf)$:
\[ S_v(f) = \sum_k w_k \Phi(Mf - \Delta - k) \]
where sequence $w_k$ has period $M$, and $\Delta$ is a frequency offset.
(Sequence $w_k$ is binary valued)

For a process $z(n)$ driving a filter $h(n)$, and output process $u(n) =
(z \ast h)(n)$, the output power is the
value of the output autocorrelation function at offset 0. The output
autocorrelation function is defined as
\[ \mathcal{R}_u(n) = E[ (h \ast z)(k)(z' \ast h')(n-k) ] \] 
Further manipulation results in 
\begin{eqnarray*}
\mathcal{R}_u(n) 
&=& \sum_\ell \sum_m h(\ell) E[z(k-\ell) z'(n-k-m)] h'(m) \\
&=& \sum_\ell \sum_m h(\ell) E[z(w) z'(n-w-\ell-m)] h'(m)
\end{eqnarray*}
Letting $w=(k-\ell)$, 
\begin{eqnarray*}
\mathcal{R}_u(n)
&=& \sum_\ell \sum_m h(\ell) \mathcal{R}_z (n-\ell-m) h'(m) \\
&=& \sum_m (h \ast \mathcal{R}_z)(n-m) h'(m) \\
&=& h \ast \mathcal{R}_z \ast h'
\end{eqnarray*}
and so the output power is
\begin{eqnarray*}
 P_u &=& \mathcal{R}_u(0) \\
&=& (h \ast \mathcal{R}_z \ast h')(0) \\
&=& \sum_\ell \sum_m h(\ell) \mathcal{R}_z(-\ell -m) h'(m) 
\end{eqnarray*}




\opt\ random processes can be added, subtracted, scaled and convolved
with impulse responses and can have their average powers obtained. Here,
\optfn{pwr}(\optvar{s}) and \optfn{pwr}(\optvar{r}) return \ml{0.75} and
\ml{0.5} respectively. If both \optvar{h} and \optvar{d} are fixed,
\begin{code}
\optvar{p} & \optfn{pwr}( \optvar{r} .* ( \optvar{h} - \optvar{d} )) ;
\\
\end{code}
sets \optvar{p} to the numeric MSE in \req{MSE} with $W(f)$ the PSD of
\optvar{r}.  There are no simulations or numerical integrals;
calculation to machine precision is possible because available basis
function exist internal to \opt\ as Fourier transform pairs, and the
required integration is actually an inverse Fourier transform.

Since possible power spectrum basis functions and their Fourier
transforms are known a priori, the autocorrelation function
$\mathcal{R}_z(\tau)$ is easy to compute.



%---------------------
\subsection{Examples}

The following creates a zero-mean, discrete-time \opt\ random process:
\begin{code}
\optvar{rp} & \optfn{Process} ( $\Phi$ , $\Delta$, $\vv{w}$ )
; \\
\end{code}
\vfill
The three arguments to \optfn{Process} specify the PSD of \optvar{s},
where $\Phi$ is the {\em basis function}, $\Delta$ is the frequency
offset and $\vv{w}$ is the weight vector. \opt\ constructs the PSD in
steps:


\begin{enumerate}
\item Weight vector $\vv{w}$ is extended periodically to create a
doubly-infinite sequence $w_k \defn \vv{w}_{1+(k\ \mathrm{mod}\ M)}$,
\item the infinite weight sequence is applied to integer shifts of the
supplied basis function $\Phi(f)$, creating a function $\sum_k
w_k \Phi(f-k)$ with period $M$,
\item this function is offset upward in frequency by $\Delta$, yielding
$\sum_k w_k \Phi(f-\Delta-k)$, and finally
\item scaling the offset function down in frequency by $M$ sets the
period in normalized frequency $f$ to unity:
\begin{equation}
 \mbox{PSD} = \sum_k w_k \Phi(Mf - \Delta - k)
 \label{MSE}
\end{equation}
\end{enumerate}

Here are two example random processes, \optvar{s} and \optvar{r}:
\begin{code}
\optvar{s} & \optfn{Process} ( \optfn{'Box'} , \ml{0} , \ml{[1 1 0 1]} )
; \\
\optvar{r} & \optfn{Process} ( \optfn{'Triangle'} , \ml{0} , \ml{[1 1
3/4 1/4 0 0 0 1/4 3/4 1]} )
; \\
\end{code}

\begin{figure}[h]
%\centerline{\epsfxsize=4in\epsfbox{art/circArray.eps}}
\begin{center}
\includegraphics[width=2.0in,height=1in]{art/prbox.eps}
\quad
\includegraphics[width=2in,height=1in]{art/prtri.eps}
\end{center}
\caption{PSD of random processes.}
\label{art.PSDs}
\end{figure}


The \optfn{'Box'} basis function has unit width, so offset and
scaled copies of it touch without overlap to create an ideal brick-wall
spectrum. To instead create a spectrum linearly interpolated between
samples, use the symmetric-triangle basis function supported on $[-1,
1]$. The PSD of \optvar{r} in \rfig{art.PSDs} approximates a 60\%
raised-cosine spectrum at an oversampling rate of two. A longer weight
vector would improve the approximation.


%-------------

\section{Processes in Several Variables}

\subsection{Average Output Power}

The generalization of an \opt\ random process to several variables is
straightforward. Let $z(\vv{n})$ be a wide-sense-stationary random
process, where $\vv{n} = \{n_1 \ldots n_N \}$.
 Let $h(\vv{n})$ be the impulse response of a linear
shift-invariant (LSI) system.
If $z(\vv{n})$ is the input to this
system, the output random process $u(\vv{n})$ is 
\begin{eqnarray}
 u(\vv{n}) &=& \sum_{k_1} \cdots \sum_{k_N} h(\vv{k}) z(\vv{n}-\vv{k})
		\nonumber \\
	&=& \sum_{\vv{k}} h(\vv{k}) z(\vv{n}-\vv{k}) \label{NDconv} \\
	&=& (h \star z)(\vv{n}) \nonumber
\end{eqnarray}
where the symbol $\vv{k}$ under the sum indicates an $N$-tuple sum over
all possible values of the vector $\vv{k}$.

The average output power, $\mathcal{R}_u(\vv{0})$, is derived as
follows, beginning with the definition
\[ \mathcal{R}_u(\vv{n}; \vv{k}) = E[u(\vv{n}) u^{\ast}(\vv{k})] . \]
Substituting \req{NDconv}, we obtain
\begin{eqnarray*}
 \mathcal{R}_u(\vv{n}; \vv{k}) 
&=& E\left[ \sum_{\vv{l}} \sum_{\vv{m}} h(\vv{l}) z(\vv{n} -
\vv{l}) h^{\ast}(\vv{m}) z^{\ast}(\vv{k} - \vv{m}) \right] \\
&=& \sum_{\vv{l}} \sum_{\vv{m}} h(\vv{l}) h^{\ast}(\vv{m})
	E\left[z(\vv{n} -\vv{l}) z^{\ast}(\vv{k} - \vv{m}) \right]
\end{eqnarray*}
since the impulse response $h(\vv{n})$ is deterministic and can be
removed from the expectation. Using the definition of an autocorrelation
function we can substitute
\begin{eqnarray*} 
  \mathcal{R}_u(\vv{n}; \vv{k})
&=& \sum_{\vv{l}} \sum_{\vv{m}} h(\vv{l}) h^{\ast}(\vv{m})
	\mathcal{R}_z(\vv{n}-\vv{l}; \vv{k}-\vv{m}) \\
&=& \sum_{\vv{l}} \sum_{\vv{m}} h(\vv{l}) h^{\ast}(\vv{m})
	\mathcal{R}_z(\vv{n}-\vv{l}-\vv{k}+\vv{m})
\end{eqnarray*}
since $\mathcal{R}_z(\vv{n}) = E[x(\vv{k})x^\ast(\vv{k}-\vv{n})]$ for
stationary $z$. Since $\mathcal{R}_u$ is shown to be a function of
$\vv{n}-\vv{k}$, it as also stationary, and can be rewritten
\begin{eqnarray*}
\mathcal{R}_u(\vv{n}) 
&=& \sum_{\vv{l}} \sum_{\vv{m}} h(\vv{l}) h^{\ast}(\vv{m}) 
	\mathcal{R}_z (\vv{n} - \vv{l} + \vv{m}) \\
&=& \mathcal{R}_z(\vv{n}) \star h(\vv{n}) \star h^\ast(-\vv{n}) \\
&=& (h \star \mathcal{R}_z \star h')(\vv{n})
\end{eqnarray*}
with $h'$ indicating a conjugate transpose of the sequence $h$.

The average output power is then expressed in a conveniently-computable
form:
\[ \mathcal{R}_u(\vv{0}) 
   = \sum_\vv{l} \sum_\vv{m} h(\vv{l}) h^\ast(\vv{m}) 
	\mathcal{R}_z(\vv{m}-\vv{l}) \]

%------------- symmetry
\subsection{Periodicity and Symmetry Issues}
In the one-dimensional case, a discrete-time LTI system $h(n)$ with
support on the integers has a frequency response $H(f)$ which is
periodic with period 1. Similarly, in multiple dimensions, if an impulse
response $h(\vv{n}), \: \vv{n} = \{n_1, \ldots, n_N\},$
 has support on the square lattice $n = \field{Z}^N$, its frequency
response $H(\vv{f}), \: \vv{f} \in \field{R}^N,$ will be periodic with its
{\em fundamental period} the hypercube with side length 1 centered at
the origin.


%--------------- various basis fcns
\subsection{Basis Functions}
\subsubsection{Box}
As before, this expression for average power is useful since we have
several PSD-autocorrelation transform pairs available. One such PSD
basis function is the multi-dimensional generalization of the scaled `Box'
function $\Box_\vv{w}(\vv{f}) = \prod_k \Box_{w_k}(f_k)$, where each 
$\Box_{w_k}(f_k)$ is a one-dimensional box function, defined as 
\[ \Box_w(f) = \left\{ \begin{array}{cc} 
1 & -w/2 \leq f \leq w/2 \\
0 & \mbox{otherwise} 
\end{array} \right. \]

Several Box functions can be used together to map out regions of the
frequency response to be used as passbands or stopbands.

Since the function $\Box_\vv{w}(\vv{f})$ is separable, its inverse Fourier
transform is the product of one-dimensional inverse Fourier
transforms. So for a WSS random process with PSD $S(\vv{f}) =
\Box_\vv{w}(\vv{f})$, we have
\[ \mathcal{R}(\vv{\tau}) = \prod_{k=1}^N \frac{\sin(\pi w_k \tau_k)}{\pi
\tau_k} \] 
For an offset in frequency $\vv{off}$, the transform pair becomes
\[ \Box_\vv{w}(\vv{f} - \vv{off}) \quad \longleftrightarrow \quad
e^{j2\pi <\vv{\tau}, \vv{off}>}
 \prod_k  \frac{\sin(\pi w_k \tau_k)}{\pi\tau_k} \]
where $<,>$ indicates an inner (dot) product. For a PSD constructed from
a superposition of shifted, scaled boxes, the autocorrelation is
\begin{equation} 
\mathcal{R}(\vv{\tau}) = \sum_\ell c_\ell e^{j2\pi <\vv{\tau}, \vv{off}_\ell>}
	\prod_k \frac{\sin\left(\pi(w_\ell)_k \tau_k\right)}{\pi \tau_k} 
\end{equation}
where $c_\ell$ is the multiplier coefficient of the $\ell$th instance of
the basis function.

\subsubsection{Triangle}
The `Triangle' basis function, $\bigtriangleup_\vv{w}(\vv{f})$, is the
multi-dimensional generalization of the one-dimensional 
\[  \bigtriangleup_w(f) = \left\{ \begin{array}{cc} 
1+f/w & -w \leq f \leq 0 \\
1-f/w & 0 \leq f \leq w \\
0 & \mbox{otherwise} 
\end{array} \right. \]
which is the convolution of two one-dimensional Boxes. As such, a PSD
made up of the superposition of shifted, scaled Triangles is
\begin{equation} 
\mathcal{R}(\vv{\tau}) = \sum_\ell c_\ell e^{j2\pi <\vv{\tau}, \vv{off}_\ell>}
	\prod_k \left( 
	\frac{\sin\left(\pi(w_\ell)_k \tau_k\right)}{\pi \tau_k} 
	\right)^2
\end{equation}

\subsubsection{Impulse}
The `Impulse' basis function, $\delta(\vv{f})$, is the
multi-dimensional generalization of the one-dimensional $\delta(f)$. It
is useful to use a grid of Impulses to map out an irregularly-shaped
band of the desired filter's frequency response. For a PSD composed of a
single Impulse, the PSD / autocorrelation pair is
\[ S(\vv{f}) = \delta(\vv{f}) \quad \longleftrightarrow \quad 
	\mathcal{R}(\vv{\tau}) = 1 \]
A PSD composed of the superposition of several shifted, scaled impulses
therefore has the autocorrelation function
\[ \mathcal{R}(\vv{\tau}) = \sum_\ell c_\ell 
	e^{j2\pi <\vv{\tau},\vv{off}>} \]

\subsubsection{Circle}
The `Circle' basis function, available in a two-dimensional version, is
defined as
\[  \bigcirc_w\left(f_1, f_2\right) = 
\left\{ \begin{array}{cc} 
1 & \sqrt{f_1^2 + f_2^2} \leq w \\
0 & \mbox{otherwise} 
\end{array} \right. \]
This function has circular symmetry and can be written as a function of
$\rho = \sqrt{f_1^2 + f_2^2}$, and therefore its inverse fourier
transform can be found using the Fourier-Bessel transform formula,
\[ g(r) = 2\pi \int_0^\infty \rho G(\rho) J_0(2\pi r \rho) \: d\rho \]
where $J_\nu(a)$ is a Bessel function of the first kind of order $\nu$.
The inverse transform, as a function of $r = \sqrt{\tau_1^2 + \tau_2^2}$,
is
\begin{eqnarray*}
 \mathcal{R}(r) 
&=& 2\pi \int_0^\infty \rho \bigcirc_w(\rho) J_0(2 \pi r \rho) \: d\rho
\\
&=& 2\pi \int_0^w \rho J_0(2 \pi r \rho) \: d\rho
\end{eqnarray*}
Letting $u = 2\pi r \rho$, and thus $du = 2\pi r \: d\rho$, we have
\[
 \mathcal{R}(r) 
= \frac{1}{2\pi r^2} \int_0^{2 \pi r w} J_0(u) \: u \: du \]
Using the identity
\[ \int_0^u u' J_0(u') \: du' = u J_1(u) , \]
we have 
\[
 \mathcal{R}(r) 
= w \frac{J_1(2\pi r w)}{r} . \]
and thus
\[ \mathcal{R}(\tau_1, \tau_2) 
= w \frac{J_1(2\pi w \sqrt{\tau_1^2 + \tau_2^2} )}
	 {\sqrt{\tau_1^2 + \tau_2^2} } \]

For a PSD composed of the superposition of shifted Circles of various
radii the autocorrelation is
\[ \mathcal{R}(\vec{\tau}) = \sum_\ell c_\ell w_\ell
	e^{j2\pi <\vec{\tau},\vv{off}>} 
   \frac{J_1\left(2\pi ||\vec{\tau}|| w_\ell \right) }{ ||\vec{\tau}|| }
\]

A useful limit is
\[ \lim_{||\vec{\tau}|| \to 0} w  
\frac{J_1\left(2\pi ||\vec{\tau}|| w \right) }{ ||\vec{\tau}|| } 
= \pi w^2 \]

\subsection{Examples}

A multi-dimensional random process is created using  
\begin{code}
\optvar{rp} & \optfn{NDProcess} ( \optvar{basis}, \optvar{param},
\optvar{freqs}, \optvar{coeff} )
; \\
\end{code}
Here \optvar{rp} is created with a PSD consisting of a superposition of
basis functions at shifts given in \optvar{freqs} and with weights given
in \optvar{coeff}. 
\optvar{freqs} is a matrix in which each row is a vector of frequency
offsets; therefore, \optvar{freqs} should contain as many rows as
shifted basis functions, and as many columns as the dimension of the
random process. \optvar{coeff} is a vector of length equal to the number
of rows of \optvar{freqs}.

If the PSD is to be made up entirely of a single type
of basis function, \optvar{basis} should be the name of the basis
function; for example, \ml{'Box'}. If the PSD is to be composed of
several types of basis functions, \optvar{basis} must be a column vector
cell array of length equal to the number of rows of \optvar{freqs}; that
is, one entry for each shifted basis function. For example:
\ml{ \{ 'Box' ; 'Impulse' ; 'Impulse' \} }.

\optvar{param} contains the parameters for each shifted copy of a basis
function. For a single-type PSD, \optvar{param} should be a matrix, each
row of which containing the parameters for another copy of the basis
function. For a multiple-type PSD, since different basis function types
require different parameters to be specified, \optvar{param} must be a
cell array column vector, in which each entry contains the parameters
for another shifted basis function. 

The basis function types \ml{'Box'} and \ml{'Triangle'}
 in $N$ dimensions require one width
parameter for each dimension. The \ml{'Impulse'} type requires no
parameters, and the empty matrix [] should be used in place of any
parameter. The \ml{'Circle'} basis function requires a radius parameter.



Here are two example random processes:
\begin{code}
\optvar{s} & \optfn{NDProcess} ( \ml{ \{'Box'; 'Circle'; 'Impulse'\} },
\ml{ \{[1 1]; .2; []\} }, \ml{[0 0; 0 0; 0 0]}, \ml{[1; -1; 1]}); \\
\optvar{r} & \optfn{NDProcess} ( \optfn{'Box'} , 
\ml{[.1 .1 .1; .1 .1 .1]}  , \ml{[0 .25 .25; 0 .25 -.25]} , \ml{[1 1]} )
; \\
\end{code}

Random process \optvar{s} is two-dimensional. It consists of three basis
functions, a Box, a Circle and an Impulse. The Box spans the entire
fundamental period; it has width 1 in each dimension and is centered at
the origin. The Circle is of radius .2 and centered at the
origin. Because the Circle's coefficient is -1 the effect is of a
rectangle with a circle cut out of it. An impulse is added back in at
the origin. Note that the \optvar{param} entries must be split into a
cell array vector because the PSD consists of different basis function
types. 

Process \optvar{r} is three-dimensional, and is composed only of
Box-type basis functions. Since it is a single-type PSD, the
\optvar{param} entry is a matrix. The \optvar{freqs} matrix implies a 3-D
process by its width, and that there are two copies of a Box by its height.
