# Confusion about ensembles and averages in autocorrelation matrices

I just noticed that until now I often don't cared about the scaling of the autocorrelation Matrix in Matlab. I then checked with the book Statistical Digital Signal Processing from M. H. Hayes. There it is defined like this: First of all, why does the author note it like $$x(0)x^*(0)$$? Written like this wouldm't this mean first element of the vector x times first element of the vector x? And the expected value is then what $$\frac{x(0)x^*(0)}{2}$$?

$$r_x(\tau)$$ can be calculated as $$E[X(t)X(t+\tau)]$$ This is more what I would expect. A pointwhise vector multiplication of a $$\tau$$ shifted vector.

Here comes my actual question. For e.g. $$x = [1 ,2, 3]^T$$, $$r_x(0)$$ should become $$\frac{1^2+2^2+3^2}{3}=\frac{14}{3}$$ in my oppinion. Matlab shows for xcorr(x)

3.0000    8.0000   14.0000    8.0000    3.0000


$$r_x(0)$$ should be the value in the middle. But why don't they take the expected value? On this page they suggest that the autocorrelation matrix can be calculated with

z=autocorr(x)


and

Rxx=toeplitz(z)


Instead of autocorr, xcorr(x) should work too stated other posts. As expected this results in

Rx =

3.0000    8.0000   14.0000    8.0000    3.0000
8.0000    3.0000    8.0000   14.0000    8.0000
14.0000    8.0000    3.0000    8.0000   14.0000
8.0000   14.0000    8.0000    3.0000    8.0000
3.0000    8.0000   14.0000    8.0000    3.0000


Again, no expected values. And this solution is also wrong because it is shifted. $$r_x(0)$$ is not in the upper left corner. I could shift it by hand, but is that how autocorrelation matrixes are calculated in Matlab? I checked with some other sources in the internet. Often the algorithm to calculate an entry of an autocorrelation contains just the sum.

• one averages over an ensemble of samples, not lags – user28715 Feb 3 '19 at 22:21

EDIT: sorry that I didn't fully get your question. Now I got your points of conufion and try to provide a shortest possible answer:

First, you should better use xcorr function when estimating ACS $$r_x[k]$$ from data sequence $$x[n]$$. The output of xcorr will return the lag-0 at the middle sample as you recognize. Look at the input argument scaleopt which provides normalization to the computed output. The default is none ! and causes a confusion for you. You should select biased so that equation 7 below and xcorr output will be identical. Like this:

rx = xcorr([1,2,3],'biased'); % identical to (1/3)*xcorr([1,2,3])


Second, do not use the raw output of xcorr as the input to the toeplitz function. As the output of xcorr will include lags for rx[-2],rx[-1],rx,rx,rx. But you only need rx,rx,rx. So you better call it like:

L = 3;
x = [1,2,3];
rxL = xcorr( x, 'biased');
Rx = toeplitz( rxL(L:end) );  % now yo uare ok ;-)


Now it could be argued that the Matlab's default xcorr scale option none is misleading. And yes, I would say so. In its default form it computes the convolution and leaves it unscaled...

Finally, please distingusih a theoretical expectation from practical estimation. When algorithms work on data, they estimate and they don't compute expectations, which is a theoretical tool. Hence, all expectations are replaced with estimations when algorithms are run for the practical applications... It seems (to me) that your confusion is about three things.

• Theoretical definition of the autocorrelation sequence $$r_{XX}[k]$$ of a random process $$X[n,s]$$
• Practical estimation of it from a sample function $$x[n]$$ via ergodicity.
• Interpretation of a random vector (including samples of a random process)

Consider the interpretation of a random process $$X[n,s]$$ as an ordered collection (set) of indexed random variables $$X_n(s)$$ for each time index $$n$$. Each random variable $$X_n$$ may, in principle, has its own unique probability distribution/density function $$F_{X_n}(x)$$. Consider the enormous implication of this on the number and type of all possible joint density functions.

To make things extremely simplified, we make the assumptions that all those random variables have the same PDF $$F_{X_n}(x)$$ for all $$n$$, and independent of each other, also further that first and second moments are independent of partciular time instants, which yields an i.i.d. , WSS random process $$X[n,s]$$ at hand.

The theoretical auto-correlation sequence (ACS) of this WSS random process is defined as:

$$r_{X_n X_m}(n,m) = E\{X_n X^*_m\} = E\{X[n]X^*[n+k] \} = r_{X X}[k] \tag{1}$$

The computation of this theoretical ACS $$r_x[k]$$ requires you to evaluate the integrals with PDF functions. This is possible with questions on the paper, but with applications in working systems, where you won't solve (mathematical) questions but perform statistical calculations, you need to estimate it from data, based on ergodicty assumption, which states that theoretical ensemble averages can bu estimates from practical sample function time-averages. Then you have

$$\hat{r}_X[k] = \lim_{N \to \infty} \frac{1}{2N+1} \sum_{n=-N}^{N} x[n]x^*[n+k] \tag{2}$$

Note that $$x[n]$$ is a particular sample function (a realization) of the hrandom process $$X[n,s]$$. Note also that the infinite limits are replaced by data observation boundaries in practice.

Now consider the random vector $$X(n) = [X(n),X(n-1),..,X(n-N+1)]^T$$ which includes $$N$$ consecutive random variables from the random process $$X[n,s]$$, beginning at index $$n$$.

From this vector build up the outer product $$X \cdot X^H$$ , which includes a matrix of pairwise multiplications of random variables $$X_i X_j$$ as in

$$RM = \begin{bmatrix} X_1 X^*_1 & X_1 X^*_2 & ... & X_1 X^*_N \\ X_2 X^*_1 & X_2 X^*_2 & ... & X_2 X^*_N \\ ... & ... & ... & ... \\ X_N X^*_1 & X_N X^*_2 & ... & X_N X^*_N \\ \end{bmatrix} \tag{3}$$

Now taking the expectation of this $$RM$$ matrix yields the theoretical auto-correlation matrix of the WSS random process as:

$$R_X = E\{ X \cdot X^H \} = E\{RM\} = E\{\begin{bmatrix} X_1 X^*_1 & X_1 X^*_2 & ... & X_1 X^*_N \\ X_2 X^*_1 & X_2 X^*_2 & ... & X_2 X^*_N \\ ... & ... & ... & ... \\ X_N X^*_1 & X_N X^*_2 & ... & X_N X^*_N \\ \end{bmatrix} \} \tag{4}$$

and this results in: $$R_X = E\{ X \cdot X^H \} = \begin{bmatrix} r_x & r_x[-1] & ... & r_x[-N+1] \\ r_x & r_x & ... & r_x[-N+2] \\ ... & ... & ... & ... \\ r_x[N-1] & r_x[N-2] & ... & r_x \\ \end{bmatrix} \tag{5}$$

for real data it becomes: $$R _x = \begin{bmatrix} r_x & r_x^* & ... & r_x[N-1]^* \\ r_x & r_x & ... & r_x[N-2]^* \\ ... & ... & ... & ... \\ r_x[N-1] & r_x[N-2] & ... & r_x \\ \end{bmatrix} \tag{6}$$

Again this theoretical ACM $$R_X$$ cannot be computed but is estimated from sample function data...

Estimation of the auto-correlation matrix $$R_x$$ can be done in a number of ways (both related) but at the basic level you see that what you should actually estimate is the auto-correlations sequence $$r_x[k]$$ for the lags $$k=0,1,...,N-1$$. Then from the Toeplitz and Hermitian symmetric property of $$R_x$$, you can construct it.

One way to practically estimate ACS $$r_x[k]$$ from the data $$x[n]$$ of length $$L$$ is :

$$\hat{r}[k] = \frac{1}{L} \left( x[k] \star x^*[-k] \right) \tag{7}$$

This yields an estimate of ACS for the lags $$-L < k < L$$, (you only need $$N$$ of those for $$k=0,1,...,N-1$$)

Alternatively, ACM can be estimated as a whole similarly to ACS from:

$$\hat{R}_x = \frac{1}{K} \sum_{k=1}^{K} X(k) \cdot X(k)^H \tag{8}$$ where data vector $$X(k)$$ includes $$N$$ consecutive samples from $$x[n]$$ (specially processed at the edges)

• Now I completed it, I hope I understood your question right... – Fat32 Feb 4 '19 at 18:36
• Regarding the scaling, all estimation equations do have a scale, as you can see too. You are somewhat right that statistical signal processing can provide puzzling equations at first but they all turn to be same later. May be you should try to focus on the very step that causes the first confusion, rather than bundling them together. – Fat32 Feb 5 '19 at 1:06
• Now finally I understood your main reason of confusion and provided an answer for it (put it at the beginning)... – Fat32 Feb 5 '19 at 1:36
• Well... I hate it. Guess I just got it. In the All-Pole Modeling were the autocorrelation matrix is used ($R_x*a_p=-r_x$) you can just factor out 1/N on both sides. That's probably why they don't bother calculating the autocorrelation correctly. I'm freaking out. Why can't people just write it down like it is :) – Mr.Sh4nnon Feb 5 '19 at 15:00
• @Mr.Sh4nnon yes! that's the reason. For computational efficiency purposes, if both sides of the linear equation are weighted by $1/N$ (which is the case for $R_x a_p = r_x$ all-pole ,Yule-Walker eqs) then yo ucan ignore them. However, be very cautious that ACM $R_x$ or ACS $r_x$ when computed alone, do need the weight in for their ergodic estimates from data. – Fat32 Feb 5 '19 at 17:19