How to simulate a delta-correlated random process I'm trying to do the simulation described in the paper attached, but there is something I don't understand. The author says that the random variables which satisfy the relation (Eq. (4) in the paper) below were used. Could someone tell me how to generate them?
$\left\langle \delta \lambda_0 \right\rangle = 0$
$ \left\langle \delta \lambda_0(z) \delta \lambda_0(z^{'}) \right\rangle = l_c \sigma_{\lambda}^2\ \delta(z-z^{'}) $
$\delta \lambda_0$ : random variable
$l_c$ : correlation length
$\sigma_{\lambda}$ : standard deviation
$z$ : coordinate along the fiber
paper_page1
paper_page2
 A: (I've revised my answer to account for the authors' use of the Dirac delta function, which I had misinterpreted as a Kronecker delta function. This is sometimes called a  delta-correlated random process.)
The authors appear to apply the following simulation procedure (for heuristics, see here):

If $\{X(z)\}_{z\in\mathbb{R}}$ is a continuous Gaussian White Noise process such that
  $$\begin{align}&E[X(z)]=0\\
&E[X(z)X(z')]=c\ \delta(z-z')\end{align}
$$ then a discrete sampling of the process, with a uniform sampling interval $\Delta$, viz.,
  $$X(z),X(z+1\Delta),X(z+2\Delta),...$$ is simulated as an i.i.d  sequence of Gaussian$(\text{mean}=0,\text{variance}=\frac{c}{\Delta})$ random variables, the quality of the simulation increasing as $\Delta\to 0$.  

(The power spectral density of the i.i.d. sequence approaches that of the continuous white noise that it simulates -- flat and equal to the same constant value $c$ -- except that for the i.i.d. sequence it is "band limited", i.e. vanishing outside of a finite-width interval.)    
Thus, the authors take $X=\delta\lambda_0$, $c=\sigma_{\lambda}^2 l_c$, and apparently equate $l_c$ to the sampling interval. E.g., for Figure 1, they took $\sigma_{\lambda}=1\ \mathbb{nm}$ for a fiber of length $500$ meters, and took $l_c=5$ meters, giving a simulation of  $100$ samples. These simulated samples would have been generated as pseudorandom numbers from a Gaussian$(\text{mean}=0, \text{variance}=\frac{c}{\Delta}=\sigma_{\lambda}^2=1\ \mathtt{nm^2})$ distribution. At another point, they take $l_c=50$ meters and remark that this magnitude produced an inaccurate simulation (consistent with the quality depending on a "sufficiently small" sampling interval).

How to generate i.i.d. Gaussian pseudorandom numbers depends partly on the software you use. One option is to use a programming language that has this functionality more-or-less built-in (e.g. in Python, import the random module and use the function gauss(mu,sigma); or, in R, use the function rnorm(n,mu,sigma), etc.). Another option is to write your own generator, using one of various algorithms to transform uniformly distributed pseudorandom numbers, such as discussed here.
