Power Spectral Density, Welch Technique

[glamotte7]

Hello,

Does anyone know how to calculate a power spectral density from discrete data using Welch's method with Mathematica?

Matlab readily does this but I need to do it with Mathematica.

Any guidance would be appreciated.

glamotte7

 

[Bill Simpson]

Where is your definition of this?

Is it equation 5 in http://www.rssd.esa.int/SP/LISAPATHFINDER/docs/Data_Analysis/PSD_Recipe.pdf
or is it something else?

If you can come up with a precise definition of this then perhaps a couple of lines of code will do it.

 

[glamotte7]

Hello Bill,

I've carefully read the article you suggested which is not entirely, but a bit, over my head. I am new to DSP and so am trying to gain my footing. I'm sure that equation (5) of the paper you referenced would be an excellent approach. I'm not at all sure where to begin to write lines of code to accomplish this as the equation itself is a bit opaque to me. Breaking the data into overlapping segments certainly seems to be a part of the approach in the article, which is what I was expecting. Polynomial curve fitting of segments as well.

Have used the ListConvolve command in Mma and am very impressed with it's flexibility and how clear it makes convolving of the impulse response(kernal) with the data.

Guess that explains where I'm at with this. The whole problem is very intriguing and I suspect there are a number of learning points along the way. Not sure how to start though.

glamotte7

 

[Bill Simpson]

Here is (5) from the paper roughly translated into Mathematica.

2 DeltaT/Nw Sum[Abs[1/Sqrt[Nd]Sum[w[[n]]x[[n]]^l E^((-I 2 Pi k(n-1))/Nd)]^2, {l, 1, Nw}], {n, 1, Nd}]

To use this you need to assign appropriate constant values to DeltaT, Nw, Nd, choose and initialize a windowing function in list w, choose and initialize some plausible time domain signal list x.

Then you should probably try using you data in both Matlab and Mathematica so you can compare the results and we start the process of trying to figure out why the results are completely different, what mistakes have been made, how to fix those, etc, etc, etc.

I would start with really really simple data. A single simple sine wave that exactly fills a single segment with complete cycles is always a nice start. Then the sum of two different sine waves. Then... you get the idea. Gradually ratchet it up with multiple windows and overlapping windows and more complicated time domain data, exterminating every small error before making the next step.

If Matlab anywhere in their documentation describes exactly what they do then we might avoid some of the confusion by trying to port the Matlab method to Mathematica, rather than me just grabbing the first random page that Google showed me and asking you if this was it.

 

[glamotte7]

Bill,

OK. Thanks very much for the code. I'll try working with it as you recommended, very simple things first.

I'll check Matlab to see if they describe their approach elsewhere or what their exact code is. It didn't come up readily but maybe I can find it with more digging around.

glamotte7

 

[milad hoseinp]

hi
I'm looking for a MATLAB code to calculate displacement psd (power spectral density) units[m2/(cycles/m)] vs spatial frequency units[cycles/m]
Could you do me a favor and help me