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!