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
(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):
(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
randommodule and use the functiongauss(mu,sigma); or, in R, use the functionrnorm(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.