While it doesn't exactly give the power spectral density formulated by you, a bimodal (when considering also negative frequencies) power spectral density is given by the time evolution of the position or displacement of a harmonic oscillator driven by Brownian motion. The only modification compared to the Ornstein–Uhlenbeck process is inclusion of inertia, which enables oscillations when the harmonic oscillator is underdamped. An Ornstein–Uhlenbeck process is governed by the differential equation:
$$\dot x(t) = -\theta x(t) + F(t) \tag{1}$$
The dot above the $x(t)$ indicates the derivative of $x(t)$. In an infinitesimal time step $dt$, displacement $x(t)$ relaxes a step $\theta\,dt\,x(t)$ towards zero and is also affected by white noise external forcing $F(t)$. There is no inertia: neither the relaxation term or the external forcing term describe "memory of past velocity". I have poor differential equation handling skills, but my hunch is that if a bit of the past velocity $\dot x(t-dt)$ was mixed into the current velocity, it would give this modified form with inertia:
$$\ddot x(t) = -\gamma_m \dot x(t) - f_m^2 x(t)+F(t) \tag{2}$$
Here $\ddot x(t)$ denotes the second derivative of $x(t)$. Now there are terms both for damping of velocity and for relaxation of the displacement. For the modified process, the power spectral density is of form (Gröblacher 2012, Coffey & Kalmykov 2012):
$$S(f) = \frac{a}{\left(f_m^2-f^2\right)^2 + f^2\gamma_m^2}, \tag{3}$$
where $\gamma_m$ is the width of the resonant peak at half maximum, and $a$ is a normalization constant which has also absorbed the power spectral density of the driving force which is white (constant in frequency) for a Brownian heat bath. The true peak frequency $f_\text{res}$ is (Gröblacher 2012):
$$f_\text{res}=f_m\sqrt{1-\frac{\gamma_m^2}{2f_m^2}}.$$
Figure 1. Power spectral density of the displacement of a massive (as opposed to massless) harmonic oscillator in a heat bath, with $a = 1$, $f_m = 1$ and $\gamma_m=\frac{1}{2}$.
$S(f)$ of Eq. 3 is identical in form to the square of the magnitude frequency response of a two-pole resistor-inductor-capacitor (RLC) filter:
$$|H(if)|^2 = \frac{a}{(f^2 - f_0^2)^2 + f^2\,(f_0/Q)^2}$$
If the filter is fed white noise, the power spectrum of the output will be identical to $S(f)$. It is characterized by a -12 dB/octave (exactly -40 dB/decade) asymptotic slope:
Figure 2. Square of the magnitude frequency response of the two-pole filter, with logarithmic axes. The constants are equivalent to those in Fig. 1.
I'll return to the "just add inertia" argument once more. Consider a discretization of the Ornstein–Uhlenbeck process defined in Eq. 1:
$$\dot x[k] = c_0 x[k] + F[k]$$
We can define velocity and acceleration as:
$$\dot x[k] = x[k] - x[k-1]$$
$$\ddot x[k] = \dot x[k] - \dot x[k-1] = x[k] - 2 x[k-1] + x[k-2]$$
If we add damped inertia to the process as a term $c_1 \dot x[k-1]$:
$$\begin{array}{rrcl}
&\dot x[k] &=& c_0 x[k] + c_1 \dot x[k-1] + F[k]\\
\Rightarrow& x[k] - x[k-1] &=& c_0 x[k] + c_1 x[k-1] - c_1 x[k-2] + F[k]\\
\Rightarrow& 0 &=& (c_0 - 1) x[k] + (1 + c_1) x[k-1] - c_1 x[k-2] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 + c_0 - 1) x[k] + (1 - c_1) x[k-1] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 + c_0 - 1) x[k] - (c_1 - 1) x[k-1] + F[k]\\
\Rightarrow& c_1 \ddot x[k] &=& (c_1 - 1) \dot x[k] + c_0 x[k] + F[k]\\
\Rightarrow& \ddot x[k] &=& \frac{c_1 - 1}{c_1} \dot x[k] + \frac{c_0}{c_1} x[k] + \frac{1}{c_1}F[k]\\
\end{array}$$
We have arrived at the discrete equivalent of Eq. 2.
References
Gröblacer, Simon (2012) "Quantum Opto-Mechanics with Micromirrors: Combining Nano-Mechanics with Quantum Optics" DOI: 10.1007/978-3-642-34955-3, chapter 2, page 7.
Coffey, William T.; Kalmykov, Yuri P. (2012) "The Langevin Equation: With Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering" 3rd ed. ISBN: 978-981-4355-66-8, page 241.