First, the output of the xcorr() function returns lag-0 of the auto-correlation sequence (ACS) estimate at its middle sample, as you recognize. The function argument scaleopt provides normalization of the output: the default value none causes your confusion. If you select biased, instead, then Eq-7 (below) and xcorr() outputs 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 $r_x[-2],r_x[-1],r_x[0],r_x[1],r_x[2]$. But you only need $r_x[0],r_x[1],r_x[2]$. Instead, call it like:
L = 3;
x = [1,2,3];
rxL = xcorr( x, 'biased');
Rx = toeplitz( rxL(L:end) ); % now you're ok!
Finally, distinguish (theoretical) expectations from their (practical) estimations. When algorithms work on data they don't calculate expectations, they estimate. Hence, all expectations are replaced with their estimates in practical applications.
NOTE: Below is some optional discussion.
It seems (to me) that your confusion is about three things.
- Theoretical definition of the ACS $r_{XX}[k]$ of an RP $X[n,s]$.
- Practical estimation of ACS, from a sample sequence $x[n]$ (one realization of the RP) via ergodicity.
- Utilization of a random-vector
Consider the interpretation of an RP $X[n,s]$ as an ordered collection (set) of indexed RVs $X_n(s)$ for each time-index $n$. Each $X_n(s)$, in principle, has its own probability distribution/density function $F_{X_n}(x)$.
Unfortunately, this yields enormously large number and type of possible joint density functions needed to characterize the RP, therefore, things are simplified by assuming that all those RVs have the same (identical) PDF $F_{X_n}(x)$ for all $n$, and that all are independent of each other, and further that first & second moments (mean,variance and aurocorrelations) are independent of partciular time-instants. This simplification leaves an i.i.d. (independent, identically distributed) WSS RP $X[n,s)$ at hand.
One more simplification (notationally) is to replace the notation $X_n(s)$ with $x[n]$ to represent the RP. It may cause some confusions, and I will continue with $X_n(s)$ below.
The (theoretical) ACS of an i.i.d. WSS RP is defined as:
$$ r_x[k] = r_{x x}[k] = r_{X_n X_{n+k}}[n,n+k]= E\{X_n X^*_{n+k}\} = E\{x[n]x^*[n+k] \} \tag{1}$$
Calculation of $r_x[k]$ requires evaluation of expectations: $E\{X_n X_m\} = \int_{-\infty}^{\infty} x_1 x_2 f_{X_n,X_m}(x_1,x_2) dx_1 dx_2 $. This is ok with questions on paper, but for practical applications it's estimated from data, based on ergodicty assumption: theoretical ensemble-averages can be estimated from sample function time-averages. Then you have
$$ \hat{r}_X[k] = \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]$ on RHS is a particular sample function (a realization) of the RP $X[n,s)$. Note also that the summation limits are determined by data window boundaries.
Now consider the random vector $ \bar{X}_N(n) = [X_n , X_{n-1} , ..., X_{n-N+1}]^T $ which includes $N$ consecutive RVs of the RP $X[n,s)$, beginning at index $n$.
From this vector, generate a Random Matrix (RM) via the Hermitian outer product $\bar{X}_N \cdot \bar{X}_N^H$:
$$
RM =
\begin{bmatrix}
X_n X^*_n & X_n X^*_{n-1} & ... & X_n X^*_{n-N+1} \\
X_{n-1} X^*_n & X_{n-1} X^*_{n-1} & ... & X_{n-1} X^*_{n-N+1} \\
... & ... & ... & ... \\
X_{n-N+1} X^*_n & X_{n-N+1} X^*_{n-1} & ... & X_{n-N+1} X^*_{n-N+1} \\
\end{bmatrix} \tag{3}
$$
Taking the expectation of this RM, yields the theoretical auto-correlation matrix (ACM) $R_x$ of the WSS RP $X[n,s)$ :
$$ R_x = R_{XX} = E\{RM\} = E\{ \bar{X_N} \cdot \bar{X_N}^H \} = E\{\begin{bmatrix}
X_n X^*_n & X_n X^*_{n-1} & ... & X_n X^*_{n-N+1} \\
X_{n-1} X^*_n & X_{n-1} X^*_{n-1} & ... & X_{n-1} X^*_{n-N+1} \\
... & ... & ... & ... \\
X_{n-N+1} X^*_n & X_{n-N+1} X^*_{n-1} & ... & X_{n-N+1} X^*_{n-N+1} \\
\end{bmatrix} \} \tag{4}$$
and this results in:
$$ R_x = E\{ \bar{X}_N \cdot \bar{X}_N^H \} = \begin{bmatrix}
r_x[0] & r_x[-1] & ... & r_x[-N+1] \\
r_x[1] & r_x[0] & ... & r_x[-N+2] \\
... & ... & ... & ... \\
r_x[N-1] & r_x[N-2] & ... & r_x[0] \\
\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[0] & r_x^*[1] & ... & r_x^*[N-1] \\
r_x[1] & r_x[0] & ... & r_x^*[N-2] \\
... & ... & ... & ... \\
r_x[N-1] & r_x[N-2] & ... & r_x[0] \\
\end{bmatrix} \tag{6}
$$
Again, this theoretical ACM $R_x$ is estimated from sample observation data.
Estimation of the ACM $R_x$ can be done in a number of related ways, but basically, it reduces into an estimation of the ACS $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 similarly to ACS as:
$$ \hat{R}_X = \frac{1}{L} \sum_{k=0}^{L-1} \bar{X}_N(k) \cdot \bar{X}_N(k)^H \tag{8}$$ where the data vector $X_N(k)$ includes $N$ consecutive samples of $x[n]$ beginning from $x[k]$, and $L$ is the length of the observation sequence $x[n]$.
