You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
matlab area calculation under fitted curve
1 view (last 30 days)
Show older comments
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?
Accepted Answer
Rik
on 25 Sep 2020
Edited: Rik
on 28 Sep 2020
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
38 Comments
Ancalagon8
on 25 Sep 2020
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
on 25 Sep 2020
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
on 25 Sep 2020
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
on 25 Sep 2020
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
on 26 Sep 2020
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
Ancalagon8
on 26 Sep 2020
Edited: Ancalagon8
on 26 Sep 2020
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')
Ancalagon8
on 28 Sep 2020
Edited: Ancalagon8
on 15 Oct 2020
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
on 30 Sep 2020
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
on 30 Sep 2020
There is a cooresponding y value for every x value.. So how can i avoid this?
Rik
on 30 Sep 2020
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
on 5 Oct 2020
So these calaculated values what units have? Is it energy?
Ancalagon8
on 5 Oct 2020
And what about energy? I thought it would be equal to the area under the smoothed curve..
John D'Errico
on 5 Oct 2020
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
on 7 Oct 2020
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 , ...)
Ancalagon8
on 7 Oct 2020
Thank you Rik! Any chance of helping me how to syntax the abovementioned?
Rik
on 7 Oct 2020
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
on 9 Oct 2020
y11 only contains real numbers, changed y11(L) with y11(L).^2) but where is the division with total energy?
Rik
on 9 Oct 2020
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
on 11 Oct 2020
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
on 11 Oct 2020
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
on 14 Oct 2020
My problem is that i cannot calculate partial energy with 1 Hz step divided by total energy. (E1/Etotal, E2/Etotal... E45/Etotal)
Ancalagon8
on 14 Oct 2020
Can you please help me implementing the code?
Ancalagon8
on 14 Oct 2020
the step changes. Nothing more.. I can make it again 5 or 10..
Ancalagon8
on 15 Oct 2020
you mean:
A=trapz(x1(L),y11(L));
A=A/sum(A);
Rik
on 15 Oct 2020
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
on 15 Oct 2020
Edited: Ancalagon8
on 15 Oct 2020
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
on 15 Oct 2020
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
on 15 Oct 2020
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:![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/381548/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/381548/image.jpeg)
Ancalagon8
on 15 Oct 2020
Thank you again!
Why you think that the legend is not correct?
More Answers (0)
See Also
Categories
Find more on Spline Postprocessing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)