Some subtleties in the PSD of digitally generated noise with $S_X(f) = A(f)$
You are only looking at the positive frequencies. The negative frequencies are also there. So yes, the spectrum is indeed (with some abuse of notation)

\begin{equation}
A(f) \approx \frac{1}{1+\left(\frac{f-\mathrm{sgn}(f)f_0}{b}\right)^2}\end{equation}

such that the PSD is mirrored with respect to $y$-axis. 

I don't speak Mathematica, but I can tell you how it works. Your FrequencySamplingFilterKernel is in fact a simple frequency sampling approach to build FIR filters with arbitrary-magnitude frequency response (similar to Matlab's fdesign.arbmag). Given you have $N$ points of the frequency response, it essentially calculates an impulse response (by taking the inverse  DFT) whose Fourier transform is equal to your desired spectrum when evaluated at the corresponding frequencies. That is:

$$h[n]=\frac{1}{N}\sum_{k=0}^{N-1}H(k)e^{j\frac{2\pi}{N}k}$$

where $H(k)$ is the desired DFT (That is, $N$ samples from $\sqrt{A(f)}$).

Notice that the default impulse response $h[n]$ in the documentation of FrequencySamplingFilterKernel (under details and options) is an odd-length FIR filter with even symmetry. From this point everything is just regular FIR filtering with $h[n]$, and when a real process is filtered it will generate a PSD with even symmetry:
$$S_Y(f)=|H(f)|^2S_X(f)$$
where $|H(f)|^2$ is equal to $|A(f)|$ at the given $N$ desired points. 



If you want to do it in Matlab:

N = 100;
f = 0:0.01:1;
b = 1;
f0 = 0.75;
A =sqrt(1./(1+((f-f0)./b).^2));
d = fdesign.arbmag('N,F,A',N,f,A);
hd = design(d,'freqsamp');
w = -pi:.001:pi;
freqz(hd,w);
x=randn(1000000,1);
y=filter(hd,x);
psd(y);




&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;



with low-order FIR approximation (N=8):



&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;