# Quantifying how good I'm producing a spectral density

I recently asked a question here on how to create noise with a specific spectral density $S_{xx}(f)$. The helpful people of this stack exchange told me one way to do it was filtering white noise with a unitary spectral density using a filter with impulse response $H(f) = \sqrt{S_{xx}(f)}$ and then suggested how to do this. The way that was proposed was using a frequency sampling based FIR filter. As I use Mathematica the idea was to use FrequencySamplingFilterKernel.

Now, this seems to work reasonably well qualitatively, but I am trying to find a way of quantifying this. So my question comes down to how I quantitaviely test how well the noise produced with my filter resembles the desired spectral density.

To take a specific example, I am trying to produce (not exclusively, but to begin with) the Ornstein Uhlenbeck spectral density given by $S_{xx}(f) = \frac{2c\tau^2}{1+(2 \pi f)^2}$. So I suppose it is the easiest if I just show my work and then get to my specific question. What I did was generate the amplitudes

Γ = 10; (* bandwidth of Lorentzian*)
P = 2; (* std dev of noise *)
τn = 1/Γ;(* relaxation time of OU process *)
c = 2 P^2 Γ;(* diffusion constant of OU process *)
S[f_] := (2 c τn^2)/(
1 + (2 Pi f  τn)^2); (* power spectral density *)
K = 100; (* number of frequency components taken into account *)
fmax = 5; (*max frequency*)
Δf = fmax/(K - 1);(* freq step *)
fk = Range[0, fmax, Δf]; (* sampling frequencies *)
FilterAmps = Sqrt[S[fk]]; (* Amplitudes for the FIR filter *)


The parameter choices are arbitrary, and definitely not necessarily what I want to do in practice, but that shouldn't matter too much for now. I then generate some white noise (the times list is not needed at this point, but I'm just doing it to keep track of what I'm doing for now):

Ns = 1000; (* number of time points generated; output time points \
after filter is n = Ns - 2K + 2*)
Δt = 1/(2 fmax); (*Nyquist rate*)
Tmax = (Ns - 1) Δt;
times = Range[0,
Tmax, Δt]; (*Using Nyquist rate as step size*)
UnitDataWhite =
RandomVariate[NormalDistribution[0, 1], Length[times]];


Okay, great, so then we filter

FilteredNoise =
ListConvolve[FrequencySamplingFilterKernel[FilterAmps],
UnitDataWhite];
FilteredTimes =
Range[0, (Length[FilteredNoise] -
1) Δt, Δt];


Giving me a dataset of $n = N - K + 2$ noise datapoints which should behave according to the spectral density of the Ornstein Uhlenbeck process. The problem is, how do I quantify how well this works? The first thing that came to mind was just constructing the power spectral density and seeing how well it agreed with theory. Well, there I ran into my first real meeting with the fact that one does not simply find the power spectral density from a dataset to some arbitrary degree of precision. Mathematica has a built in function periodogram which is supposed to do this, but while the shape is definitely as desired, is not exactly.

Okay, so then I thought perhaps it is too much to ask to reconstruct the spectral density in this way. Should I look instead at other quantities that carry the same information, such as the autocovariance? For the OU process this is given by $C_{x}(t) = \frac{c\tau}{2} e^{-t/\tau}$ so that's also relatively nice. But calculating the autocovariance with CovarianceFunction again runs into amplitude trouble, so that doesn't seem like a good method either.

So in the end my question boils down to this. Creating noise with a spectral density using the way I described here, how do I reliably quantify how well I am succeeding?

*Note that I know that there are nicer ways of creating OU noise, but I am trying to construct a method for general spectral densities

• You are getting the correct shape even by using FilterAmps = S[fk] instead of FilterAmps = Sqrt[S[fk]]? – SleuthEye Mar 4 '16 at 3:23
• That's an embarrassing mistake. The question still holds, but that does already fix it a little. – user129412 Mar 4 '16 at 14:40
• I can probably add to this in the following way. Although perhaps with carefully picking the specific parameters of periodogram (averaging over partitions, averaging, smoothing functions) I can eventually get close to the correct power spectral density, but this is only because I know where I want to go. I haven't really been able to find a reliable guide for how to pick your overlap/averaging sizes depending on the data size, most that I can find tells you that you should adjust this according to your specific situation, which doesn't say that much. – user129412 Mar 4 '16 at 15:04
• Keep in mind that a periodogram is only an estimation method, and you are generating colored noise which is by nature unknown (you could generate something that has exactly the spectrum you want in a deterministic way, but then it would no longer be 'noise'). If that periodogram converges to what you expect as you increase the number of samples (so you can get a better estimate) then your signal generation is probably fine. – SleuthEye Mar 4 '16 at 15:16
• That does make sense, but I'm aiming to use this for an experiment where it will be quite crucial to accurately control the parameters of the noise. In the end I'll be moving on to a physical spectrum analyser so perhaps that will help. Still a little confused on why the autocorrelations are misbehaving however, but perhaps it is simply some vertical offset that I need to understand. – user129412 Mar 4 '16 at 15:20