How to measure power spectral density in matlab?

Question

I am trying to measure the PSD of a stochastic process in matlab, but I am not sure how to do it. I have posted the exact same question here, but I thought I might have more luck here.

The stochastic process describes wind speed, and is represented by a vector of real numbers. Each entry corresponds to the wind speed in a point in space, measured in m/s. The points are 0.0005 m apart. How do I measure and plot the PSD? Let's call the vector V. My first idea was to use

[p, w] = pwelch(V);
loglog(w,p);

But is this correct? The thing is, that I'm given an analytical expression, which the PSD should (in theory) match. By plotting it with these two lines of code, it looks all wrong. Specifically it looks as though it could need a translation and a scaling. Other than that, the shapes of the two are similar.

UPDATE: The image above actually doesn't depict the PSD obtained by using pwelch on a single vector, but rather the mean of the PSD of 200 vectors, since these vectors stems from numerical simulations. As suggested, I have tried scaling by 2*pi/0.0005. I saw that you can actually give this information to pwelch. So I tried using the code

[p, w] = pwelch(V,[],[],[],2*pi/0.0005);
loglog(w,p);

instead. As seen below, it looks much nicer. It is, however, still not perfect. Is that just something I should expect? Taking the squareroot is not the answer, by the way. But thanks for the suggestion. For one thing, it should follow Kolmogorov's -5/3 law, which it does now (it follows the green line, which has slope -5/3). The function I'm trying to match it with is the Shkarofsky spectral density function, which is the one-dimensional Fourier transform of the Shkarofsky correlation function. Is it not possible to mark up math, here on the site?

UPDATE 2: I have tried using [p, w] = pwelch(V,[],[],[],1/0.0005); as I was suggested. But as you can seem it still doesn't quite match up. It's hard for me to explain exactly what I'm looking for. But what I would like (or, what I expected) is that the dip, of the computed and the analytical PSD happens at the same time, and falls off with the same speed. The data comes from simulations of turbulence. The analytical expression has been fitted to actual measurements of turbulence, wherein this dip is present as well. I'm no expert at all, but as far as I know the dip happens at the small length scales, since the energy is dissipated, due to viscosity of the air.

Answer 1

Presumably you need to rescale the frequency (wavenumber) to units of 1/m. The frequency units from pwelch should be rescaled, since as the documentation explains:

W is the vector of normalized frequencies at which the PSD is estimated. W has units of rad/sample.

Off the cuff my guess is that the scaling factor is

 scale = 1/0.0005/(2*pi);

or 318.3 (m^-1).

As for the intensity, it looks like taking a square root might help. Perhaps your equation reports an intensity, not PSD?

Edit

As you point out, since the analytical and computed PSD have nearly identical slopes they appear to obey similar power laws up to 800 m^-1. I am not sure to what degree you require exponents or offsets to match to be satisfied with a specific model, and I am not familiar with this particular theory.

As for the apparent inconsistency at high wavenumbers, I would point out that you are entering the domain of very small numbers and therefore (1) floating point issues and (2) noise are probably lurking. The very nice looking dip in the computed PSD on the other hand appears very real but I have no explanation for it (maybe your noise is not white?).

You may want to look at this submission at matlab central as it may be useful.

Edit #2

After inspecting documentation of pwelch, it appears that you should pass 1/0.0005 (the sampling rate) and not 2*pi/0.0005. This should not affect the slope but will affect the intercept.

Answer 2

What about using the standard equation for obtaining a PSD? I'd would do this way:

Sxx(f) = (fft(x(t)).*conj(fft(x(t))))*(dt^2);

or

Sxx = fftshift(abs(fft(x(t))))*(dt^2);

Then, if you really need, you may think of applying a windowing criterium, such as

• Hanning
• Hamming
• Welch

which will only somehow filter your PSD.

Answer 3

The dip in PSD in your simulation results looks similar to aliasing artifacts that I have seen in my data when the original data were interpolated with a low-order method. To make this clearer - say my original data was spaced at 0.002m, but in the course of cleaning up missing data, trying to save space, whatever, I linearly interpolated those data onto a 0.005m spacing. The frequency response of linear interpolation is not well-behaved, and will introduce peaks and valleys at the high wavenumber end of your spectrum.

There are different conventions for spectral estimates - whether the wavenumber units are 1/m, or radians/m. Single-sided spectra or double-sided spectra.

help pwelch

shows that the default settings return a one-sided spectrum, i.e. the bin for some frequency ω will include the power density for both +ω and -ω. You should double check that the idealized spectrum to which you are comparing is also a one-sided spectrum. Otherwise, you'll need to half the values of your one-sided spectrum to get values representative of the +ω side of a two-sided spectrum.

I agree with Try Hard that it is the cyclic frequency (generally Hz, or in this case 1/m) which should be specified to pwelch. That said, the returned frequency vector from pwelch is also in those units. Analytical spectral formulae are usually written in terms of angular frequency, so you'll want to be sure that you evaluate it in terms of radians/m, but scale back to 1/m for plotting.