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

First, the output of xcorr() function will return the lag-0 of the ACS $$r_x[k]$$ of $$x[n]$$, at the middle sample as you recognize. Then, look at the function argument scaleopt which provides normalization to the computed output. The default being none causes the 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 full output of xcorr as the input to the toeplitz function, since the xcorr() output will include $$r_x[k]$$ estimates for all lags 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 ;-)

Finally, distingusih a theoretical expectation from its practical estimation. When algorithms work on data, they estimate; they don't calculate expectations, which is a theoretical tool. Hence, all expectations are replaced with their estimates when algorithms are run in 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, it's assumed 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}$$

Calculation of this theoretical ACS $$r_x[k]$$ requires you to evaluate expectations of the kind: $$E\{X\} = \int_{-\infty}^{\infty} x f_X(x) dx$$ . This is typically done with questions on the paper. But in practical applications you perform statistical calculations, you estimate ACS from available data, based on ergodicty assumption: theoretical ensemble-averages can be estimated from 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 window 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}$$

where $$E\{X_n X^*_m\}= r_x[n-m]$$, and becuse of Hermitian symmetry of Auto-Correlations, 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)

• That was a very good answer. However, I'm still confused about Matlabs output. If I use your second equation for $\hat{r}_X[k]$ for lag 0 (the upper left most entry) then I would not get $1+2^2+3^2=14$. From intuition and my daily use of convolution I know thats true for lag 0. But I can't connect it with the newest theory I'm learning about autocorrelation matrices. All those formulas do not result in what I used to know. I still make mistakes regarding ensembles and means. But again, with all I've read above, this doesn't fit. There should at least be a 1/N or (1/K) factor for the estimation Feb 4 '19 at 23:39
• 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. Feb 5 '19 at 1:06
• Yes, thats perfectly what I was confused about. I'm sorry I didn't explained it well the first time. It's actually pretty hard to say what you don't understand yet. On page 393 spectrum estimation Hayes states the same formula to estimate $\hat{r}_x(k)$ as you. However, in our exam cheat sheet ACS is calculated without 1/N. This is either a mistake or one of those damn conventions. I hope this gets better with time. At the moment I'm freaking out every time a formula doesn't look exactly like I'm used to. Feb 5 '19 at 14:56
• 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 :) 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. Feb 5 '19 at 17:19