Discussion:
Numerical integration/ calculation of spectral moments
(too old to reply)
Robbie Brady
2012-06-26 15:54:07 UTC
Permalink
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
2012-06-26 17:41:06 UTC
Permalink
Post by Robbie Brady
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
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
How did you compute frequency3d from frequency?
Barry
Robbie Brady
2012-06-26 18:27:07 UTC
Permalink
Post by Barry Williams
Post by Robbie Brady
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
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
How did you compute frequency3d from frequency?
Barry
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
2012-06-27 10:49:07 UTC
Permalink
Post by Robbie Brady
Post by Barry Williams
Post by Robbie Brady
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
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
How did you compute frequency3d from frequency?
Barry
"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
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
2012-06-27 20:11:08 UTC
Permalink
Post by Robbie Brady
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
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
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
2012-07-03 09:46:07 UTC
Permalink
Post by TideMan
Post by Robbie Brady
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
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
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?
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.
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
2014-03-09 16:43:08 UTC
Permalink
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
2014-03-09 17:29:08 UTC
Permalink
In other words, am I calculating negative moments in a wrong way?
dpb
2014-03-09 17:55:32 UTC
Permalink
Post by Halvar
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
2014-03-09 18:44:08 UTC
Permalink
Thanks dpb! My first frequency element on the x axis was 0 Hz. I removed that and then it worked fine !
TideMan
2014-03-09 18:44:42 UTC
Permalink
Post by Halvar
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
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
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
2014-03-09 20:53:08 UTC
Permalink
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
2022-06-27 18:39:51 UTC
Permalink
Post by Halvar
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
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
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 ?
Loading...