Calculate the area under a section of the curve?

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!

Answer 1

[Star Strider]

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:

.

Answer 2

[Luca Merolla]

Hey, thanks for your reply!

The purpose was exactly to get the area under the curve of just those intervals. Could you help me understand this part of the code please?

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

The purpose was to extract the AUC between three specific intervals, which you stated in the ws matrix. So I guess you solve it! Thank you again very much!

[Star Strider]

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.