Hi everyone,
it is possible that I already posted a similar question in the past 
, but this one is a slightly different problem. The physical background is the fatigue calculation derived from random vibration with a known PSD. There are many literature sources for that, I will just mention a good practical summary by Tom Irvine (hope the link will work), or just google "Estimating Fatigue Damage from Stress Power Spectral Density Function".
Now my problem is: the input for these approaches is always a continuous integrable function like for example this example from that paper:
Obviously the units of these function are the square of the input quantity divided by Hz (for stress for example ksi**2/Hz or MPa**2/Hz). In this way the square root of the integral of the function will provide the RMS of the input quantity (in this example 6.64 ksi).
Now let's say I have a measured time signal (it is actually a vibration velocity, but it would be no problem to convert it in stress, so that is not the point). I use pspectrum (from the Signal Analysis Extension Pack) to compute its "discrete PSD". Using /cfft**2/ gives similar results, so the computation should be ok.
I say discrete, because while the sum of the "spectral lines" divided by their number is constant (and equal to the sum of the square of the input values, also divided by their number, see illustration from the "spectral analysis chapter of the e-book "Signal Analysis Extension Pack"),
well, now depending on the number of intervals that I chose, I get different "heights". Here an example:
This is obvious, because the sum must remain constant, and so the function, though shown with a track, like it was continuous, is actually made from many vertical "lines" separated one from the other. In other words, the units of this psd output are (I believe) the square of the input quantity, and not, as in the analytical/numerical PSD example above, the quantity divided by Hz.
Now I would like to be able to fit some kind of envelope to the measured PSD, in order to do some computations, something like this:
Possibly the integral of the fitted function should be only slightly above the RMS of the input signal, to avoid overconservativity, but that is another subject.
Now finally my question: how do I get the "right" units [MPa**2/Hz] or [ksi**2/Hz] on the vertical axis? It may be so a dummy question that the answer is: just divide the peaks by their frequency and that's it, but somehow I am unsure if that's all what I need to do.
Thanks a lot for any hints
Best regards
Claudio
ACCEPTED SOLUTION
Accepted Solutions
Hello Claudio Pedrazzi,
you could use the genfit function, as in Help of the Prime and that I present briefly, here below:
Greetings
F. M.
Hi F.M., thanks a lot. I probably wrote too much and so the real question "disappeared" in the flood .
The "real question" is: how do I get from a pspectrum generated spectrum that has values in (for example) [ksi2] to a "normal" function having units of (for example) [ksi2/Hz].
But thanks for suggesting genfit, I will surely try it. But I have to get to the right numerical values before I can even start to fit.
Best regards
Claudio
So I am posting here, as usual, what is apparently working for me. I wished to use the "pspectrum" function because it is very general and has a few interesting input options, missing on cfft. For example the number of intervals, the optional windowing (you know, Hanning, Hamming, Rectangular, Tapered rectangular) and the fractional overlap.
Having said this, pspectrum provides the "power spectrum", per definition [Units: Signal**2]. So "obviously" 
to obtain the "power spectral density" [Units: Signal**2/Hz] it has to be divided by the sampling frequency fs (and multiplied by two if one wants the "one sided" spectrum that makes sense in engineering applications.
psd = 2 * pspectrum / fs
(and use only the first half of the pspectrum output)
With this transformation I obtain a vector of psd points, that if integrated numerically for example with the trapeze integration rule, no matter what the frequency resolution, always correctly reproduces the variance of the original signal. I tested this with different values of n and r in the pspectrum input. With lowering frequency resolution the peaks become flatter and wider, but the integral remains constant, as wished. Here a few examples based on the same signal. Note the ordinate scaling is changing. The integral of the curve is approximately constant and numerically equal to the variance of the signal (sigma**2).
So I believe this should be the right approach. I would be grateful for further hints and of course if someone can point to an error in my reasoning.
Also, thank you again to all answerers, they are highly appreciated.
Best regards
Claudio
PS: LucMeekes I believe that the division by sampling frequency fs is correct, because if it was wrong I would not obtain the correct integral. Dividing by the frequency of each sample would "attenuate" the higher frequencies. Also cfft, after all, has to be squared and divided by fs, in order to produce the psd.
