Thread Subject:

Numerical integration/ calculation of spectral moments

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

"Robbie Brady" <robbie.brady@postgrad.plymouth.ac.uk> wrote in message <jscluv$nfa$1@newscl01ah.mathworks.com>...
> 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
How did you compute frequency3d from frequency?
Barry

"Barry Williams" <barry.r.williamsnospam@saic.com> wrote in message <jscs7i$n5e$1@newscl01ah.mathworks.com>...
> "Robbie Brady" <robbie.brady@postgrad.plymouth.ac.uk> wrote in message <jscluv$nfa$1@newscl01ah.mathworks.com>...
> > 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
> 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

"Robbie Brady" <robbie.brady@postgrad.plymouth.ac.uk> wrote in message <jscutr$6jk$1@newscl01ah.mathworks.com>...
> "Barry Williams" <barry.r.williamsnospam@saic.com> wrote in message <jscs7i$n5e$1@newscl01ah.mathworks.com>...
> > "Robbie Brady" <robbie.brady@postgrad.plymouth.ac.uk> wrote in message <jscluv$nfa$1@newscl01ah.mathworks.com>...
> > > 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
> > 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

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

On Wednesday, June 27, 2012 3:54:07 AM UTC+12, Robbie Brady wrote:
> 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

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.

TideMan <mulgor@gmail.com> wrote in message <1c19ad81-c86a-4216-87fa-ca4cc958caab@googlegroups.com>...
> On Wednesday, June 27, 2012 3:54:07 AM UTC+12, Robbie Brady wrote:
> > 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
>
> 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.

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.