Some subtleties in the PSD of digitally generated noise with $S_X(f) = A(f)$

Question

I understand that the title is not the most descriptive, but I didn't know how to describe my problem in a more literal way. The situation is as follows:

I wanted to digitally generate noise with an arbitrary power spectral density $S(f) = A(f)$. To do this, I generated a white noise series $X$ with $S_X(f) = 1$ by simply taking $n$ independent samples from a Gaussian distribution with $\sigma = 1$ and then digitally filter them with filter amplitudes $F(f) = \sqrt{A(f)}$. One can easily show that this leads to $S_Y(f) = A(f)$, where $Y$ is the series resulting from filtering $X$ with $F(f)$.

I have done this in practice using Mathematica's FrequencySamplingFilterKernel and ListConvolve. For this to work one needs to know the desired $A(f)$ of course, which in my example I wanted to be a Lorentzian of center frequency $f_0$ and full-width-at-full-maximum b given by \begin{equation} A(f) = \frac{1}{1+\left(\frac{f-f_0}{b}\right)^2}\end{equation}

So I took this equation, sampled it, created filter amplitudes and filtered the white noise. The resulting PSD $S_Y(f)$ looked as desired when analysed with Welch's Method.

But now I am interested in a subtlety in the above. As has been pointed out the me, the above Lorentzian is not a proper power spectral density of a real signal; $A(f) \neq A(-f)$. But the method still worked, in the sense that I did not get any complex valued data, and that the periodogram showed the correct spectrum at positive frequencies. So my question is, using this method, what happened at the negative frequencies? Did I essentially create a PSD for which \begin{equation} A(f) = \frac{1}{1+\left(\frac{f-f_0}{b}\right)^2} + \frac{1}{1+\left(\frac{f+f_0}{b}\right)^2}\end{equation} or perhaps instead it is more like \begin{equation} A(f) = \frac{1}{1+\left(\frac{f-\textrm{sgn}[f]f_0}{b}\right)^2}\end{equation} where $\mathrm{sgn}[\dots]$ is the sign function.

Those are the only reasonable scenario's I can think of, resulting from the above. They are very similar in behaviour, except around $f=0$, Here the second one has a discontinuous derivative at $f=0$; is that even physically possible?

My question is thus what PSD I ended up with; one of the two I mentioned above or even a third one?

Answer 1

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);

            

---

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