matlab area calculation under fitted curve

I have a signal and want to calculate the area under the produced fft with step of 5Hz. So far i created an fft, normalized and smoothed it. The problem is that i want to make a loop with decimal step, because the signals' length is 1070 so 1 Hz is 1070/45= 23,77 and step (5Hz) = 118.58

So far my code is:

[x1Data, y1Data] = prepareCurveData(f1, p1norm);
% Set up fittype and options.
ft1 = fittype( 'smoothingspline');
opts1 = fitoptions( 'Method', 'SmoothingSpline');
opts1.SmoothingParam = 0.95; 
% Fit model to data.
[fitresult1, gof1] = fit(x1Data, y1Data, ft1, opts1);
% Plot fit with data.
figure
h1 = plot(fitresult1, x1Data, y1Data,'b' );
h = get(gca, 'Children');

Any ideas how to proceed with the loop?

[Rik]

Comment posted as answer (and erroneously flagged by the spam filter) by Sirius8:

Noone? Come on guys a little help!

Answer 1

[Rik]

You can use trapz for a numerical approximation of an integration. By selecting your bounds you can use this to calculate the area under a curve.

Edit:

Using the code you posted and adding the loop I'm getting something similar to what you describe. If you don't want the legend entries, you should use explicit handles in your call to legend.

S=load('test1.mat');
[x1,y1,y11]=deal(S.x1,S.y1,S.y11);
figure(1),clf(1)

subplot(2,1,1)
plot(x1,y1)
grid on
hold all 
plot(x1,y11,'m')
grid on
hold all 
legend ('normalized fft','smoothed')
subplot(2,1,2)
plot(x1,y11,'m')
grid on
hold all 
legend ('smoothed')

cmap=colormap('prism');n_col=0;
freq_range=5;
for x_start=0:freq_range:max(x1)
    x_end=x_start+freq_range;
    L=x1>=x_start & x1<=x_end;
    if sum(L)<=1,continue,end%skip if there is only 1 point
    X=[x_start x1(L)  x_end];
    Y=[0       y11(L) 0];
    A=trapz(x1(L),y11(L));
    n_col=n_col+1;C=cmap(n_col,:);
    patch(X,Y,C,'DisplayName',sprintf('A(%d-%dHz)=%.2f',x_start,x_end,A))
end

[Ancalagon8]

Rik how to syntax this loop when step is a decimal number? i get an error:

for i=1:23,77:1069

trapz(f1(i), p1norm(i))

end

[Rik]

Decimal values in Matlab follow the American notation: 3/10 is written as 0.3 not as 0,3.

You should be aware that you can't use decimal values as indices, if either of those are arrays.

It also looks like you're inputting scalars in trapz, instead of a range of values.

[Ancalagon8]

I tried also the following but did not worked for me, i got a message of index exceeds array bounds.

step=(length(f1)/45)*5; %so (1070/45)*5=118,58
for i=1:step:(length(f1)-1)
    trapz(f1(i:i+step), p1norm(i:i+step))
end

[Rik]

What are you trying to do? Why did you write your code this way? I have the feeling it wouldn't help if I make the edits to this code to make it run without errors. And how did you create these two variables or functions?

[Ancalagon8]

i have a normalized fft and after i smoothed it I need to calculate the area under the smoothed line within the range 0-45 Hz with step of 5 Hz

[Rik]

What code did you use to plot your lower line chart? If you used variables, can you attach a mat file with the amplitude and frequency?

[Ancalagon8]

the upper graph is normalized fft and the smoothed result, while the lower graph is only the smoothed result:

load test1.mat
subplot(211)
plot(x1,y1)
grid on
hold all 
plot(x1,y11,'m')
grid on
hold all 
legend ('normalized fft','smoothed')
subplot(212)
plot(x1,y11,'m')
grid on
hold all 
legend ('smoothed')

[Rik]

Thank you for posting the data. I will get back to you on this later.

[Ancalagon8]

Rik it works thank you! I have 2 questions:

can i export the legend values in a variable or a matrix?

And when i zoom in it leaves some area out of calculation..

Thanks in advance!

[Rik]

The wedge is there because of how the x-values are selected in the loop. If you want to avoid them, you have to make sure there is a y-value for those x-values.

If you want to store the calculated area: nothing is stopping you. If you read the code you can already see how it is calculated. You can simply store it in a vector.

[Ancalagon8]

There is a cooresponding y value for every x value.. So how can i avoid this?

[Rik]

Yes, but there isn't an x-value of 10 with the corresponding y-value. That means that my code will leave this gap. If you can't resample the curve, you could also round the boundaries to the nearest limit value. Note that this will actually increase the difference between the true value and the calculated value.

%untested coded
for n=1:(max(x1)-freq_range)/freq_range
    x_start=(n-1)*freq_range
    x_end=n*freq_range;
    
    [~,ind]=min(abs(x_start-x1));
    x_start=x1(ind);
    [~,ind]=min(abs(x_end-x1));
    x_end=x1(ind);
    
    %rest of the loop
end

[Ancalagon8]

So these calaculated values what units have? Is it energy?

[Rik]

That question is not about Matlab. The unit would be amplitude/s (i.e. the product of the two axes).

[Ancalagon8]

And what about energy? I thought it would be equal to the area under the smoothed curve..

[John D'Errico]

I would add that a highly osccilatory curve like that is usually poorly integrated using a spline.

This is because the spline is trying to fit noise. Worse, if the spline actually predicts a negative result in some places, you have negative area being summed in there, when it is surely not appropriate.

In the end, a spline, even a smoothing spline, just adds variance to the result.

These are the curves where you want to use trapz - a big reason for the presence of MATLAB. trapz is NOT just another pretty toy, put there for no good reason except for some students to use. In fact, trapz is a minimal variance estimator for an integral on noisy data. So a very useful tool.

If you insist on using the smoothing spline though, there is no reason why you could not have just used fnint, which does come with the curve fitting toolbox and should work on a smoothing spline.

[Ancalagon8]

And what about total energy?

Energy = [sum(abs(x1(1:length(x1))).^2)]

If this is correct, how can i compute the ratio of the energy with 1Hz step/total energy ?

(r1=energy(1Hz:2Hz)/ total energy , r2=energy(2Hz:3Hz)/ total energy , ...)

[Rik]

That looks like you still need to divide by the interval width. If that is the case, you can simply square your data and use tools like trapz to calculate the sum divided by width.

[Ancalagon8]

Thank you Rik! Any chance of helping me how to syntax the abovementioned?

[Rik]

Exactly like what I posted in my answer, but replacing this

A=trapz(x1(L),y11(L));

with this

A=trapz(x1(L),y11(L).^2);

I'm assuming y11 only contains real numbers, otherwise you would still need abs to get the absolute value.

[Ancalagon8]

y11 only contains real numbers, changed y11(L) with y11(L).^2) but where is the division with total energy?

[Rik]

The division by width of the frequency interval (NB: not the total energy) is happening inside trapz, as that is the whole point of using it. You aren't dividing by any total in the formula you posted either, so I assumed that wasn't required. I'm not familiar enough with the mathematics involved in your specific application, so you need to check the math on your own (or with peers). I can help you with implementing it in Matlab.

[Ancalagon8]

Ok understood. For the code below :

for n=1:45
    x_start=(n-1)*freq_range
    x_end=n*freq_range
    
    [~,ind]=min(abs(x_start-x1));
    x_start=x1(ind);
    [~,ind]=min(abs(x_end-x1));
    x_end=x1(ind);
    L=x1>=x_start & x1<=x_end;
%    if sum(L)<=1,continue,end%skip if there is only 1 point
    X=[x_start x1(L)  x_end];
    Y=[0       y11(L) 0];
    ar_ac(n)=trapz(x1(L),y11(L));
    n_col=n_col+1;C=cmap(n_col,:);
    patch(X,Y,C,'DisplayName',sprintf('ar_ac(%d-%dHz)=%.2f',x_start,x_end,ar_ac))
end

how can i change this line in order not only to calculate the area for each of 45 steps

ar_ac(n)=trapz(x1(L),y11(L));

to divide for each of the 45 steps with the total area?

[Rik]

You really shouldn't comment that line, since it is meaningles to calculate the area of a line.

I'm not sure I understand what your problem is. If you want to divide by the total area, what is stopping you? You're storing the partial areas is a vector, so it trivial to divide that vector after the loop by its sum.

[Ancalagon8]

My problem is that i cannot calculate partial energy with 1 Hz step divided by total energy. (E1/Etotal, E2/Etotal... E45/Etotal)

[Rik]

First calculate E(1), E(2), etc, the divide E by sum(E).

[Ancalagon8]

Can you please help me implementing the code?

[Rik]

What happens when you set freq_range=1;?

[Ancalagon8]

the step changes. Nothing more.. I can make it again 5 or 10..

[Rik]

I would have assumed that the area this code calculates also changes. That is what happens in my copy of Matlab. Just store those in an array and devide that array by its sum after filling all elements.

[Ancalagon8]

you mean:

A=trapz(x1(L),y11(L));
A=A/sum(A);

[Rik]

No, because that will overwrite A in every loop iteration. In a previous comment you already showed you understand how to index a variable.

After the loop you can indeed use A/sum(A).

[Ancalagon8]

you mean:

for n=1:45
    x_start=(n-1)*freq_range
    x_end=n*freq_range
        [~,ind]=min(abs(x_start-x1));
    x_start=x1(ind);
    [~,ind]=min(abs(x_end-x1));
    x_end=x1(ind);
    L=x1>=x_start & x1<=x_end;
    if sum(L)<=1,continue,end%skip if there is only 1 point
    X=[x_start x1(L)  x_end];
    Y=[0       y11(L) 0];
    A=trapz(x1(L),y11(L).^2);
    n_col=n_col+1;C=cmap(n_col,:);
    patch(X,Y,C,'DisplayName',sprintf('A(%d-%dHz)=%.2f',x_start,x_end,A))
end
A_new=A/sum(A);

[Rik]

Since you don't index A anywhere in that code: no.

for n=1:45
    x_start=(n-1)*freq_range
    x_end=n*freq_range
        [~,ind]=min(abs(x_start-x1));
    x_start=x1(ind);
    [~,ind]=min(abs(x_end-x1));
    x_end=x1(ind);
    L=x1>=x_start & x1<=x_end;
    if sum(L)<=1,continue,end%skip if there is only 1 point
    X=[x_start x1(L)  x_end];
    Y=[0       y11(L) 0];
    A=trapz(x1(L),y11(L));
    %^^^^^^
    %store this result in a vector and don't square your y here, because then it is not longer actually the area
    A_new(n)=trapz(x1(L),y11(L).^2);
    n_col=n_col+1;C=cmap(n_col,:);
    patch(X,Y,C,'DisplayName',sprintf('A(%d-%dHz)=%.2f',x_start,x_end,A))
end
A_new=A_new/sum(A_new);

[Ancalagon8]

Rik really thank you very much!

i changed:

 A=trapz(x1(L),y11(L));

with:

  A(n)=trapz(x1(L),y11(L));

in order to keep all 45 values from A also. My figure looks like this:

[Rik]

You're welcome.

I doubt your legend is correct now, but that shouldn't be important for the rest of the calculation.

[Ancalagon8]

Thank you again!

Why you think that the legend is not correct?

[Rik]

Because A is a vector once you put it in sprintf to create the legend entry. You should also change the frequency part from '%d' to something like '%.1f'.