Having trouble working understanding a noise shaping filter

As context, I'm trying to understand and implement this Von Karman turbulence model in a simulation I'm working on. In order to generate turbulence velocities on the fly using the model, a factorization of an empirically determined power spectrum density is given which supposedly serves as a shaping filter which will produce colored noise with the right PSD when driven with white noise.

I'm trying to understand how this works with a simple example so I know how to implement it in python using scipy. I posit a simple PSD given by $$\Phi(s) = \frac{1}{1+s^2}$$ and factor it into the filter $$G(s) = \frac{1}{1+s}$$ such that $$\Phi(\omega) = G(i\omega)G(-i\omega)$$.

I drafted a python script to run this filter on a set of random Gaussian samples.

from scipy.signal import welch, lfilter, bilinear, filtfilt
from numpy.random import normal
from numpy import linspace
import matplotlib.pyplot as plt

w = linspace(0, 5, 1000)
Sxx = 1/(1 + w**2)
plt.plot(w, Sxx)

z, p = bilinear(, [1, 1])
x = normal(size=1000000)
y = filtfilt(z, p, x)
f, Pxx = welch(y)
plt.plot(f, Pxx)
plt.show()

However, the output PSD and the analytical PSD do not match up, and I'm struggling to see what I'm missing here. I suspect its something to do with the discretization that I'm not understanding.

Apart from that, you get frequency warping from the bilinear transform. Increase the sampling frequency (which probably defaults to $$1$$) to alleviate that problem.