We present a novel method for stochastic interpolation of sparsely sampled time signals based on a superstatistical random process generated from a multivariate Gaussian scale mixture. In comparison to other stochastic interpolation methods such as Gaussian process regression, our method possesses strong multifractal properties and is thus applicable to a broad range of real-world time series, e.g.~from solar wind or atmospheric turbulence. Furthermore, we provide a sampling algorithm in terms of a mixing procedure that consists of generating a~$1+1$-dimensional field~$u(t,\xi)$, where each Gaussian component~$u_\xi(t)$ is synthesized with identical underlying noise but different covariance function~$C_\xi(t,s)$ parameterized by a log-normally distributed parameter~$\xi$. Due to the Gaussianity of each component~$u_\xi(t)$, we can exploit standard sampling alogrithms such as Fourier or wavelet methods and, most importantly, methods to constrain the process on the sparse measurement points. The scale mixture~$u(t)$ is then initialized by assigning each point in time~$t$ a~$\xi(t)$ and therefore a specific value from~$u(t,\xi)$, where the time-dependent parameter~$\xi(t)$ follows a log-normal process with a large correlation time scale compared to the correlation time of~$u(t,\xi)$. We juxtapose Fourier and wavelet methods and show that a multiwavelet-based hierarchical approximation of the interpolating paths, which produce a sparse covariance structure, provide an adequate method to locally interpolate large and sparse datasets.