I'm having some trouble with the function Periodogram, mainly with getting it to produce correct amplitudes for the spectrum it returns. My question is how to fix this in order to get the correct amplitudes. I understand that this is closely related to signal processing, but as the trouble I'm having comes from Mathematica's function I thought I'd ask here.
The specific example is a little bit involved, but I suppose I will have to introduce it to get to the core of my problem. What I am trying to do is simulate noise that follows the Ornstein Uhlenbeck process; this is noise with a spectral density given by $S_{xx}(f) = \frac{2c\tau^2}{1+(2 \pi f)^2}$. Now, there are a few ways to do this; I'm personally following Rice, 1944. I won't go into the theory too much, but the idea is that you use a sum of cosines with random phases and amplitudes given by gaussians with zero mean and standard deviation given by the square root of the spectral density.
The code in which I do this is given by
Γ = 10; (* bandwidth of Lorentzian in MHz*)
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 *)
fmax = 10;
Δf = fmax/(K - 1);(* freq. step *)
fk = Range[0, fmax, Δf]; (* sampling frequencies *)
ck = Sqrt[2 S[fk] Δf]; (* fourier coefficients *)
ϕ = RandomReal[{0, 2 π}, K]; (* iid phases *)
X[t_] := Total[ck Cos[2 Pi fk t - ϕ]] ;(* fourier rep. of noise *)
times = Table[i, {i, 0, 1/Δf, 1/(2 fmax)}];
Noise = Table[X[times[[i]]], {i, 1, Length[times]}];
Periodogram[Noise,
PlotRange -> All,
SampleRate -> 2 fmax,
ScalingFunctions -> "Absolute",
Joined -> True]
The first few lines just define some parameters for the spectral distribution (which are rather arbitrary at this point) and after that I built up the coefficients for the cosine amplitudes. I then generate some time values with a sampling rate given by the Nyquist rate ($\Delta t = 1/(2f_{max})$) and evaluate the noise.
If I plot the noise, it looks rather good. If I however plot the periodogram, things get a little shaky. While the general trend is where it should be, the amplitudes are not; they are off by say, a factor of ten, when compared to the analytic expression. I'm trying to figure out how to fix this.
I've been going over the documentation of periodogram but I'm not entirely sure why this is happening. Now I understand there are different types of periodograms. You can use Barletts method, Welch's method, and so on, by averaging over certain lengths of the data that do or do not overlap, and using smoothing functions. While I definitely admit that I do not understand how one optimally chooses whether or not to have data overlap and how big each section/overlap should be (for example the documentation gives Periodogram[data, 64, 32, HammingWindow, PlotRange -> All] for Welch's method, but the 64 and 32 feel arbitrary to me and I wouldn't know how to pick them myself) experimenting with this does not seem to solve the problem of the amplitude.
To clarify the above, this illustration shows the analytic spectral density versus the numerically calculated one
