0

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.

enter image description here

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?

enter image description here

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.

enter image description here

Community
  • 1
  • 1
torbonde
  • 2,459
  • 1
  • 15
  • 17
  • Duplicate of http://math.stackexchange.com/questions/468290/how-to-measure-power-spectral-density-in-matlab – Veridian Aug 15 '13 at 15:13
  • Yes. Just as I stated in the first paragraph. – torbonde Aug 15 '13 at 15:18
  • It would not hurt to see: (1) the analytical expression (in particular units) (2) the data. Even better post a minimal working example of your problem. – Buck Thorn Aug 15 '13 at 16:18
  • Seems to me the shortcoming might be with the analytical theory which does not account for the viscosity of air and does not come up with the nice dip? – Buck Thorn Aug 16 '13 at 13:43
  • No, it does. You can see the dip in the analytical PSD as well. Unless you by "dip" means the sort of "upside down peak" (or whatever)? Then it doesn't. But the part of the computed PSD where it grows actually surprised me, since it is not supposed to do that. – torbonde Aug 16 '13 at 14:06
  • Ahh, that is what I meant by dip! – Buck Thorn Aug 16 '13 at 15:00
  • In any case what I am calling a "dip" is probably your biggest problem at this point, I would look for a clear explanation for that... – Buck Thorn Aug 16 '13 at 15:03
  • You probably want to look into function `welch` (not spectrum.welch) and see what it's doing. It could be that dip is due to a filter? Default window appears to be `hamming` – Buck Thorn Aug 16 '13 at 15:31

3 Answers3

1

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.

fpe
  • 2,700
  • 2
  • 23
  • 47
  • Thanks for your response. What would `dt` be in my case? The powerspectrum I get from these two are symmetrical - is that supposed to happen? If I wanted to plot these - say by `loglog(w,Sxx(f))` - what would `w` then need to be? – torbonde Aug 16 '13 at 12:29
  • @cimrg.joe: I suppose that the time series are accompanied by a time-stamp or similar. Am I wrong? Besides, the `PSD` for a real signal is supposed to be symmetric and it's generally called `double-sided PSD`. – fpe Aug 19 '13 at 07:41
0

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.

Buck Thorn
  • 5,024
  • 2
  • 17
  • 27
  • Realized I took the wrong value to compute the scaling factor and changed my answer accordingly. – Buck Thorn Aug 15 '13 at 16:15
  • I have tried passing the sampling rate to `pwelch` instead, but what is produces is actually even further from what I expected. I have updated my question with an explanation of the dip. – torbonde Aug 16 '13 at 12:25
  • Hmmm, but I would argue that inputing the sampling rate is consistent with both the documentation and simple numerical tests of the behavior of pwelch. – Buck Thorn Aug 16 '13 at 12:35
0

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.

rusty
  • 21
  • 1