Wavelet analysis has turned out to be a quite useful tool for estimating Hölder exponents as well as spectral densities of fractal processes, typically fractional Brownian motion (fBm) and related processes. In real life applications, the process is observed at discrete times and its (theoretical) wavelet coefficients have therefore to be replaced by their discretizations. This paper studies the discretization error of a wavelet coefficient of a stationary increments fractal Gaussian process belonging to a class which includes fBm and multiscale fBm. We obtain a harmonizable representation formula for the latter error and then optimal lower and upper bounds of its mean square. Our results remain valid under very mild assumption on the analyzing wavelet ψ: we need no vanishing moment condition (which is unusual in the wavelet setting), we only have to assume that the function ψ be compactly supported and Lipschitz continuous.