4
$\begingroup$

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 enter image description here

$\endgroup$
1
$\begingroup$

The Periodogram is the squared magnitude of the discrete Fourier transform, thus you can get an equivalent answer using

ListPlot[Take[Transpose[{2 times,Abs[Fourier[Noise]]^2}], {1, (Length[Noise] + 1)/2}], Joined -> True]

The amplitudes are different due to the normalisation. If you multiply by 1/Pi then everything seems to be alright.

ListPlot[Take[Transpose[{2 times,Abs[(1/Pi)Fourier[Noise]]^2}], {1, (Length[Noise] + 1)/2}], Joined -> True]

In order to achieve this using the Periodogram function, simply scale Noise with 1/Pi i.e.

Periodogram[(1/Pi) Noise, PlotRange -> All, SampleRate -> 2 fmax, ScalingFunctions -> "Absolute", Joined -> True]

In order to tune further the normalisation play with the FourierParameters option.

$\endgroup$
  • $\begingroup$ Hmm, that is interesting. You're saying that there is a factor of pi missing from the normalization when you use the default Fourier settings? I see that it works in this specific example, but it doesn't make too much sense to me. Without the factor of Pi unit gaussian data (mean zero, sd 1) gives me a nice spectral density of around one. $\endgroup$ – user129412 Mar 4 '16 at 15:07
  • $\begingroup$ To clarify what I mean, take for example UnitDataWhite = RandomVariate[NormalDistribution[0, 1], 1000]; Periodogram[UnitDataWhite, 64, 32, ScalingFunctions -> "Absolute", Joined -> True] $\endgroup$ – user129412 Mar 4 '16 at 15:18
  • $\begingroup$ Nothing is missing. The default definition of the discrete Fourier transform does not have a Pi term. But the FourierTransform (where the formula S_{xx}, that you provide, comes from) has a term 1/Sqrt(2Pi) in front of it which after squaring yields the 1 over Pi. Check the details for the normalisation on Fourier and FourierTransform. $\endgroup$ – demm Mar 4 '16 at 17:17
  • $\begingroup$ I see what you mean by that, but I don't understand how that agrees with the example of white noise I give in the comment above here. There the factor of 1/pi gives you a completely wrong spectrum, amplitude wise. Also not sure what happens to the factor of two in your comment, but I suppose that's not the crucial part. $\endgroup$ – user129412 Mar 4 '16 at 17:21
  • $\begingroup$ For a white noise process the the S_{xx} is sigma^2/(2Pi) not just sigma^2. Thus since you do not scale properly with the Pi term the Periodogram is around sigma^2=1 as expected. $\endgroup$ – demm Mar 4 '16 at 21:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.