Calculate the area under a section of the curve?
17 次查看(过去 30 天)
显示 更早的评论
Hello everyone!
I am trying to extract the power of certain frequency bands from my EDA signal. To this purpose, I computed Welch's power spectral density on the signal, and then I wanted to extract the power of the frequencies of interest by calculating the area under the curve between the intervals w = [0, 0.045], w = [0.045, 0.15], and w = [0.15, 0.25]. Does anyone know how I can extract the AUC just within those windows?
I tried to use trapz( ) on just some portions of the signal, but the results look too weird to be accurate. Here is what I did:
% Total signal
EDA_auc = trapz(EDA_pxx); % total power = 0.1063
% AUC for x = [0 0.045]
% cut the signal and do trapz()
y = EDA_pxx(EDA_pxx >= 0 & EDA_pxx <= 0.045);
VLF_auc = trapz(y); % output = 0.1063
% AUC for x = [0.045 0.15]
y1 = EDA_pxx(EDA_pxx >= 0.045 & EDA_pxx <= 0.15);
LF_auc = trapz(y1); % output = 0
The code here should give you a better understanding of the portions of area I would like to calculate:
% w is frequency and is plotted on the x-axis; EDA_pxx is the result of the Welch's PSD analysis and is plotted on
% the y-axis.
figure;
subplot(2,1,1)
plot(w, EDA_pxx)
title("Welch's Power Spectral Density")
xlabel('Normalised frequency')
ylabel('Power')
subplot(2,1,2)
plot(w, EDA_pxx);
title('Portions of interest')
xlabel('Normalised frequency')
ylabel('Power')
hold on
idx = w >= 0 & w <= 0.25;
H = area(w(idx), EDA_pxx(idx));
set(H(1),'FaceColor',[1 0.5 0]);
xlim([0 0.4])
xline(0.045, ':r')
xline(0.15, ':r')
hold off
I hope it's clear what I'm asking for.
Thank you very much in advance for any suggetsion!
0 个评论
采纳的回答
Star Strider
2020-7-8
I am not certain what you want.
Try this:
D1 = load('EDA_pxx.mat');
EDA_pxx = D1.EDA_pxx;
D2 = load('w.mat');
w = D2.w;
ws(1,:) = [0, 0.045];
ws(2,:) = [0.045, 0.15];
ws(3,:) = [0.15, 0.25];
for k = 1:size(ws,1)
lv = w >= ws(k,1) & w <= ws(k,2);
wv{k} = w(lv);
pxxv{k} = EDA_pxx(lv);
AUC{k} = trapz(wv{k}, pxxv{k});
end
pchc = 'rgb';
figure
plot(w, EDA_pxx)
hold on
for k = 1:size(ws,1)
patch([wv{k}; flipud(wv{k})], [pxxv{k}; zeros(size(pxxv{k}))], pchc(k))
end
hold off
grid
fprintf('Area %d = %.3E\n', [1:3; [AUC{:}]])
producing:
Area 1 = 7.149E-05
Area 2 = 8.660E-05
Area 3 = 6.038E-05
and:
.
0 个评论
更多回答(1 个)
Luca Merolla
2020-7-9
1 个评论
Star Strider
2020-7-9
As always, my pleasure!
The ‘ws’ matrix contain the intervals over which you want to calculate the integrals. The ‘lv’ vector is a logical vector that selects the region of ‘w’ defined by ‘ws’ for each section. The ‘wv’ vector are the elements of ‘w’ that correspond to the limits set by ‘ws’, and ‘pxxv’ is the matching vector for ‘EDA_pxx’. The ‘AUC’ vector contains the area for each segment. I kept ‘wv’ and ‘pxxv’ as cell arrays in order to plot them using the patch call.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!