Random processes can be fully characterized within a probabilistic framework at each time instant by their probability density function or, equivalently, by the characteristic function. Although statistical moments of integer order are related to the coefficients of the Taylor expansion of the characteristic function, they are insufficient for a complete description of random processes. This limitation arises because high-order moments may diverge, and the Taylor expansion of the characteristic function converges only near zero, diverging at points farther from zero. Recently, complex fractional moments, corresponding to the Mellin transform of the probability density function (or of the characteristic function), have been introduced. Unlike the Taylor series, the inverse Mellin transform recovers the probability density function (or the characteristic function) as a Fourier series on a logarithmic scale. A key advantage of this approach is that the real part of the domain remains constant in the inverse Mellin transform, effectively addressing divergence issues. Recent research has shown that the loss of Markovianity in fractional differential equations subjected to normal white noise can be mitigated if the response process is self-similar. This property significantly reduces the computational effort required to compute the statistics of the response process. However, to date, self-similarity has only been defined for integer-order statistics. In this paper, the concept of self-similarity is extended from classical statistical moments to complex fractional moments, without loss of generality, for the case of fractional Brownian motion. Analytical solutions of both the complex fractional moments of the fractional Brownian motion’s probability density function and the complex fractional moments of the fractional Brownian motion’s characteristic function have been calculated and an innovative approach based on the self-similarity of complex fractional moments and on the aforementioned analytical solutions has been proposed in order to quickly reconstruct the probability density function and the characteristic function from numerical data. Validation through Monte Carlo simulation is presented and the results are discussed in detail showing that the proposed approach greatly reduces the computational burden required to calculate the complex fractional moments.