Motivated by the construction of the It\^o stochastic integral, we consider a step function method to discretize and simulate volatility modulated L\'evy semistationary processes. Moreover, we assess the accuracy of the method with a particular focus on integrating kernels with a singularity at the origin. Using the simulation method, we study the finite sample properties of some recently developed estimators of realized volatility and associated parametric estimators for Brownian semistationary processes. Although the theoretical properties of these estimators have been established under high frequency asymptotics, it turns out that the estimators perform well also in a low frequency setting.