# 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, 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.

The function filtfilt filters the signal twice (forward and backward) in order to eliminate phase distortions. As a side effect, your signal is not filtered by the transfer function corresponding to the numerator and denominator polynomial coefficients supplied to the routine, but by its squared magnitude. You should use an ordinary filter routine, such as lfilter.
Apart from that, you get frequency warping from the bilinear transform. Increase the sampling frequency (which probably defaults to $$1$$) to alleviate that problem.