Simulating atmospheric turbulence This problem comes from simulation of atmospheric turbulence. 
It starts with the following transfer function in the frequency domain:
\begin{equation}
G_v(s)=\sigma k_1 \frac{1 + k_2 s}{\left(1+k_3 s\right)^2}
\end{equation}
This comes from the so-called Dryden turbulence model, which says that the speed of the turbulence perturbation on the spatial dimension $v$ is $G_v(s)$ applied to band-limited white noise. As usual in block-based system simulation, there are black-box assumptions which I cannot blindly accept and use, so I am trying to replicate the model with my own code. 
Still in the frequency domain, I completed the equation with the meaning of the model, that is
\begin{equation}
y_v(s)=G_v(s)sW(s)
\end{equation}
where $y_v(s)$ is the turbulence speed and $sW(s)$ is the derivative of a Wiener process in the frequency domain, which should be equivalent to a white noise. Treating $W(t)$ as a generalized random process allows the existence of its $n$-th order derivatives and, by consequence, their equivalents in the frequency domain. Substituting the formula of $G_v(s)$ from above in the latter, we obtain
\begin{equation}
y_v(s)\left(1+k_3 s\right)^2=\sigma k_1 \left(1 + k_2 s\right)sW(s)
\end{equation}
Doing the inverse Laplace transformation to the time domain and abstracting the constants, we obtain a stochastic differential equation of the following form:
\begin{equation}
a_2 y_v^{\prime\prime}(t)+a_1 y_v^{\prime}(t)+a_0 y_v(t)+b_2W^{\prime\prime}(t)+b_1W^{\prime}(t)=0
\end{equation}
where $W(t)$ is the Wiener process in the time domain and, by consequence, $W^{\prime}(t)$ should be white noise and $W^{\prime\prime}(t)$ should be the derivative of white noise.
Now, to simulate this equation numerically, it would be relatively easy if not for the $W^{\prime\prime}(t)$. With a We may interpret $W^{\prime}(t)$ either in the straightforward sense, 
\begin{equation}
\frac{dW}{dt} \approx \frac{1}{h} \Big[ W(t+h) - W(t) \Big]
\end{equation} 
and we can recur to Euler-Maruyama to solve this numerically. However, how to solve numerically for $W^{\prime\prime}(t)$?  
 A: Am I correct that you are asking how one could interpret a second-order derivative of a white noise?
You really need to consult a mathematician. I have an idea how you can define a second order derivative of white noise, but rely on it only if you are desperate, I don't know how correct this is.
First, recall how you view the first derivative through a finite difference:
$$
  \frac{dW}{dt} = \frac{W(t+\Delta t) - W(t)}{\Delta t} = \frac{N_t(0,1) \sqrt{\Delta t}}{\Delta t} = N_t(0,1) \frac{1}{\sqrt{\Delta t}},
$$
where $N_t(0,1)$ is a standardized Gaussian random variable with zero mean and variance of $1$. The subscript $t$ is just to remind that at each time $t$ you calculate a different random variable.
Now we can go by analogy and use a forward finite difference approximation for the second derivative:
\begin{align}
  \frac{d^2W}{dt^2}
  &= \frac{W(t+2\Delta t) - 2W(t+\Delta t) + W(t)}{\Delta t^2}\\[1ex]
  &= \frac{\bigl[W(t+2\Delta t) - W(t+\Delta t)\bigr] - \bigl[W(t+\Delta t) - W(t)\bigr]}{\Delta t^2}\\[1ex]   &= \frac{N_{t+\Delta t}(0,1)\sqrt{\Delta t} - N_{t}(0,1)\sqrt{\Delta t}}{\Delta t^2} \\[1ex]
  &= N_t(0,2) \frac{1}{\Delta t^{3/2}} \\[1ex]
  &= N_t(0,1) \frac{\sqrt 2}{\Delta t^{3/2}}.
\end{align}
Here I used that the sum of two standard Gaussian random variables is again Gaussian with variance $2$ (standard deviation $\sqrt 2$).
I am using forward differences due to the use of Ito's calculus.
Now you can take you differential equation and turn it into a finite difference scheme.
