I have a random process which is defined according to the model, \begin{equation*} Y_i=Bf_i+W_i, \quad i=1,\ldots,N \end{equation*}
Where $B,W_1,\ldots,W_N$ are i.i.d. with, \begin{gather*} W_i\thicksim\mathcal{N}(0,\sigma_w^2)\\ B\thicksim\mathcal{N}(0,\sigma_b^2)\\ \end{gather*}
The covariance matrix was determined to be, \begin{equation*} \mathbf{C}= \begin{bmatrix} f_1^2\sigma_b^2+\sigma_w^2&f_1f_2\sigma_b^2&\cdots&f_1f_N\sigma_b^2\\ f_1f_2\sigma_b^2&f_2^2\sigma_b^2+\sigma_w^2\\ \vdots& &\ddots\\ f_1f_N\sigma_b^2&\cdots& &f_N^2\sigma_b^2+\sigma_w^2 \end{bmatrix} \end{equation*}
And therefore, \begin{equation*} \mathbf{y}\thicksim\mathcal{N}\left(\mathbf{0}_N,\mathbf{f}\,\mathbf{f}^T\sigma_b^2+\mathbf{I}\,\sigma_{w}^2\right) \end{equation*}
I am now trying to show that the derivation is correct in Matlab. I created $R=500$ realizations of a random signal. Each realization took a value from $B$. Each random signal was $N=100$ items long and the $i_{th}$ item within the random signal took a value from $W_i$. For simplicity, I assumed that the vector $\mathbf{f}$ is a vector of ones.
If I autocorrelate each random signal and take the average, I get, ${\mathbf{f}\,\mathbf{f}^T\sigma_b^2+\sigma_{w}^2}$. If I autocorrelate the values for a particular time-slice, $i$, across all realizations, I also get, ${\mathbf{f}\,\mathbf{f}^T\sigma_b^2+\sigma_{w}^2}$.
I am struggling however to simulate the covaraince terms. For example, I thought that if I cross-correlate realization 1 with realization 2, then I would be able to get ${\mathbf{f}\,\mathbf{f}^T\sigma_b^2}$, but no success.