An fBm is characterized by a power spectrum $P(f) = Cf^{-(2H + 1)}$ with $0 < H < 1$ being the Hurst parameter. Why can't I just take the square root of the power spectrum $P(f) = Cf^{-\alpha}$, multiply with $e^{i\theta_n}$ ($\theta_n$ being $N/2$ random phases in $[0, 2\pi]$), and then do an inverse Fourier transform? This can produce arbitrarily long fBm's with a $1/f^\alpha$ spectrum.
A popular book on chaos theory by A.A. Tsonis (pp. 62-53) suggests the exact same method for generating an fBm trace:
How To Generate an fBM Trace with a Desired Spectral Density Function. A time series $X_i$, $i = 1, ..., N$, whose spectral density function satisfies the relation $S(f) = Cf^{-a}$ where $a = 2H + 1$, is obtained via the relation
$$x_i = \sum^{N/2}_{k=1} \left[Ck^{-a}\left(\frac{2\pi}{N}\right)^{1-a}\right]^{1/2}\cos{\left(\frac{2\pi ik}{N} + \phi_k\right)}$$
where $C$ is a constant and $\phi_k$ are $N/2$ random phases randomly distributed in $[0, 2\pi]$.
While this method is very simple, almost no papers (for e.g. this review or the fArma R package) on fBm generation discuss it. Am I missing something?
