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: $ dy=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?