Numerical integration/ calculation of spectral moments

Robbie Brady

unread,

Jun 26, 2012, 11:54:07 AM6/26/12

to

Hello,

I am currently trying to calculate spectral moment m-1 within matlab by integrating frequency^-1 with respect to power spectral density.

frequency = 64 x 1 x 727 matrix
PSD = 64 x 1 x 727 matrix

I have written the equation as following:

m1 = simpson (frequency3d.^-1.* PSD)

But what is returned is a matrix where each value is infinity, which makes little sense as when done on a calculator I get a numerical value.

If anyone has any idea where I am going wrong I'd be extremely grateful. I should also mention that I have tried other methods of integration (trapz, cumtrapz etc) with no success.

Many thanks,

Robbie

Barry Williams

unread,

Jun 26, 2012, 1:41:06 PM6/26/12

to

"Robbie Brady" <robbie...@postgrad.plymouth.ac.uk> wrote in message <jscluv$nfa$1...@newscl01ah.mathworks.com>...

How did you compute frequency3d from frequency?
Barry

Robbie Brady

unread,

Jun 26, 2012, 2:27:07 PM6/26/12

to

"Barry Williams" <barry.r.wil...@saic.com> wrote in message <jscs7i$n5e$1...@newscl01ah.mathworks.com>...

Sorry, my bad. My original post should read :

"frequency3d = 64 x 1 x 727 matrix"

It is just a 3d matrix of 64 different frequencies, repeated 727 times so that the dimensions would agree with the PSD matrix.

Many thanks,

Robbie

Barry Williams

unread,

Jun 27, 2012, 6:49:07 AM6/27/12

to

"Robbie Brady" <robbie...@postgrad.plymouth.ac.uk> wrote in message <jscutr$6jk$1...@newscl01ah.mathworks.com>...

Without at least a little of your data, it's a bit difficult to see what's going wrong. I'm assuming you are using simpson.m from the File Exchange. By the way, I had to search the FE to find that. It would have been helpful to mention.
When you say that you have run this on a calculator, do you mean that you have run the full set of your data simultaneously, i.e., 64 X 1 X 727?
I'm assuming frequency3d and PSD are double precision. If you are gettin infinite values, obviously something somewhere is exceeding those limits. You could try assigning the result of frequency3d.^-1.* PSD to a variable and examine those values. Do you have a feel for at least order of magnitude of your expected results?
Barry

TideMan

unread,

Jun 27, 2012, 4:11:08 PM6/27/12

to Robbie Brady

I don't know which "simpson" you're using, but don't you need to include either the frequency vector or the interval in frequency if it's equispaced?
Also, you should do:
f=squeeze(f);
PSD=squeeze(PSD);
to remove the surplus dimension.

In contrast to what you think, I think it makes a lot of sense to get infinity as the answer. Have a look at f(1) and think about what that becomes when you invert it.

Robbie Brady

unread,

Jul 3, 2012, 5:46:07 AM7/3/12

to

TideMan <mul...@gmail.com> wrote in message <1c19ad81-c86a-4216...@googlegroups.com>...

My apologies for the late reply. Tideman- your response made me look again at my frequency variable, which I realised was written incorrectly.

Thanks to all who replied for the help.

Halvar

unread,

Mar 9, 2014, 12:43:08 PM3/9/14

to

Hi, I am a little troubled myself. I also would like to calculate the minus first moment but i only get infinite as answer. which is weird, since m-1 / m0 should be the energy period of my spectrum (should be between 1 and 10) of ocean waves, which obviously cannot be infinite. I have no problem calculating m0 with the equation on for example this second page
http://repositorio.lneg.pt/bitstream/10400.9/410/1/1091.pdf
So i tried:
trapz(f,mx.*(f.^-1))
mx is my fft of the surface elevations recorded over time and f is my frequency array (Hz, equally spaced).

Am I typing wrongly or thinking wrongly? m0 is 0.407

Halvar

unread,

Mar 9, 2014, 1:29:08 PM3/9/14

to

In other words, am I calculating negative moments in a wrong way?

dpb

unread,

Mar 9, 2014, 1:55:32 PM3/9/14

to

On 3/9/2014 11:43 AM, Halvar wrote:
> Hi, I am a little troubled myself. I also would like to calculate the
> minus first moment but i only get infinite as answer. which is weird,

> since m-1 / m0 should be the energy period of my spectrum ...
...

Most likely there's a zero in the spectrum you're not eliminating...

--

Halvar

unread,

Mar 9, 2014, 2:44:08 PM3/9/14

to

Thanks dpb! My first frequency element on the x axis was 0 Hz. I removed that and then it worked fine !

TideMan

unread,

Mar 9, 2014, 2:44:42 PM3/9/14

to

Where did you get m-1/m0 as the period?
According to my oceanography textbook it is either:
M0/M1 (mean period) or sqrt(M0/M2) (mean apparent period)
where the nth spectral moment is:
Mn=integral(PSD*f^n)

Halvar

unread,

Mar 9, 2014, 4:53:08 PM3/9/14

to

Hi TideMan. If you want to calculate the wave power, you need the m-1/m0 period. You can see the equation on the link i just posted.

skander Abroug

unread,

Jun 27, 2022, 2:40:01 PM6/27/22

to

Hi,
I'm trying to calculate the 0th and 1st spectral moment order in the case of ocean waves. However, when using the function trapz, i obtain a vector and not a scalar. Am i typing something wrongly ?