Comment posted as answer (and erroneously flagged by the spam filter) by Sirius8:
Noone? Come on guys a little help!
matlab area calculation under fitted curve
9 次查看（过去 30 天）
显示 更早的评论
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?
1 个评论
采纳的回答
Rik
2020-9-25
编辑：Rik
2020-9-28
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 个评论
Ancalagon8
2020-9-25
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
2020-9-25
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
2020-9-25
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
2020-9-25
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
2020-9-26
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
2020-9-26
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
2020-9-26
编辑：Ancalagon8
2020-9-26
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
2020-9-28
编辑：Ancalagon8
2020-10-15
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
2020-9-30
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.
Rik
2020-9-30
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
Rik
2020-10-5
That question is not about Matlab. The unit would be amplitude/s (i.e. the product of the two axes).
Ancalagon8
2020-10-5
And what about energy? I thought it would be equal to the area under the smoothed curve..
John D'Errico
2020-10-5
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
2020-10-7
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
2020-10-7
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.
Rik
2020-10-7
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
2020-10-9
y11 only contains real numbers, changed y11(L) with y11(L).^2) but where is the division with total energy?
Rik
2020-10-9
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
2020-10-11
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
2020-10-11
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
2020-10-14
My problem is that i cannot calculate partial energy with 1 Hz step divided by total energy. (E1/Etotal, E2/Etotal... E45/Etotal)
Rik
2020-10-14
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.
Rik
2020-10-15
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
2020-10-15
编辑：Ancalagon8
2020-10-15
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
2020-10-15
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
2020-10-15
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
2020-10-15
You're welcome.
I doubt your legend is correct now, but that shouldn't be important for the rest of the calculation.
Rik
2020-10-15
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'.
更多回答（0 个）
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spline Postprocessing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)