This is just an application of Ito's formula. Normally, Brownian motion $W(t)$ is only defined for nonnegative times $t$, so I'm assuming that $W(t)=0$ for $t\le0$, although it doesn't really matter much to the answer.
For a fixed $t$, set $V_s=F(t-s)$, which is a deterministic finite variation process. As $V$ is of finite variation, the quadratic covariation $[V,W]$ vanishes. So, by Ito's formula,
$$
\int_{-\infty}^TF(t-s)\,dW(s)=\int_{-\infty}^TV_s\,dW(s)
=V_TW(T)-\int_{-\infty}^TW(s)dV_s.
$$
Using the change of variables $u=t-s$ in the final (Stieltjes) integral,
$$
\int_0^TF(t-s)\,dW(s)=F(t-T)W(T)+\int_{t-T}^\infty W(t-u)dF(u)
$$
Using Ito's isometry, the integral on the left converges as $T\to\infty$ iff $\int_{-\infty}^tF(s)^2\,ds$ is finite for all $t$ (convergence is in $L^2$, although it can be shown to be in $L^p$ for all $1\le p <\infty$ and almost sure). The first term on the right vanishes in $L^2$ as $T\to\infty$ iff $F(-t)^2t\to0$ as $t\to\infty$. So, I'm assuming that these conditions hold. Taking the limit in $L^2$ as $T\to\infty$ gives $F\otimes W=W\otimes F$.
To compute this, $F\otimes W$ is just a standard stochastic integral, and a Riemann sum approximation method should work fine. By Ito's Isometry, $F\otimes W(t)$ is normally distributed with mean $0$ and variance $\int_0^t F(s)^2\,ds$. If you want to get faster convergence in distribution for $F\otimes W$, you could do something like use an improved approximation to the integral for the variance (e.g., Simpson's rule) and then use the corresponding weights to adjust the Riemann sum approximation to obtain the correct variance.
