The ideas behind Autoregressive model are considered in this digest. Examples of describing the autoregressive model are evaluated.
Autoregressive model of a stochastic signal in a nutshell
Assume we have a noise-alike discrete sequence $y[n]$ that we want to study. We can describe the signal in terms of autocorrelation or power spectrum density, but we can also describe it in the other way. We can actually describe the noised signal $y[n]$ as a result of passing a white noise $x[n]$ through a filter $h[n]$. The filter $h[n]$ can be described in Z-domain as following:
$H(z) = \frac{B(z)}{A(z)} = \frac{\sum\limits_{k=0}^{q} b_k z^{-k} }{1+ \sum\limits_{k=0}^{p} a_k z^{-k} }$
that is a stable shift-invariant linear system with $p$ poles and $q$ zeros. The autoregressive (AR) model is a case of the ARMA and is the following:
$H(z) = \frac{b_0}{1+ \sum\limits_{k=0}^{p} a_k z^{-k} }$
where $a_1, \dots, a_p$ are the parameters of the model. We can think\cite{ramadanadaptivefiltering} of the signal $y[n]$ being studied as a result of passing the white noise $x[n]$ with known variance $\sigma_x^2$ through the filter $H[z]$.
The autocorrelation function $y[n]$ and the parameters of the filter $H[z]$ are related via Yule-Walker equation:
$R_{yy}[k] + \sum_{m=1}^p a_m R_{yy}[k-m] = \sigma_x^2 b_0^2 \delta_{m,0},$
where $m = 0, ... , p,$ yielding $p+1$ equations. Those equations are usually written in matrix form:
$\left[ \begin{matrix}
R_{yy}(0) & R_{yy}(-1) & \dots R_{yy}(-p) \\
R_{yy}(1) & R_{yy}(0) & \dots R_{yy}(-p+1) \\
\vdots & \vdots & \vdots \\
R_{yy}(p) & R_{yy}(p-1)& \dots R_{yy}(0) \\
\end{matrix} \right]
\left[ \begin{matrix}
1 \\ a_1 \\ \vdots \\ a_p
\end{matrix}\right]
= \sigma_x^2 b_0^2
\left[ \begin{matrix}
1 \\ 0 \\ \vdots \\ 0
\end{matrix}\right]$
and we need to solve this matrix equation for coefficients $a_1 \dots a_p$. The above equations provide the way for estimation of the parameters of an AR(p) model. In order to solve the matrix equation, we must replace the theoretical covariances with estimated values such as $R_{yy}(0)$.
A note of Toeplitz matrices
When we study a signal and calculate its autocorrelation sequence, we got a vector. In order to use these data in solving the matrix equation, we need to convert the autocorrelation vector into matrix. That can be performed by Toeplitz matrices. A Toeplitz matrix is a matrix in which each descending diagonal from left to right is constant:
$\begin{bmatrix} a & b & c & d & e \\ f & a & b & c & d \\ g & f & a & b & c \\ h & g & f & a & b \\ i & h & g & f & a \end{bmatrix}.$
A Toeplitz matrix is defined by one row and one column. A symmetric Toeplitz matrix is defined by just one row. In order to generate the Toeplitz matrix, one can use \textbf{toeplitz} function in MATLAB that generates Toeplitz matrices given just the row or row and column description.
In the example below we will generate the autocorrelation matrix from the vector using Toeplitz matrices. Consider we have the autocorrelation vector of the observed signal $y[n]$ over the lag range \textit{[-maxlags:maxlags]}. The output vector $r_{yy}$ will have the length \textit{2*maxlags+1}:
\begin{verbatim}
[r_yy, lags] = xcorr(y,y,xcorr_maxlags,'biased');
\end{verbatim}
The output will be autocorrelation vector like this:
$r_{yy} = 0.7066 \,\,\,\, 0.7116 \,\,\,\, 0.7176 \,\,\,\,0.7147 \,\,\,\, 0.7018 \,\,\,\, 0.6954 \,\,\,\, 0.7029 \,\,\,\, 0.7208 \,\,\,\, 0.7323 \,\,\,\, 0.7245$
Then we need to make a autocorrelation matrix - use the Toeplitz matrices:
\begin{verbatim}
R = toeplitz(r_yy(1,1:acorr_matrix_size)); %% create the autocorr. matrix
\end{verbatim}
The result will be like:
$R =
\left[ \begin{matrix}
0.7066& 0.7116& 0.7176& 0.7147& 0.7018& 0.6954 \\
0.7116& 0.7066& 0.7116& 0.7176& 0.7147& 0.7018 \\
0.7176& 0.7116& 0.7066& 0.7116& 0.7176& 0.7147 \\
0.7147& 0.7176& 0.7116& 0.7066& 0.7116& 0.7176 \\
0.7018& 0.7147& 0.7176& 0.7116& 0.7066& 0.7116 \\
0.6954& 0.7018& 0.7147& 0.7176& 0.7116& 0.7066 \\
\end{matrix} \right]$
That is how we can make the autocorrelation matrix from the autocorrelation vector.
Example of the Autoregressive model
Consider the example in MATLAB from the book \cite{ramadanadaptivefiltering}. First of all, let's generate the signal $y[n]$ from the white noise (it is a simulation of the real signal).
\begin{verbatim}
x = 1.0*(rand(1,signal_length) - 0); % x is a White Gaussian Noise.
y = filter(1, [1 -0.9 0.5], x); % here y is observed (filtered) signal.
\end{verbatim}
The \textbf{filter} function filters \textbf{y = filter(b,a,X)} the data in vector \textbf{X} with the filter described by numerator coefficient vector \textbf{b} and denominator coefficient vector \textbf{a} that is:
$y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) - a(2)*y(n-1) - ... - a(na+1)*y(n-na)$
Any WSS process $y[n]$ can be represented as the output of filter $h[n]$ that is driven by noise x[n]. Using that filter coefficients (arbitrary chosen - we just need to generate the signal), we design the signal which we are going to study. Then we describe the signal $y[n]$ in terms of filter function defined by AR coefficients.
In order to solve Yule-Walker equations, we need to find the autocorrelation function of the signal $y[n]$:
\begin{verbatim}
[r_yy, lags] = xcorr(y,y,xcorr_maxlags,'biased'); %% autocorrelation vector of the observed signal y[n]
R = toeplitz(r_yy(1,1:acorr_matrix_size)); %% create the smaller autocorrelation matrix
\end{verbatim}
We used the Toeplitz matrix, as described above, to obtain the autocorrelation matrix. The variable \verb|acorr_matrix_size| is for setting the size of the autocorrelation matrix, which is also setting the order of the filter H(z):
$H(z) = \frac{b_0}{1+ \sum\limits_{k=0}^{p} a_k z^{-k} }$
OK, now we need to find the coefficients $a_1 \dots a_p$ so we need to find the inverse autocorrelation matrix:
\begin{verbatim}
R_inv = pinv(R); %% R_inv is the inverse of R
\end{verbatim}
First element of the $R_{yy}$ is the $b_0^2$ coefficient that is needed for the filter $H(z)$:
\begin{verbatim}
b0_squared = 1/R_inv(1,1); %% this is a variance of the process R_yy(0)
b0 = sqrt(b0_squared); %% find b0 that is \sqrt{R_yy(0)}.
\end{verbatim}
Using these data, we can find the AR coefficients $a_1 \dots a_p$:
\begin{verbatim}
a = b0_squared*R_inv(2:acorr_matrix_size,1); %%% find the AR coefficients
\end{verbatim}
Then we should denote the order of the filter that is $H(z) = \frac{b_0}{1+a_1z_{-1}+...+a_kz^{-k}}$ by denoting the number of autocorrelation coefficients used. Finally, we can calculate the filter:
\begin{verbatim}
H = b0./(fft([1 a'],filter_fft_length)); %% H is the frequency response of the filter
\end{verbatim}
Now we can describe the process y[n] in terms of the autoregressive model coefficients $a_1 \dots a_p$.
No comments:
Post a Comment