Estimate area under a complex curve

I would like to calculate the area under a complex curve whose function is not known.

The graph shows the power spectral density versus frequency.

As part of the design, I need to estimate the total energy which is the integral of the psd vs frequency curve.

I have used trapz (x,y) function to estimate the area under the curve but it seems it is not accurate.

Is there a better way to estimate the area in the matlab?

[KSSV]

Why do you think trapz is not accurate?

Answer 1

[Star Strider]

The trapz function is likely the best option. The ‘accuracy’ may reflect the way you called it. (Without the data and your code, we will never know.) The image shows the power spectral density from about 1490 (Hz?) to about 4000 (Hz?), so it would be necessary to determine the indices of that part of the spectrum and integrate only that part of it.

[Julian M]

Please see the data points.

I need to calculate the whole area under the curve.

I just loaded the data points and used the trapz (x,y) function.

As the the sampling frequency is 57000 Hz and no unit spacing is used, is it correct to use trapz (x,y) function?

[Star Strider]

There is a huge spike at abouit 1892 that may be the reason you suspect the trapz result. It is neverhteless accurate.

To Illustrate —

D = readmatrix('points.xlsx');
f = D(:,1);
psd = D(:,2);

AUC = trapz(f,psd);                                     % Total Area Under The Curve
AUCc = cumtrapz(f,psd);                                 % Cumulative Area Under The Curve

figure
yyaxis left
semilogy(f, psd)
ylabel('PSD(f)')
yyaxis right
plot(f, AUCc)
ylabel('Cumulative Integral Of PSD(f)')
xlabel('Frequency')
grid

Also, the frequency axis begins at 24.2502, not 0, so this also may contribute to what you consider to be an inaccurate result. The trapz output (as illustrated by the cumtrapz plot in yyaxis right) nevertheless appears to be correct. (You can use plot instead of semilogy for yyaxis left, however the spike obscures all the detail of the rest of the PSD.)

.

[John D'Errico]

In fact, when you have a noisy messy curve, the best estimator of the integral is a tool like trapz. You seriously do NOT want to use a higher order approximation on thse problems, because those other tools, anything from say Simpson's rule on up, make assumptions about the smoothness of the curve they are integrating. That is an actively BAD thing to do when the curve is a noisy one. Other tools, for example, a spline fit, followed by integrating the spline are similar in the way they will behave. Avoid them all in the presence of heavy significant noise.

Trapz is the correct tool. If the spacing is a unit one, then only y need be supplied, where stride between points will be assumed to be 1. If the spacing is not 1, then you need to supply the spacing as a scalar if it is a constant. If the spacing is non-uniform, then x needs to be supplied as a vector of the same length as y.

The reason trapz is good here is because the nodal weights are essentially constant across the span. Higher order methods have explicitly (or implicitly) variable point weights across the span of the integration domain. (For an appreciation of this, consider the nodal weights for Simpson's rule, which is about the lowest order rule you can find that is "better" than trapezoidal rule. You will see they vary, bouncing by a factor of 2 back and forth. In fact, that amplifies any noise in the problem, making these higher order rules ones where they can amplify the noise more than they decease the integration error. What happens is higher order rules actually end up costing you more than they gain. I even recall writing this idea up long ago, showing the tradeoff in methods versus the gain.)

You just don't want to use those "better" schemes for integration.

Other people seem to think they can use a tool like integral to integrate such a problem. You cannot do so. Use trapz.

Finally, in case someone suggests you use an interpolation to "smoothly" create more points under the curve, this is also not a good idea. It will not improve the estimate of the integral, because it is implicitly the same thing as using a higher order rule.

So I'm sorry, but there really is no better way to infer the area under a curve that is highly noisy, except to find data that is less noisy. In the presence of noise, the lowest variance estimator will be a tool like trapz. And that is only one of several reasons why trapz is such a valuable tool in MATLAB. It is easy to use, it is efficient, and it is as well-behaved as possible in the presence of noise.

[Star Strider]

John — Thank you!