I want to numerically study the behavior of a system described by a set of differential equations in the presence of colored noise. It seems that the standard procedure is to use the Langevin equation: $ y=fdt+gdX; dX=-aXdt+bdW$, where $W$ is the Wiener process.

The spectral density function of $X$ is $1/f^2$ type. What about the spectral density function of $dX$? If my MATLAB calculations are correct, in this case the power density increases with the frequency. If this is true, how can I introduce $1/f$ noise in my system?