🚀 go-pugleaf

The easiest way to approximate this is by taking white noise through a number
of evenly-spaced low-pass filters and adding the results, weighted according
to the frequency distribution you wish.  You can use the same noise source for
each filter, or a different one for each.  The analysis is slightly different
but both work out in the end.

The response of the filters isn't important, as long as they are identical,
with frequency scaled of course.  Simple filters will give less ripple, so a
simple RC is what you want.  That's a lossy integrator for you digital folks.
Sharper is not better in this application.

The appoximation is good because the errors cancel each other very well.  The
ripple is 0.1% even when the filters are a decade apart!  Clearly fringing
effects will be dominant here.  The filters should cover the frequency band of
interest, plus an additional one or two on each end to keep errors reasonable
at the ends.

This method works for any exponent, although there are some constraints on the
filters.  Specifically at any given frequency the output from each filter,
including weighting, should go to zero as the filter frequency goes higher or
lower.  For example for f^-3 power you need at least a second-order low-pass
filter.
--
Joe Keane, amateur mathematician