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

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?

• I don't think this can be answered without a careful look at your code. – MBaz Oct 18 '16 at 21:35
• @MBaz Hm. Really? I somehow thought it would be implicit in using filter amplitudes sampled from a PSD that does not include any overlap between the positive or negative frequencies, or in other words, that is just cut off at $f=0$. But if it cannot be said from these considerations, would it be worthwhile to post the actual code or is this not something that is done on this site? – user129412 Oct 18 '16 at 22:51
• Well, the thing is that $S_Y(f)=S_X(f)|A(f)|$ is not symmetric, because $Y$ is complex. If you get a real $Y$, or a symmetric $S_Y(f)$, something is wrong. It's hard to speculate about what happened without looking at what you did exactly. If the code is long, consider posting it somewhere else (like a github gist) and providing a link. Somebody fluent in Mathematica could then provide more help. As an alternative, provide a more detailed account of what you did. – MBaz Oct 19 '16 at 0:58

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): • not sure i agree with you about the "simple frequency sampling approach to build FIR filters with arbitrary frequency response." the OP has "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)}$"  the question is: how, precisely, did the OP digitally filter the samples? – robert bristow-johnson Oct 19 '16 at 5:02
• doesn't answer my question. how does $F(f)$ get turned into an impulse response? – robert bristow-johnson Oct 19 '16 at 7:28
• Interesting response! That does seem to answer my question. And you have, much more eloquently than I can, indeed explained my procedure. One question that remains is about this discontinuity at $f = 0$, or more explicitly, the discontinuity in the derivative. Is that not problematic? – user129412 Oct 19 '16 at 10:47
• @user129412 Although it looks very much to be, it is not a discontinuity. You will get a very smooth edge at $f=0$ with a low filter order (N=8 for instance). So it is only the high filter order that makes it sharp. The effect looks similar to when you truncate the Fourier series of a function with discontinuity. The higher the number of terms you consider, the sharper the transition becomes. I will also add low-order FIR result to the answer. – msm Oct 19 '16 at 12:44